想发个贴2-记录裁剪dem、在matlab中提取对应数据、回到GIS绘图。

Bubble_liu

想发个贴2-记录裁剪dem、在matlab中提取对应数据、回到GIS绘图。

缘起：我有全省的DEM.tif，500m×500m，及与其配套的栅格气象数据（nc），但是我需要对其中一个市的气象数据做一些处理之后再绘图。因此我决定，从DEM.tif中截出这个市的tif，利用市外的区域为缺省值的方法，在matlab中提取配套的栅格气象数据（全省.nc）中这个市的部分。然后做计算。然后把计算结果套上头行回到GIS去绘成地图。

一、在GIS中裁剪单个市

遇到的问题是，发现可以用clip，也可以用mask。这里记录下对比：

1、用 SpatialAnalyst Tools-> Extraction->Extraction by Mask 工具路径下的mask工具，需要输入三个参数：Input raster选全省DEM.tif，Feature mask data选市.shp，Output raster选市_mask.tif。

file:///C:/Users/admin/AppData/Local/Temp/msohtmlclip1/01/clip_image002.jpg

图1（黑灰为全省DEM.tif，彩色为新栅格）可知，mask出来的栅格，是在原栅格上重新计算的值。因为它的小格格与原来的小格格不是同四个边的小格格。更像是新栅格的顶点取了原栅格的数值。

2、用DataManagement Tools->Raster ->Raster Processing->Clip工具路径下的clip工具，这个工具要分4次来讲，因为，它有2个optional的选项，上面是“Use Input Feature for Clipping Geometry”，下面是“MaintainClipping Extent”。

2.1 两个optional 都不选。

file:///C:/Users/admin/AppData/Local/Temp/msohtmlclip1/01/clip_image004.jpg

图2可知，两个都不选的clip，第一：新栅格与原栅格的四边一致。第二，对比图1、2可知，我的shp肯定不是长方形，但是下面有个长长的直直的边。图3表述更清楚，直接有个直角，表明，两个都不选时，如网上可以查到的，会clip出来一个矩形形状的结果，而不是原shp形状的。

file:///C:/Users/admin/AppData/Local/Temp/msohtmlclip1/01/clip_image005.png

图3

2.2 两个optional ，上选下不选。

file:///C:/Users/admin/AppData/Local/Temp/msohtmlclip1/01/clip_image007.jpg

图4可知，是我想要的结果，就是我自己的shp要再挑个好点的。因为它第一新栅格与原栅格的四边一致，第二clip出来的形状就是我的shp的形状，没有多余的部分。

2.3 两个optional ，上不选下选。

file:///C:/Users/admin/AppData/Local/Temp/msohtmlclip1/01/clip_image009.jpg

图5换了个更能代表的位置截图。第一新栅格与原栅格的四边不一致，第二clip出来的形状是矩形。

2.4 两个optional ，上下都选。

file:///C:/Users/admin/AppData/Local/Temp/msohtmlclip1/01/clip_image010.png

图6可知，第一新栅格与原栅格的四边不一致，第二clip出来的形状是市.shp。

一的总结：

Mask新栅格与原栅格的四边不一致，不是我要的。

Clip有2个optional，上面的选了出来的就是市.shp形状的结果，不选就是大矩形。

下面的选了新栅格与原栅格的四边就不一致，不选就一致。

因为我是要用来在matlab里面选nc的格点，所以我希望他们四边一致，市.shp形状。我用的是clip+上选下不选。在matlab用imread(“绝对路径”)，然后用image画出来看看，好巧。一致。这部分over。

二、在matlab中做计算并回GIS绘图

验证很快乐。我是最后才知道，matlab中可以用geotiffread读取含有投影或地理坐标的tif文件，因此此贴中未有对这个命令做深入探讨。我用了很简单粗暴的方法，哈哈，说起来还挺好笑。

首先，我用imread直接把tif给读了，因缘巧合，还真就是我要那个数据，至于这个怎么验证，简单，在matlab中读到的那个矩阵的对应的点，在gis中直接identify它的信息，你就会看到它的dem高程，看下是不是跟矩阵的数值一样，就可以了。

然后，imread读到的区域对不对？有没有旋转？0点在哪里，如何对应我的nc中的点呢？

这个复杂了。主要是脑袋里面要转过来这个弯。这里用到的重要matlab函数是image。

Image这个函数，直接image(矩阵)，左上角就是原点，关于哪里是原点，建议看官仔细看下横纵坐标。然后我就知道了在我的matlab里面，dem矩阵和nc数据的原点分别在哪里。

备注下：如果image出来的颜色太接近，方法一可以把缺省设为nan，然后把非缺省全部设为500，或者别的你喜欢的数值。Image就只要看轮廓就好了嘛，数值大小不重要。方法二：image(矩阵,'CDataMapping','scaled')和colorbar，可以把颜色框在数值的范围上，看起来会方便点。

备注下二：我检查了我imread出来的dem数值，它在非本市的区域里，定义为缺省。

备注下三：为了脑子少转几个弯，可以在matlab里面使用rot90函数，把dem和nc的原点转到同一个方向。

第三步：自己手绘一下，看看，你的dem的原点，在nc数据中对应行列值，然后找到你的dem区域，在nc数据中的行列值。写个for循坏，把需要的数据提出来。提出来的这个数据就是你要的dem区域的nc数值了。

这里要特别注意的是，matlab没有（0，0）这个位置，它只有（1，1），所以脑子要转过来，（1，1）代表的是哪个格点，（end，end）代表的又是哪个格点。他们分别在dem和nc里面又是哪个格点，从dem到nc，行列分别是加上多少才对，这里要多检查几遍。以免偏移了。

第四步：会不会已经忘了自己把nc数据提出来是要干嘛的？求均值？求最大？求初日终日？爱干嘛干嘛。处理到自己最后想要的那个数据矩阵。需要在地图上画出来的那个矩阵之后，我们进入第五步。

第五步：把矩阵（不含投影、地理坐标信息，但是又与dem长相一致的矩阵），在GIS中画出来。

哭一下，这一步让我痛苦了很久。

我提供一个方案，不代表仅此一个方案：

1、在GIS中，toolbox里面有个ConversionTools—From Raster—Raster to ASCII。用它把dem.tif转换为txt，用笔记本打开找到头行，也就6行左右吧。

2、在matlab中把你第四步的计算结果写到txt里，把dem.tif的txt的头行粘贴到这个txt里。

3、写好txt后，回到GIS，toolbox里面有个ConversionTools—To Raster—ASCII to Raster。

备注1：matlab写的txt，间隔符号要跟GIS里面导出来的txt一致。

备注2：如果后续还有很多要画图的，可以写个matlab程序把头行自动写道你的txt前面，因为头行有string，可以用用fprintf，一个一个写进去。先写头行，然后用dlmwrite('S1.txt',S1,'-append','delimiter',' ')来写。

备注3：第四步的矩阵可能跟dem.tif的txt的矩阵的方向不一样，如果在GIS里画出来的是横条横条，不跟原dem一样，那你就在matlab里面rot90，直到旋转到一样为止，最多也就试4次必出正确角度，对吧。

————————————————

此处罗列参考网址和文章

【1】       https://blog.csdn.net/weixin_39591632/article/details/109691330 ArcGIS中mask(掩膜提取)和clip工具的一点思考

【2】       https://blog.csdn.net/Cheese_pop/article/details/71192114matlab cell类型数组存至txt文件

【3】       https://ww2.mathworks.cn/help/index.html强大的matlab帮助中心。最好用没有之一。