关于smos数据的处理--读取，投影转换，裁剪

        smos是目前唯一一个提供不同入射角的被动微波L波段亮温，目前主要用来做海水盐分和土壤水分，也有一些研究用来做地表冻融状态监测，最近老板让我试试smos数据，正好处理的时候遇到一些问题，顺便就记录下来，希望能对别人有所帮助。（主要用的L3数据）

 1. 关于数据

   smos数据是netCDF4格式的，最主要的特点是数据分层存储的，数据纬度为（15，584,1388），即15层，584行，1388列；不同于一般的nc格式数据，smos的数据并不等间距分布，所以Arcgis等一般GIS软件是打不开的，目前来说能打开这个数据的好像只有smosview，整个数据的介绍如下：  

                        

2.投影

       smos L3数据使用的是Ease-grid 2.0投影，关于该投影，具体可见NSIDC http://nsidc.org/data/ease/versions.html。 后期处理需要裁剪，所以这里我试着将该投影转换为常用的wgs1984，具体的思路就是，将4个角的经纬度坐标经投影转换为wgs1984的坐标，然后经度差、纬度差分别处以格网数，这样就可以将整个影像转换为wgs1984投影，且为等间距的，这样就比较方便转换为tiff格式，方便下一步处理。

3.裁剪

       将上面的影像转为tiff影像后，可以根据自己的要求选择tiff影像的波段，例如我需要入射角为52.5°时的TBv和TBh，那我的tiff影像就只有两个波段；最多可以有30个波段，裁剪的方法主要有两种：

（1）利用Arcgis提供的方法，按掩膜裁剪，这个比较方便，因为Arcgis，可以用python批量裁剪，缺点就是数据量大的时候，Arcgis会崩溃，我裁剪了半年左右崩溃了3次；

（2）python的gdal和rasterio包都提供了裁剪的函数，可以参考网上的这些教程的代码

4.部分python代码：

ncfile  = Dataset(FILE_NAME)  #打开nc文件

lat_nc = ncfile.variables['lat'][:]

lon_nc = ncfile.variables['lon'][:]

inc_nc = ncfile.variables['inc'][:]

 

BTh_nc = ncfile.variables['BT_H'][:]

BTv_nc = ncfile.variables['BT_V'][:]

 

lat = np.array(lat_nc)   #[584,]

lon = np.array(lon_nc)   #[1388,]

inc = np.array(inc_nc)   #[15,],入射角

BTh = np.array(BTh_nc)   #[15,584,1388]  

BTv = np.array(BTv_nc)

p1 = Proj(init = 'epsg:6933')   #自己定义的，在epsg文件里面增加了ease-grid 2.0的投影信息

p2 = Proj(init = 'epsg:4326')   #WGS1984

        y3, x3 = p1(np.min(lon_nc), np.max(lat_nc))   #以下是将nc转为tiff

y0, x0 = transform(p1, p2, y3, x3)

tx = np.linspace(y0, -y0, 1388) 

ty = np.linspace(-x0, x0, 584)   #等间隔

basename = os.path.basename(FILE_NAME[19:27])   #文件命名

filename = "{0}.tiff".format(basename)

 

im_data1 = BTv

im_data2 = BTh

im_data = np.append(im_data1,im_data2,axis=0)

print np.shape(im_data)

im_proj  = p2

datatype = gdal.GDT_Float32

im_bands = 30

#获取栅格数据的geotrans

XI, YI = np.meshgrid(tx, ty)

im_geotrans = (XI[0][0], XI[0][1] - XI[0][0], 0, YI[0][0], 0, YI[1][0] - YI[0][0])  #仿射变换参数

driver = gdal.GetDriverByName("GTiff")

dataset = driver.Create(filename, 1388, 584, im_bands, datatype)

dataset.SetGeoTransform(im_geotrans)

 

sr2 = osr.SpatialReference()

sr2.ImportFromEPSG(4326)

dataset.SetProjection(sr2.ExportToWkt())  #设置tiff投影

for i in range(im_bands):

dataset.GetRasterBand(i + 1).WriteArray(im_data[i])

 

5.最终的结果：

          由nc裁剪出的全球的亮温情况，下图是数天合成的

6.问题

 目前SMOS数据用的人非常少，数据的相关处理也基本没有，所以也不确定是否正确，按照官方的说法是smos数据在高原地区质量不太好，下一步想试试L1C的数据，看看L1C在中国的覆盖情况。SMAP数据在中国的覆盖情况很好，就是时间序列不够，只有2015年4月之后的数据，接下来想先看看L1C的数据，然后再试着处理一下SMAP数据，这俩数据比较类似，方法也比较类似。