ERA5-land中total precipitation、2mtemperature、Surface solar radiation downwards的下载和月度数据处理

蓦亓

于 2024-08-10 17:04:53 发布

阅读量447

点赞数 10

分类专栏：遥感数据下载及处理文章标签：算法 python

本文链接：https://blog.csdn.net/weixin_56184060/article/details/141092901

版权

遥感数据下载及处理专栏收录该内容

2 篇文章 0 订阅

订阅专栏

一、数据产品介绍

ERA5-Land monthly averaged data from 1950 to present** 是一个由欧洲中期天气预报中心（ECMWF）通过Copernicus Climate Change Service (C3S) 提供的再分析数据集。这个数据集提供了自1950年至今的月平均陆地表面气象数据，广泛用于气候研究、环境监测和决策支持。

1. 数据集概述：

数据格式：nc
开发机构: 欧洲中期天气预报中心（ECMWF）。
数据覆盖范围: 全球陆地表面（不包括海洋）。
时间跨度: 从1950年1月至今，且数据不断更新。
空间分辨率: 9公里（0.1°）。
时间分辨率: 月平均数据，每个变量的月平均值是基于该月内的小时数据进行平均计算得出的。

2. 主要变量:

该数据集提供了一系列与陆地表面相关的关键气象和气候变量，这些变量被按月平均值提供，包括但不限于以下内容：
2米气温: 近地表温度（月平均、日最高、日最低）。
降水量: 月降水总量。
地表蒸发: 总蒸发量。
土壤湿度: 多个深度层的月平均土壤湿度。
风速和风向: 10米高度的月平均风速和风向。
地表净辐射: 月平均短波、长波净辐射。
雪深: 月平均雪深。
地表压力: 月平均地表气压。

3. 数据集特点:

长时间跨度: 数据集从1950年起覆盖至今，使其成为研究长期气候变化的有力工具。
高空间分辨率: 9公里的空间分辨率比典型的再分析数据要细致，适用于局部和区域尺度的分析。
高时间分辨率: 基于小时数据计算的月平均值提供了细致的时间尺度信息，适合分析季节性和年度变化。

4. 数据应用:

气候变化研究: 用于分析长期气候变化趋势，如温度、降水和极端天气事件的变化。
农业和水资源管理: 提供与土壤湿度、降水、蒸发等相关的数据，支持农业管理和水资源规划。
-生态系统研究: 支持对生态系统健康状况的长期监测，如分析植被的季节性变化、干旱和降水对生态系统的影响。
气候模型验证：用于验证和校准气候模型，确保模型预测的准确性。

5. 数据获取:

Copernicus Climate Data Store (CDS): 数据集可以通过Copernicus Climate Data Store免费下载和访问。用户可以通过CDS API编程接口或直接在在线平台上选择时间段和变量来下载数据。

https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=form

二、python下载和处理代码

1.nc转tif

#以获取气温为例
# -*- coding: utf-8 -*-
import netCDF4 as nc
from osgeo import gdal, osr
import numpy as np
import os
#1.2定义写图像文件的函数
def write_img(filename, im_proj, im_geotrans, im_data):
    # 判断栅格数据的数据类型
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32
    # 判读数组维数
    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1, im_data.shape
    # 创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
    dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    dataset.SetProjection(im_proj)  # 写入投影
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset

def nc_totif(input_path, output_path):
    # 读取nc文件
    tep_data = nc.Dataset(input_path)
    # 获取nc文件中对应变量的信息
    lon_data = tep_data.variables["longitude"][:]
    lat_data = tep_data.variables["latitude"][:]
    # 影像的左上角&右下角坐标
    lonmin, latmax, lonmax, latmin = [lon_data.min(), lat_data.max(), lon_data.max(), lat_data.min()]
    # 分辨率计算
    num_lon = len(lon_data)  # 269
    num_lat = len(lat_data)  # 157
    print(num_lon)
    print(num_lat)
    lon_res = (lonmax - lonmin) / (float(num_lon) - 1)
    lat_res = (latmax - latmin) / (float(num_lat) - 1)
    # 定义投影
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(4326)  # WGS84
    proj = proj.ExportToWkt()  # 重点，转成wkt格式
    # print(prj)     字符串
    geotransform = (lonmin, lon_res, 0.0, latmax, 0.0, -lat_res)
    # 获取2m温度
    t2m = tep_data.variables["t2m"][:]  #降水：tp 气温：t2m 辐射：ssrd
    t2m_arr = np.asarray(t2m)
    #年份
    yearlist = [2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020]#下载数据的年份
    for i in range(len(yearlist)):
        year = yearlist[i]
        for j in range(12*i, 12*(i+1)):
            month = (j % 12) + 1
            outputpath = output_path + os.sep + str(year) + "_" + str(month) + "_tm.tif"
            write_img(outputpath, proj, geotransform, t2m_arr[j])
