整个过程分为两步,第一步用python将smap L1B的数据转为csv,实际上就是提取原数据中的经纬度和亮温值,生成一个(lon,lat,data)这样的csv文件,第二步是将csv转tiff,用的r语言。
第一步:hdf转csv(需要自己设置路径)
import numpy as np
import matplotlib.pyplot as plt
import numpy.ma as ma
import glob, os
def printname(name):
print name
import h5py
save_path = " "
Dir_path = " " #path
files = glob.glob(Dir_path+"*.h5")
print files
for file_index in files:
filename = os.path.basename(file_index)[20:41]
print filename
f = h5py.File(file_index,'r')
#f.visit(printname)
#sm_D = f["/Soil_Moisture_Retrieval_Data_AM/soil_moisture"][...]
#lat_D = f["/Soil_Moisture_Retrieval_Data_AM/latitude"][...]
#lon_D = f["/Soil_Moisture_Retrieval_Data_AM/longitude"][...]
sm_A = f["/Brightness_Temperature/tb_h"][...]
sm_B = f["/Brightness_Temperature/tb_v"][...]
sm_c = (sm_B - sm_A)/(sm_B + sm_A)
lat_A = f["/Brightness_Temperature/tb_lat"][...]
lon_A = f["/Brightness_Temperature/tb_lon"][...]
#print sm.shape
#sm = ma.masked_where(sm == -9999, sm)[::].filled(np.nan)
#txt_D = open(save_path_D+filename+'.csv','w')
#txt_D.write("lon,lat,value"+'\n')
txt_A = open(save_path+filename+'.csv','w')
txt_A.write("lon,lat,value"+'\n')
#for lon_content_D,lat_content_D,sm_content_D in zip(lon_D,lat_D,sm_D):
#for lon_index_D,lat_index_D,sm_index_D in zip(lon_content_D,lat_content_D,sm_content_D):
#if sm_index_D > 0.01:
#txt_D.write(str(lon_index_D)+','+str(lat_index_D)+','+str(sm_index_D)+'\n')
for lon_content_A,lat_content_A,sm_content_A in zip(lon_A,lat_A,sm_c):
for lon_index_A,lat_index_A,sm_index_A in zip(lon_content_A,lat_content_A,sm_content_A):
if sm_index_A > 0 and sm_index_A < 300:
txt_A.write(str(lon_index_A)+','+str(lat_index_A)+','+str(sm_index_A)+'\n')
#txt_D.close()
txt_A.close()
f.close()
第二步:csv转tiff(需要自己设置路径和投影)
library("rgdal")
#library("gstat")
library("raster")
library('maptools')
library('sp')
files<-list.files(" ", pattern = '*.csv$',full.names = TRUE)
files_num<-length(files)
save_path<-" "
target_prj <- CRS(" ") #wgs84
#Set the spatial extent of SMAP
extent.geo<-data.frame(x = c(-180,-180,180,180), y=c(86.4,-86.4,86.4,-86.4))
coordinates(extent.geo) <- ~x+y
proj4string(extent.geo)=CRS("+proj=longlat +datum=WGS84 +no_defs ")
extent.prj<-spTransform(extent.geo, target_prj)# Convert latlon extent to spatial projection
#extent.prj<-c(-17334174.8036164, 17340825.1963836, -7344751.36559205, 7355248.63440795)
#SMAP 36km
ncol<-
nrow<-
for (k in 1:files_num){
#Convert csv to point shpefile
filename <- paste(substr(basename(basename(files[k])),1,15),'tif',sep='.')
print (filename)
pts <- read.csv(files[k])
coordinates(pts) <- ~lon+lat
proj4string(pts)=CRS("+proj=longlat +datum=WGS84")
setwd(" ")
if (file.exists("t1.shp"))
file.remove("t1.shp","t1.dbf","t1.prj","t1.shx")
pts<-spTransform(pts, target_prj)
#print (pts)
#setwd('../') #Return last dir
writeOGR(pts, dsn=getwd(),layer="t1",driver="ESRI Shapefile")
writeSpatialShape(pts, "t1")
#Convert point shapefile to asc file
pts<-readOGR('.','t1')
rast <- raster(ncol = ncol, nrow = nrow)#Set the row and col of raster file
extent(rast)<-extent(extent.prj)# Set the spatial extent of raster file
rast2 <- rasterize(pts, rast, pts$value,na.rm=TRUE) # Perform rastersize
#Transform to lon-lat
proj4string(rast2)=target_prj
#rast2<-projectRaster(rast2, crs="+proj=longlat +datum=WGS84")
#writeRaster(rast2, paste(save_path,filename,sep="\\"),overwrite=TRUE) # Write the .asc file
writeRaster(rast2, paste(save_path,filename,sep="\\"), format="GTiff", overwrite=TRUE)
print (paste(save_path,filename,sep="\\"))
}