使用MATLAB 将EASE-Grid 2.0投影坐标系下的NC文件转换为相同坐标系下的geotiff文件

最新推荐文章于 2024-08-20 11:03:44 发布

路飞的橡皮拳

最新推荐文章于 2024-08-20 11:03:44 发布

阅读量1w

点赞数 19

文章标签：科研日记

本文链接：https://blog.csdn.net/weixin_38953602/article/details/101158084

版权

以SMOS L3 土壤水分产品数据为例

SMOS L3 级土壤水分产品是 .nc 格式的，分辨率为25km，要将其转换为geotiff格式的，这里有两种方法。
第一种方法是最简单的，在网上找一个EASE-Grid2.0 25km的 .tif 模板数据，将nc文件里的数据读进来，按照对应经纬度和对应行列号，直接替换掉原来模板中的数据。

楼主因为没有找到这个模板，自己摸索出了一种方法。

EASE-Grid2.0 不同分辨率网格参数

首先弄清楚每种分辨率下EASE-Grid 2.0 的行列号数，和起始 x,y 坐标，这个网址给出了详细说明：
https://nsidc.org/ease/ease-grid-projection-gt 网址链接
截图如下：
在这里插入图片描述

使用maprefcells创建对应网格参照系

为此使用以下代码创建一个参考坐标系

%%EASE-Grid 2.0  25 km 投影参数，x y坐标范围
xworldlimits=[-17367530.45,17367530.45];
yworldlimits=[-7307375.92,7307375.92];
% EASE-Grid 2.0  25km 行列号数
rastersize=[584,1388];
% 创建参考坐标系
R= maprefcells(xworldlimits,yworldlimits,rastersize,'ColumnsStartFrom','south');

R的属性信息如下图：
在这里插入图片描述
这里相当于替大小为584*1388的矩阵创建了一个参考坐标系，而且属性表中的Columns相当于对应着纬度，Rows 相当于对应着经度，参数 ‘ColumnsStartFrom’ 要设置为‘south’的原因是因为下载的SMOS nc 文件中，记录的数据，纬度是从南往北排列的，相当于上下颠倒了，使用 Panoply软件查看的数据记录方式如下：

在这里插入图片描述
再读取nc文件中的数据，以土壤湿度为例，

sm=ncread(InFile,'Soil_Moisture');

然后使用geotiffwrite函数写入tiff文件

info=geotiffinfo('D:\EASE_Grid2.0\EASE_Grid2_9km.tif');
geotiffwrite('output.tif',sm',R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
%% ncread读进来的sm数据需要转置

这里需要注意的是 geotiffwrite 函数中，第一个参数为输出文件名，第二个参数为输出数据矩阵，楼主将NaN值全部替换为了0，以便后续再Arcgis中显示；第三个参数为EASE-Grid2.0 25km 对应的参考坐标系，前面已经设置好。最后一个参数 ‘GeoKeyDirectoryTag’特别特别重要，特别重要，特别重要，重要的事情说三遍，geotiffwrite本身是针对地理参考坐标系（即常用的WGS-84坐标）所设计，而对于像EASE-Grid2.0 这样的平面投影坐标系，输出为tif文件时，‘GeoKeyDirectoryTag’ 这个参数的值一定要设置好。
来自MATLAB帮助文档中geotiffwrite函数里的一段介绍
楼主 ‘GeoKeyDirectoryTag’ 这个变量取值来自于一个EASE-Grid2.0 9km的tiff模板，可能你会疑问为什么有EASE-Grid2.0 9km的tiff模板了，怎么不用arcgis直接重采样为25km的，这里也需要强调的是EASE-Grid2.0 1km, 3km, 9km, 36km分辨率是一个体系，它们的x和y起始坐标是一样的；而3.125km，6.25km, 12.5km, 25km的又是一个体系，它们的x和y起始坐标是一样的，但与9km的不一样，使用arcgis直接将9km重采样为25km时，会多出一行。前面挂的链接有对不同分辨率下EASE-Grid2.0 起始x,y坐标详细的介绍，9km截图如下，和前述对比，y坐标的起始值大一点。这也正是楼主为什么要绕了一大圈要这么辅以程序来生成25km网格的tif文件的原因。而关于9km的网格，可以在earthdata 上下载SMAP数据时（SMAP数据也为EASE-Grid2.0），选择‘Customize’,将输出文件格式转为Geotiff，就可以获取9km 的模板tiff文件。
在这里插入图片描述
‘GeoKeyDirectoryTag’ 变量是EASE-Grid2.0投影坐标系的一种固有属性，与分辨率无关，楼主对比了1km和9km下这个变量的值，都是一样的，因而认定25km的也一样。

转换后的tiff文件检查

将转换好的tif文件用arcgis打开，如下图：

在这里插入图片描述
使用arcgis 查看属性信息，如下：

虽然行列号数对应上了，但是Cell size有一点小小的问题，长宽尺寸差个0.00001，严格上说应该都是相等的。

注： EASE-Grid2.0全称是Equal-Area Scalable Earth Grid 2.0，依据适用范围的不同分为3种子投影：Global，Northern Hemisphere，Southern Hemisphere。三者都是基于WGS84椭球体，但投影方式有所差异，其中，“EASE-Grid2.0 Global”属于Cylindrical_Equal_Area（等积圆柱投影），是一种切投影，标准纬线为赤道；“EASE-Grid2.0 Northern Hemisphere”和“EASE-Grid2.0 Southern Hemisphere”采用Lambert_Azimuthal_Equal_Area（兰伯特等积方位投影），切点分别为南北极点。

以上纯属科研小白的学习日记。