其他分享
首页 > 其他分享> > ERA5数据处理,简单记录。Arcpy

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