Google Earth Engine谷歌地球引擎GEE地理坐标系与投影坐标系变换与重投影

  本文主要对GEE中地理坐标系与投影坐标系转换、重投影等操作加以介绍。本文是谷歌地球引擎（Google Earth Engine，GEE）系列教学文章的第十三篇，更多GEE文章请参考专栏：GEE学习与应用（https://blog.csdn.net/zhebushibiaoshifu/category_11081040.html）。

  在第十二篇GEE教学博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/119644114）中，我们介绍了GEE中图层投影信息（包括基准参考坐标系及其空间转换参数）的获取方法；本文则在这一基础之上，进一步介绍GEE中图层的重投影操作。

  在第九篇GEE教学博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/119545059）中，我们介绍了用代码导入Landsat系列遥感影像的方式；为了更为全面地展示GEE的应用，且突出本文即将介绍的投影转换操作，本文则改用MODIS系列遥感影像中的第006版MOD13A1全球16日植被指数Vegetation Indices 16-Day Global 500m产品为例，进行后续操作的介绍。

var modis=ee.Image("MODIS/006/MOD13A1/2020_06_09");
print(modis);
Map.addLayer(modis);
Map.setCenter(116.36863,39.961029,9);

  执行代码，可以看到遥感影像已经出现在交互式地图中，且地图已经自动缩放到北京市新街口外大街附近区域。当然，在这里有一点需要注意：一般做植被指数相关研究，肯定很少会将研究区域的中心设立在城市内，更多是对森林、农田进行探究；而本文由于重点在于讲解GEE的相关操作而非真正进行具体研究，因此就直接沿用了之前教学博客中的区域。

  接下来，依据第十二篇GEE教学博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/119644114）中内容，打印这一景MODIS遥感影像的投影信息。

print(modis.projection());

  执行代码，得到结果如下。

  其中，得到的SR-ORG:6974是许多MODIS遥感影像产品使用的一种全球正弦投影，是投影坐标系的一种。

  随后，依据第四篇GEE教学博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/117573181）中内容，在代码中进行图层可视化配置，并重新将遥感影像导入地图中；这样使得遥感影像更加直观，方便我们后续的操作。

var vis_parameter={
  min:2039.2255686140938,
  max:8445.630153788503,
  bands:["NDVI"],
  palette:["63ff14","41f6ff","5f37ff","e033ff","ff3f7d"]
};
Map.addLayer(modis,vis_parameter,"NDVI");

  执行代码，得到结果如下。

  如果大家在ArcMap等软件中处理过MODIS遥感影像数据，会发现若不对其自身投影信息加以修改，其像元就会呈现出一种“扭曲”的倾斜状态；而我们将GEE地图放大，同样可以看到，在GEE中导入MODIS遥感影像亦会出现这样的情况。

  接下来，我们就将这一景遥感影像转换为墨卡托投影（Mercator projection），从而消除上述像元倾斜的状态。在GEE中，坐标系变换的代码其实非常简单：

var new_modis=modis.reproject("EPSG:4326");
print(new_modis);

  其中，.reproject()函数用以变换指定图层的坐标系，"EPSG:4326"参数就是我们常用的World Geodetic System（WGS）1984地理坐标系，是地理坐标系的一种。这里需要注意，在GEE中进行投影变换操作时，其默认通过最邻近插值方法进行重采样，从而实现新图层的生成；如果不希望用这一默认的重采样方法，可以通过.resample()函数或.reduceResolution()函数实现，具体请见第十四篇GEE教学博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/119721332）。

  执行代码，结合右侧所显示的两景遥感影像的信息，可以看到投影转换对于影像的基本信息似乎并没有很大影响。

  但是，当我们打印投影转换后遥感影像的投影信息，并与投影转换前图层的投影信息进行对比，便可以看到原有遥感影像的SR-ORG:6974投影坐标系已经变为了"EPSG:4326"地理坐标系，说明坐标变换已经完成。

print("New projection information is:",new_modis.projection());

  执行代码，得到结果如下。

  随后，我们将坐标变换后的新的图层导入地图中显示，查看其效果。

Map.addLayer(new_modis,vis_parameter,"NDVI_NEW");

  执行代码，得到结果如下。

  可以看到，坐标变换后的图层尽管像元已经不再是倾斜的状态，但是其空间分辨率发生显著变化，降低了很多。依据第十二篇GEE教学博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/119644114）中介绍的.projection().nominalScale()函数，查看新图层的空间分辨率。

print("New resolution is:",new_modis.projection().nominalScale());

  执行代码，得到结果如下。

  可以看到，坐标变换后的图层空间分辨率已经达到111319m，即111 km左右；而我们所使用的是MOD13A1全球16日植被指数Vegetation Indices 16-Day Global 500m产品，其原本的空间分辨率是500 m左右。由此可知，经过了坐标变换，遥感影像的空间分辨率也随之发生变化。

  出现这个问题的原因在于，我们在利用.reproject()函数进行坐标变换时，没有指定变换后图层的空间分辨率，那么系统将自动以转换后的坐标系默认空间分辨率作为所得新图层的空间分辨率；为了保证图层的空间分辨率在进行坐标变换时不发生改变，我们可以在执行.reproject()函数时将新图层的期望空间分辨率作为参数放入函数内部。

var new_modis=modis.reproject({crs:"EPSG:4326",scale:500});
//print(new_modis);
print("New projection information is:",new_modis.projection());
Map.addLayer(new_modis,vis_parameter,"NDVI_NEW");
print("New resolution is:",new_modis.projection().nominalScale());

  其中，我们在.reproject()函数中增添了对scale参数的赋值，从而保证新图层的空间分辨率为500 m。这里还用到了第十篇GEE教学博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/119568274）中介绍的将函数的参数修改为字典的形式的方法，从而跳过.reproject()函数中一些不需要的参数。

  执行代码，得到结果如下。

  可以看到，经过上述操作，所得到的坐标变换后图层空间分辨率为500 m，恢复了正常的情况。

欢迎关注公众号/CSDN/知乎/微博：疯狂学习GIS