干货分享 | 基于Python和ArcGIS的空间插值和空间Mann-Kendall(M-K)分析_mask

相关文章推荐

重感情的熊猫 · 关于公布北京协和医学院2024博士研究生复试 ...· 2 月前 ·

面冷心慈的酱牛肉 · 美第奇家族︱出了两位法国王后，把艺术和时尚带入巴黎· 4 月前 ·

儒雅的领结 · 大庆市委书记李世峰到肇州县、肇源县调研-黑龙 ...· 7 月前 ·

绅士的炒饭 · 我国各省产油排名去年我曾经统计过我国的各省 ...· 11 月前 ·

果断的风衣 · 魂兵之戈（魂兵之戈漫画版）_百度百科· 1 年前 ·

原标题：干货分享 | 基于Python和ArcGIS的空间插值和空间Mann-Kendall(M-K)分析

Mann-Kendall分析 在地学研究中是很常见的一种研究手段，我们可以在很多场合看到这个方法的身影。 ArcGIS 提供的 空间数据插值 之前也学习过了，是可以通过站点，也就是具体的经纬点来插值的。那么接下来就先做ArcGIS对CRU降水格点数据做提取和插值美化工作，好好地出一张地图出来，然后使用Python对下载的CRU降水格点数据做Mann-Kendall趋势分析。

然后打开ArcMap的多维工具箱从nc创建一个栅格图层，维度选择time维。

加载完成之后，顺便把shp也加进去，然后在数据管理工具箱中合成全部band，由于我做的时候已经全部做好了，这感觉似乎不是很好写，我也只能将一个大概的步骤了。把数据提取到point去，就可以在shp上得到数据了，统计均值，如果你的数据原本就被方大了十倍的话，那么可以使用地图代数去除以10。。。。我发现这部分真不好写，特别是这个图，一大长条的太影响美观了。

最后，通过提取到point的格点，也就是shp范围内的格点来做ID插值，需要说明的是从CRU下载的数据是全球尺度的，于是提取到point去之后我就clip出了shp青藏高原的来了，没有mask，因为这里两个数据都已经是矢量了。

对提取到point的降水数据进行IDW插值之后再进行简单的美化修饰操作就得到了如下的结果图了，然后打印出地图来看看

这样的话，一张简单的地图就完成了，这个测试的目的就是提取栅格点出去插值。

第二部分就是使用Python就是对下载来的CRU降水数据做空间Mann-Kendall测试分析。为了让测试看起来比较理想，我选择下载另一个CRU文件，时间范围比较长，1901-2020年的，前面的网址打开之后可以去下载到。测试环境为Jupyter Notebook。需要的一个关键包是pymannkenall，之外还需要xarray，numpy等。

这个nc数据大小为2.78G，我也不想选择一个子区域了，直接就做全球的

如果还没安装pymannkendall的话，先安装与喜爱，然后导入必要的包并读取，这部分代码零零散散的我就以截图形式展示

trends = ['decreasing','no trend','increasing'] # 使用布尔索引把str变到float去 outdatac[outdatac==trends[2]]=1 # up outdatac[outdatac==trends[0]]=-1 # down outdatac[outdatac==trends[1]]=0 # 没啥趋势的格点

# outdatac = np.int_(outdatac) # tren = xr.DataArray(outdatac, dims=('lat','lon'), coords={'lat':ext_data.lat, 'lon':ext_data.lon}) # cs = tren.where(tren!=np.nan).plot.contourf

由于数据有点大，计算要花点时间的，算完之后就可以使用cartopy去绘制地图了

import cartopy.crs as ccrs import cmaps import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

lon, lat = ext_data.lon, ext_data.lat nx, ny = np.meshgrid(lon, lat)

plt.figure(figsize=(16, 8))

#让colorbar字体设置为新罗马字符 plt.rcParams['font.family'] = 'Times New Roman' plt.rcParams['font.size'] = 20

ax = plt.subplot(projection=ccrs.PlateCarree) ax.coastlines(lw=0.4,color='k') ax.set_global

cb = ax.contourf(nx,ny,outdatac, cmap=cmaps.GMT_jet, levels=[-2,-1,0,1],transform=ccrs.PlateCarree) cbar = plt.colorbar(cb,shrink=0.6,pad=0.01)

cbar.set_ticks([-1.5, -0.5,0.5]) cbar.set_ticklabels(['decreasing', 'no trend', 'increasing'])

plt.title('CRU 1901-2020 yr Precipitation Mann-Kendall Test',pad=15, fontdict={'family' : 'Times New Roman', 'size' : 20}) plt.grid ax.set_xticks(np.arange(-180, 180, 30),crs=ccrs.PlateCarree) plt.xticks(fontproperties = 'Times New Roman',size=16) plt.yticks(fontproperties = 'Times New Roman',size=16) ax.set_yticks(np.arange(-90, 90, 30),crs=ccrs.PlateCarree) ax.xaxis.set_major_formatter(LongitudeFormatter)#经度0度不加东西 ax.yaxis.set_major_formatter(LatitudeFormatter)

然后就得到了结果图：

最开始的时候，colorbar并不是有文字在上面的，需要使用如下脚本添加一下：

cbar.set_ticks([-1.5, -0.5,0.5]) cbar.set_ticklabels(['decreasing', 'no trend', 'increasing'])

转载自 Python干货铺子

经作者授权转载

文章仅代表作者观点，与本公众号无关，版权归原作者所有

原文标题：基于Python和ArcGIS的空间插值和空间Mann-Kendall(M-K)分析

排版：许多

审编：黄莘绒

终审：颜子明黄宗财鲁嘉颐