网站:Copernicus Climate Data Store | Copernicus Climate Data Store

数据目标:2013-2020年温度日数据

一、数据下载

在网站上注册后,点击Datasets,以“ERA5-Land”搜索,选择“ERA5-Land hourly data from 1950 to present

1. 手动下载每小时数据

点击Download data,选择需要的时间和地理范围,点击Submit Form,可以在Your Requests中查看请求进程,数据量越大等待时间越久。

2. python实现自动下载

(1)cdsapi库的安装

cdsapi库的安装可以移步官方教程How to install and use CDS API on Windows - Copernicus Knowledge Base - ECMWF Confluence Wiki

如果安装了Pycharm,可以点击File->Settings->Python Interpreter,点击➕加号,搜索【cdsapi】点击安装即可

(2)保存.cdsapirc文件

点击这个链接 How to use the CDS API | Copernicus Climate Data Store

复制下图右方的url和Key,新建文本文件,放入,再修改文件拓展名为【.cdsapirc】

 

 将该文件放入C盘->用户->用户名,文件夹内。

(3)python批量下载

代码主体可从官网上获取,下方是实现循环批量下载月数据

# coding=utf-8
import cdsapi

c = cdsapi.Client()
for year in range(2013, 2021):
    for month in range(1, 13):
        outpath = 'E:\\test\\' + str(year) + str(month).zfill(2) + '.nc'
        print(outpath)
        c.retrieve(
            'reanalysis-era5-land',
            {
                'variable': '2m_temperature',
                'year': str(year),
                'month': [
                    str(month).zfill(2)
                ],
                'day': [
                    '01', '02', '03',
                    '04', '05', '06',
                    '07', '08', '09',
                    '10', '11', '12',
                    '13', '14', '15',
                    '16', '17', '18',
                    '19', '20', '21',
                    '22', '23', '24',
                    '25', '26', '27',
                    '28', '29', '30',
                    '31',
                ],
                'time': [
                    '00:00', '01:00', '02:00',
                    '03:00', '04:00', '05:00',
                    '06:00', '07:00', '08:00',
                    '09:00', '10:00', '11:00',
                    '12:00', '13:00', '14:00',
                    '15:00', '16:00', '17:00',
                    '18:00', '19:00', '20:00',
                    '21:00', '22:00', '23:00',
                ],
                'area': [
                    60, 70, 0,
                    140,
                ],
                'format': 'netcdf',
            },
            outpath)

批量下载参考:

https://blog.csdn.net/qq_34734252/article/details/108781538

http://bbs.06climate.com/forum.php?mod=viewthread&tid=100616&highlight=ERA5

3. 手动下载日数据

这里用到官网上自带的应用:

Daily statistics calculated from ERA5 data

不过每次只能计算一个月的

二、数据处理

1. 每小时数据取平均得到每日数据

# !usr/bin/env python
# -*- coding: utf-8 -*-

import os, sys
import netCDF4 as nc
import numpy as np
import arcpy as ap
from arcpy.sa import *
import datetime
import gdal,osr
def dayextract_nc(nc_file, output_dir):
    dataset = nc.Dataset(nc_file)
    Lat = dataset.variables["latitude"][:]
    Lon = dataset.variables["longitude"][:]

    # 影像的左上角和右下角坐标
    LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()]

    # 分辨率计算
    N_Lat = len(Lat)
    N_Lon = len(Lon)
    Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1)
    cols = N_Lon
    rows = N_Lat

    # 设置影像的显示范围
    geo = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)

    # 构造projection
    src_srs = osr.SpatialReference()
    src_srs.ImportFromEPSG(4326)  # 定义输出的坐标系为"WGS 84",AUTHORITY["EPSG","4326"]
    src_srs_wkt = src_srs.ExportToWkt()   # 给新建图层赋予投影信息

    st = datetime.datetime(1900, 1, 1, 0, 0)
    time1 = []

    for n in time:
        time1.append(str(st + datetime.timedelta(hours=int(n))))

    key = np.array(dataset.variables['t2m'][0, :, :])

    ## 测试是否需要缩放和偏移,结果表明不需要
    # print(dataset.variables['t2m'].set_auto_maskandscale(True))
    # print(dataset.variables['t2m'][0, :, :])  # data1
    # print(dataset.variables['t2m'].set_auto_maskandscale(False))  # 设置缩放和掩膜数组关闭
    # print(dataset.variables['t2m'][0, :, :])  # 关闭后打印数据,得出的数据全部为整数
    # 将关闭后的数据应用缩放和偏移转换
    # print(dataset.variables['t2m'][0, :, :] * dataset.variables['t2m'].scale_factor + dataset.variables['t2m'].add_offset)  # data2

    dayarray = np.zeros(key.shape, key.dtype)
    dayarray1 = np.zeros(key.shape, key.dtype)

    hour = 0
    a = 0
    for i in range(0, len(time)):
        dayarray = dayarray + np.array(dataset.variables['t2m'][i, :, :])
        hour = hour + 1
        if hour % 24 == 0:
            numarray = np.ones(dayarray.shape, dayarray.dtype) * 24
            dayarray = dayarray / numarray
            out_file = os.path.join(output_dir, time1[a][0:10] + ".tif")
            print(out_file)
            write_tiff(out_file, geo, src_srs_wkt, rows, cols, dayarray)
            a = a + 24
            dayarray = dayarray1

这里有一个点,我用matlab查看Nc数据时发现,有scale_factor和offset_factor,就在想这个数据是否是真实值还是需要处理的,于是测试了一下是否需要偏移和缩放。

一般来说如果值是整数有可能需要,如果是小数则不需要

测试是否需要缩放和偏移参考:

https://blog.csdn.net/Will_Zhan/article/details/105388609

2. 转tif

def nc_to_tif(nc_file, output_dir):
    nc_file_name = os.path.basename(nc_file)

    # 提取nc文件
    dayextract_nc(nc_file, output_dir)

def write_tiff(tiff_file, geo, proj, rows, cols, data_array):

    driver = gdal.GetDriverByName("Gtiff")

    outdataset = driver.Create(tiff_file, cols, rows, 1, gdal.GDT_Float32)
    outdataset.SetGeoTransform(geo)
    outdataset.SetProjection(proj)

    band = outdataset.GetRasterBand(1)
    band.WriteArray(data_array)

Logo

为开发者提供学习成长、分享交流、生态实践、资源工具等服务,帮助开发者快速成长。

更多推荐