if __name__ == "__main__":
    # nc文件输入输出路径
    input_path = r"D:\test\data.nc"
    output_path = r"D:\result"
    # 读取nc文件，转换为tif文件
    nc_totif(input_path, output_path)

2.投影+重采样+裁剪研究区

import arcpy
import os
from arcpy.sa import ExtractByMask
Coordinate_System = "PROJCS['WGS_1984_Albers',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',105.0],PARAMETER['Standard_Parallel_1',25.0],PARAMETER['Standard_Parallel_2',47.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0],AUTHORITY['EPSG',32651]]"
inputdir=r"D:\test"
outFolder = r"D:\result"
mask = r"D:\data\mask.tif" # 参考栅格
arcpy.env.parallelProcessingFactor = '0'
arcpy.env.workspace = inputdir
rasterList = arcpy.ListRasters("*", "TIF")
for raster in rasterList:
    outname = outFolder + os.sep + raster[5:12] + "_proj.tif"
    if os.path.exists(outname):
        print(outname + "have exist")
    else:
        arcpy.ProjectRaster_management(raster, outname, Coordinate_System, "BILINEAR",  "#", "#", "#", "#") # 设置投影
        print (outname + "have been done")
        outExtractByMask = ExtractByMask(outname, mask) #裁剪图像
        outExtractByMask.save(outFolder + os.sep + raster)

3.单位换算+月度数据处理

降水：

# -*- coding: UTF-8 -*-
import arcpy
arcpy.CheckOutExtension("Spatial")
import os
output_path = r"D:\tpmonth"
# 降水根据每个月天数求月总降水量
inputpath = r"D:\tp"
arcpy.env.workspace = inputpath
rasters = arcpy.ListRasters("*tp*","TIF")
print(rasters)
for raster in rasters:
    outFolder = output_path + os.sep + str(raster[0:4])
    isExists = os.path.exists(outFolder)
    if not isExists:
        os.makedirs(outFolder)
    else:
        outname = output_path + os.sep + str(raster[0:4]) + os.sep + raster
        if raster[5:7] in ['1_', '3_', '5_', '7_', '8_', '10', '12']:
            out = -31000*arcpy.Raster(raster)
            out.save(outname)
        if raster[5:7] in ['4_','6_','9_','11']:
            out = -30000*arcpy.Raster(raster)
            out.save(outname)
        if raster[5:7] in ['2_']:
            if int(raster[0:4])/4 == 0:
                out = -29000 * arcpy.Raster(raster)
                out.save(outname)
            else:
                out = -28000 * arcpy.Raster(raster)
                out.save(outname)
    print(outname + "完成")

气温：

# -*- coding: UTF-8 -*-
import arcpy
arcpy.CheckOutExtension("Spatial")
import os
output_path = r"D:\tmmonth"
inputpath = r"D:\tm"
arcpy.env.workspace = inputpath
rasters = arcpy.ListRasters("*tm*","TIF")
print(rasters)
for raster in rasters:
    outname = output_path + os.sep + str(raster[0:4]) + os.sep + raster
    out = arcpy.Raster(raster) - 273.15  # 温度单位换算
    out.save(outname)
    print(outname + "完成")

辐射：

# -*- coding: UTF-8 -*-
import arcpy

arcpy.CheckOutExtension("Spatial")
import os

output_path = r"D:\srdmonth"
inputpath = r"D:\srd"
arcpy.env.workspace = inputpath
rasters = arcpy.ListRasters("*rd*", "TIF")
print(rasters)
for raster in rasters:
    outname = output_path + os.sep + str(raster[0:4]) + os.sep + raster
    out = arcpy.Raster(raster) / 86400   # 辐射单位换算
    out.save(outname)
    print(outname + "完成")

蓦亓

关注

10
点赞
踩
9

收藏

觉得还不错? 一键收藏
0
评论
ERA5-land中total precipitation、2mtemperature、Surface solar radiation downwards的下载和月度数据处理

数据格式：nc开发机构: 欧洲中期天气预报中心（ECMWF）。数据覆盖范围: 全球陆地表面（不包括海洋）。时间跨度: 从1950年1月至今，且数据不断更新。空间分辨率: 9公里（0.1°）。时间分辨率: 月平均数据，每个变量的月平均值是基于该月内的小时数据进行平均计算得出的。
复制链接

扫一扫