SMAPL1B转tiff详细代码

        整个过程分为两步，第一步用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="\\"))

}