ERA5数据处理，简单记录。Arcpy
作者：互联网
第一次写这个，还有点小激动，新手，简单记录一下，地学同一问题方法上可能有不同，代码和方法如果有改进的，还请多多指出！ hhh
其实代码还在跑着，不过目前没问题，先写上，写的是ERA5数据下载完经过mosaic后的处理过程，下次有时间再写下载和拼接的过程吧。
以提取中国2015年小时级的lai_lv（Leaf area index, low vegetation）变量为例（ERA5的数据提取的时候是按变量的简写来提取的，详见ERA5-Land: data documentation - Copernicus Knowledge Base - ECMWF Confluence Wiki）
思路也包含在代码里，直接上代码咯。
说一下思路哈。拼接完的tif面积比中国大。ERA5hourly数据0.1°左右，也就是10km左右（大概，别喷我）。目标是：得到中国1km的数据
涉及到投影，重采样，裁剪啥的（这个先后随便写的）
首先，由于是陆地数据，海洋无值，10km的格子和1km的格子对于海岸线边界来说差别挺大的。所以一开始重采样不行。
1.先设空值，把异常值，本文中是<0的设为0。然后你会发现，海岸线有些格子为0。所以中国面积小了，那怎么行，得给它补上hhhh。
2.咋办呢，用邻域分析呗。 气象要素总体来说是平稳的。于是我们用附近的值给它补上。我用的是半径为2个unit的圆。总体效果还可以。这样你会发现，中国的面积大了，因为只要是圆内有一个格子有值，它就会做邻域分析。等会儿裁剪重采样就行了。
3.记住了，邻域分析完是另一张平滑过的图，和设置了空值的那张图是两张图，我们要的是空值那张图中缺失的海岸线的部分格子就行，也就是这些格子的值从邻域分析的图取。千万别把邻域分析的图直接用到后面，因为这样内陆地区的值都不是原始数据。 这个步骤就是很长的那句，fill_raster = 。。。。。。。。。。。。。。。。
4.投影，就是调用一下，方法自己选哈。
5.重采样到一公里，方法自己选哈。
6. 用中国边界掩膜，格子对应上。
建议在原数据上先做设空值，邻域， 先做重采样啥的沿海边界会出问题，因为海上原本是没有值的，重采样还填充值。。。。。 处理顺序要注意一下
下面是代码
# -*- coding: utf-8 -*-
#导入包
import arcpy
from arcpy import env
from arcpy.sa import *
import os,datetime
#获取程序开始时间
start = datetime.datetime.now()
#设置处理环境
variable = "2015_lai_lv" #这个就是个要处理文件夹的名称而已，文件夹下有一堆要处理的数据
final_Dir = "D:/第一次文件夹/ERA5/" +variable #这个是最终的结果输出文件夹，就是可以用的数据。一般不用中文路径，这个路径因为有姓名，随便打的哈。
if not os.path.exists(final_Dir):
os.makedirs(final_Dir) #就是个普普通通的判断语句
arcpy.env.workspace = final_Dir #设置工作空间，还挺重要的！！
arcpy.env.parallelProcessingFactor = "21" #调用电脑的多核运行，现在很多电脑都是多核的，但是有没有用呢，官方给出的是部分gis工具有并行处理的功能
elevRaster = 'D:/随便文件夹/Boundary.tif' # 这个是个裁剪的tiff，就是中国边界。 至于路径嘛，写文章时随便取的名字，跑代码还是英文好一点。
#打开工具箱，获得License
arcpy.CheckOutExtension("Spatial")
#覆盖输出
arcpy.env.overwriteOutput = True
#获取投影
Projection = arcpy.Describe(elevRaster).spatialReference
#打开文件
# variables = os.listdir(r"D:\Data_collect\ERA5\ERA5_tiff")
# for variable in variables: #获取输入文件夹文件名，当时为了循环很多个变量写的，最终数据量太大太慢了，放弃了。复制代码的时候自己删除吧。
inDir = "F:/Data_collect/ERA5/ERA5_tiff/" +variable #打开文件夹，需要处理数据的文件夹
# prj_Dir = "D:/Data_collect/ERA5/try/ERA5_fill_prj/" +variable
# resam_Dir = "D:/Data_collect/ERA5/try/ERA5_fill_prj_resam/"+variable
# if not os.path.exists(prj_Dir):
# os.makedirs(prj_Dir)
# if not os.path.exists(resam_Dir):
# os.makedirs(resam_Dir) #这些注释嘛，没用。因为arcpy有些spatial analysis 工具用的时候需要保存一个位置，所以生成了一些中间数据的文件夹。
inList = [name for name in os.listdir(inDir) if name.endswith(".tif")] # 一个简单的循环，获取所有tif文件
# 遍历所有tif
for file in inList:
# 输入文件的完整路径
mem_Dir = "in_memory/" #这是个内存路径
input_raster = os.path.join(inDir, file) #join获取完整路径
final_raster = os.path.join(final_Dir,file)
#判断文件是否存在，若存在则跳过
if os.path.exists(final_raster):
continue
set_null_ra = SetNull(input_raster, input_raster, "Value < 0") #设置小于O的值为空值
fill_raster = arcpy.sa.Con(arcpy.sa.IsNull(set_null_ra),
arcpy.sa.FocalStatistics(set_null_ra, arcpy.sa.NbrCircle(2, "cell"), "MEAN", "DATA"),
set_null_ra) # 邻域计算并替代空值
file = file[:-4] #获取文件的文件名但不包括扩展（即.tif），因为中间数据放在内存空间时是不允许有文件拓展的
prj_file = "prj_"+file #中间数据名
resam_file = "resam" +file #中间数据名
prj_raster = os.path.join(mem_Dir,prj_file)
resam_raster = os.path.join(mem_Dir,resam_file) #join 同上，不过是内存的路径， 要不然直接命名“in_memory/.......”应该也可以
arcpy.ProjectRaster_management(fill_raster, prj_raster, Projection, "NEAREST", ) #重投影
arcpy.Resample_management(prj_raster,resam_raster, "1000", "BILINEAR") # 执行重采样
arcpy.gp.ExtractByMask_sa(resam_raster, elevRaster,final_raster) # 执行掩膜，注意输出文件名长度
arcpy.management.Delete(prj_raster) #删除内存的文件，怕内存不够用
arcpy.management.Delete(resam_raster)
#获取程序的结束时间
end = datetime.datetime.now()
SpendTime = end-start
print ("%s is ok!"%variable)
print("spend time :%s"%SpendTime)
我也不知道这玩意咋修改，也太难用了。代码段传上去就改不了了。 经常看到保姆级啥的，这怎么也得小学级了。代码可能有点乱，刚学习，见谅哈。
几个问题解释一下：
1.因为我有21个变量，当时本来想写个循环，遍历所有，太慢了。而且硬盘空间啥的不够。后面我就用一个变量一个python文件来跑了，反正就是改variable = 某个变量嘛，改也挺快。
2.设置处理环境时，为什么把创建输出文件夹啥的语句放前面呢。1.中说的那样， 是要注意的是arcpy执行的时候，多个python的workspace如果一样，因为下面调用arcpy.某个工具的时候会生成一些临时文件，用完了会自动删除。但是如果workspace都一样，有可能会因为其他python文件在用这个文件而删不掉，而且这些临时文件还没办法命名，反正就是可能出现啥unable to delete .....啥的，然后程序就gg了，我也不知道是不是因为这个，但是我把多个python文件的workspace改了之后就不出错了。 代码就是上面那样就行
3.关于调用多核进行啥的，arcpy.env.parallelProcessingFactor = "21"就是这句，目前官方说的是有的工具可以加快，21就是你要用的核数，但我这个程序没起到啥效果，想用查查，或者找其他人博客看看，反正就是企图加速用的，官方说明Parallel Processing Factor (Environment setting)—ArcGIS Pro | Documentation
4.我虽然写了arcpy.env.overwriteOutput = True 输出覆盖，但是前面跑过的很多输出的文件，已经输出了，代码出错了，要再跑一遍可太揪心了。ERA5用过都知道，一年8760个小时，要是跑了几千个，再跑一遍没必要。所以先加个判断语句，判断一下文件在不在if os.path.exists(final_raster):
continue
至于arcpy.某个工具啥的语法应用，复制一下，搜搜就出来了。我就不多说了。
说一下我最开心的一点就是 将中间数据存到内存空间，然后再删除就行。速度直接加了好几倍。最开始的时候arcpy.ProjectRaster_management、arcpy.Resample_management、arcpy.gp.ExtractByMask_sa这些工具都要求你指定输入和输出路径，就出现了我傻傻的创建多个中间数据的文件夹。先输出到，再读入tif给下一步用，慢的........ 后面将中间数据输出到内存，后面步骤调用完就删除，快了好几倍，hhhh，改代码的乐趣。 一定要记得删除，这个内存是有限的。也就是delete的那个！！！！！
代码写的有点乱，新手，多见谅。然后很多注释都是写这边文章的时候写上去的。还有些是之前尝试的，删了就行。
最后放张图，研究区是中国大陆地区哈，当然最后研究作完出图还得注意南海九段线啥的，国家领土完整嘛。
有点啰嗦，第一次写，别嫌弃哈。。。。。。
标签：raster,os,Arcpy,prj,ERA5,数据处理,arcpy,Dir 来源： https://blog.csdn.net/weixin_43927082/article/details/121379301