本文主要是介绍Arcpy / Matlab / Arcgis处理CMIP6数据,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
目前找到了直接使用matlab处理的方法,关键是使用matlab的interp2或griddata函数,正在试验,目前结果还可以,待后续更新
本文主要是对下载的CMIP6数据进行重采样及区域裁剪。
对于CMIP6数据,看它的网格格式是什么样的,如果是gn格式,那就代表该数据不是等间距投影,不能直接用Arcgis的创建NetCDF栅格格式会报错,也不能直接用matlab提取,会存在投影误差(下文我会提到,如果不想看可以直接跳过)。如果是等间距投影可以直接用matlab编程提取(我感觉matlab提取比较方便。。下图是用matlab提取研究范围的思路,但不能用于gn格式的数据提取)
ncdisp([ncpath,ncname])
lon = ncread([ncpath,ncname],'lon');
lat = ncread([ncpath,ncname],'lat');
time = ncread([ncpath,ncname],'time');
[lat1,~] = find(lat >= 1.9 & lat <= 24.1);
lat_min = min(lat1); lat_max = max(lat1); % 列的范围
[lon1,~] = find(lon >= 105 & lon <= 121);
lon_min = min(lon1); lon_max = max(lon1); % 行的范围
start_loc = [lon_min,lat_min,1];
count = [numel(lon1),numel(lat1),numel(time)];
sst_data = ncread([ncpath,ncname],'sst',start_loc,count);
sst_3D = rot90(sst_data);
文章目录
- gn格式的CMIP6数据存在的问题
- arcgis处理CMIP6数据
- arcpy批量处理
处理思路:Arcgis提供了创建NetCDF要素图层(Make NetCDF Feature Layer)的工具,这个工具可以处理gn格式的CMIP6数据。随后再用要素转栅格得到研究区域的数据。
为什么要用Arcpy呢?因为我的CMIP6数据是多时间段的(基本都是多时间),在Arcgis里做很麻烦我想批量化处理。
gn格式的CMIP6数据存在的问题
利用matlab读取下载的nc文件,查看一下这个数据的变量、经纬度等信息
clear all; clc
[ncname,ncpath] = uigetfile('*.*','请选择nc文件');
ncdisp([ncpath,ncname])
lat = ncread([ncpath,ncname],'latitude');
lon = ncread([ncpath,ncname],'longitude');
我下载的数据是SSP126情景下2100年1月1日起逐日全球海表温度数据,空间分辨率为25km(实际有投影误差,可以看到经纬度都是一个二维矩阵)
打开经纬度矩阵(lat和lon实际需要逆时针旋转90度才是正常的地理分布,但我这里没变所以lat看每列,lon看每行),可以发现,同一纬线下经度变形不大,而同一经线下纬度变化很大。如果我直接提取研究区范围(南海区域),可以看一下下图纬度的变形程度
可以看到同一经线下纬度已经不是0.25°了,在逐渐缩小,所以这样提取出来的区域要比实际0.25°等间距区域大。
arcgis处理CMIP6数据
我先尝试着用arcgis处理了一下这个数据,Make NetCDF Feature Layer,具体参数设置:
variables只选择我要的变量海温tos
row dimensions选择i和j
dimension value选择time
随后会出现巨多的点点,每个点代表了该位置的海温,空白区域即无数据。
这个数据默认显示的是第一天的,可以从该图层属性里修改你想看到的日期(但是感觉日期不对啊我是2100年的数据这里显示2096???,但是经过对比实际上还是从2100年1月1日开始的。研究了一下是因为arcgis不会读取nc文件time的日历calendar,后面arcpy批量操作那里可以看到我对时间time做了操作,times_true会计算出正确的时间)
随后我要对这个要素图层转栅格,使用Conversion Tools里的Feature to Raster
output cell size设置为0.25(也就是~25km)
注意这里可以直接设置我想要的区域,但是设置的前提是你已经有一个裁剪好的区域栅格数据了。
也就是说,第一次用要素转栅格,先做全区域的转,然后用clip裁剪出你想要的区域,后续再用要素转栅格,你就可以在Environments里Processing Extent设置研究区域,就不用再clip裁剪。
最后就能得到我研究区域的海温数据了
arcpy批量处理
其实上面arcgis做那么多都是为了在arcpy里做批量处理打基础,因为acrpy那些函数的参数设置就是我们刚才设置的那些参数
强烈建议使用BY INDEX
import arcpy# Execute NetCDFtoFeatureLayer
arcpy.env.workspace = "F:/0TCLI/30_Years_Daily_SST/CMIP6/ncresample"
inNetCDFFile = "F:/0TCLI/30_Years_Daily_SST/CMIP6/tos_SSP_data/tos_Oday_HadGEM3-GC31-MM_ssp126_r1i1p1f3_gn_21000101-21001230.nc"
inVariables = "tos"
inXVariable = "longitude"
inYVariable = "latitude"
outFeatureLayer = "tos_layer1"
rowDimensions = ["i", "j"]
ZVariable = ""
MVariable = ""
dimensionValues = "time"
valueSelectionMethod = "BY_VALUE"# Execute MakeNetCDFFeatureLayer
arcpy.MakeNetCDFFeatureLayer_md(inNetCDFFile, inVariables, inXVariable,inYVariable, outFeatureLayer, rowDimensions,ZVariable, MVariable, dimensionValues,valueSelectionMethod)# Set local variables
inFeature = outFeatureLayer
outRaster = "F:/0TCLI/30_Years_Daily_SST/CMIP6/ncresample/SCS_tos_1_arcpy.tif"
cellSize = 0.25
field = "tos"
# Set Extent Range: Right Top Left Bottom
arcpy.env.extent = "121.125000 24.137497 104.875000 1.887497"# Execute FeatureToRaster
arcpy.FeatureToRaster_conversion(inFeature, field, outRaster, cellSize)
用arcpy提取的和arcgis逐步操作结果一致
这里只做了一个时间的,其实可以用循环提取所有时间,等后续更新吧
2022-07-19更新:想不到啥好方法直接在python里计算平均值,主要是这个数据属是无语它是gn格式,然后arcpy这些函数生成的栅格图层只能导出不能放在工作空间里当作变量(或者是我不知道栅格变量在哪),所以我只能把每个月的数据给导出然后再在matlab里处理吧。。。。
import arcpy
import netCDF4 as nc# Calculate NetCDF data time
ncfilepath = "E:/2022_Summer/Marine_Environment_Data/CMIP6_data/"
ncfilename = "tos_Omon_GFDL-CM4_ssp585_r1i1p1f1_gn_203501-205412.nc"
ncfileid = ncfilepath + ncfilename
ncfile_obj = nc.Dataset(ncfileid)
time = ncfile_obj.variables['time']
times_true = nc.num2date(time[:], time.units, time.calendar)# Execute NetCDFtoFeatureLayer
for i in range(len(time)):outpath = "E:/2022_Summer/Marine_Environment_Data/CMIP6_data/NOAA_CMIP6_tos/"arcpy.env.workspace = ncfilepathinVariables = "tos"inXVariable = "lon" # ATTENTION HEREinYVariable = "lat" # ATTENTION HEREoutFeatureLayer = "tos_layer"rowDimensions = ["x", "y"] # ATTENTION HEREZVariable = ""MVariable = ""dimensionValues = "time" + " " + str(i)valueSelectionMethod = "BY_INDEX"# Execute MakeNetCDF Feature Layerarcpy.MakeNetCDFFeatureLayer_md(ncfileid, inVariables, inXVariable,inYVariable, outFeatureLayer, rowDimensions,ZVariable, MVariable, dimensionValues,valueSelectionMethod)# Set local variablesinFeature = outFeatureLayerfiletime = times_true[i].strftime('%Y-%m-%d')outRaster = outpath + filetime + ".tif"cellSize = 0.25field = "tos"# Set Extent Rangereftif = r"E:\2022_Summer\Marine_Environment_Data\CMEMS_data\CHL_L4_test.tif"arcpy.env.extent = reftif# Execute FeatureToRasterarcpy.env.pyramid = "NONE" # No Pyramidarcpy.FeatureToRaster_conversion(inFeature, field, outRaster, cellSize)print("Date: %s No.%d / Total: %d" % (filetime, i, len(time)))arcpy.Delete_management(outFeatureLayer)
运行的时候可以看到还剩多少数据
这篇关于Arcpy / Matlab / Arcgis处理CMIP6数据的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!