Newer
Older
Jeremy Auclair
committed
# -*- coding: UTF-8 -*-
# Python
Jeremy Auclair
committed
04-07-2023
Jeremy Auclair
committed
@author: rivallandv, heavily modified by jeremy auclair
Jeremy Auclair
committed
Jeremy Auclair
committed
Download ERA5 daily weather files for modspa
Jeremy Auclair
committed
"""
import sys # for path management
import os # for path exploration
import geopandas as gpd # to manage shapefiles
from psutil import cpu_count # to get number of physical cores available
from p_tqdm import p_map # for multiprocessing with progress bars
Jeremy Auclair
committed
import modspa_pixel.preprocessing.lib_era5_land_pixel as era5land # custom built functions for ERA5-Land data download
from modspa_pixel.config.config import config # to load modspa config file
Jeremy Auclair
committed
from modspa_pixel.preprocessing.parcel_to_pixel import convert_dataframe_to_xarray
Jeremy Auclair
committed
Jeremy Auclair
committed
def request_ER5_weather(config_file: str, ndvi_path: str, raw_S2_image_ref: str = None, shapefile: str = None, mode: str = 'pixel') -> str:
Jeremy Auclair
committed
Download ERA5 reanalysis daily weather files, concatenate and calculate ET0
to obtain a netCDF4 dataset for precipitation and ET0 values. Weather data
reprojection and conversion can take some time for large spatial windows.
Jeremy Auclair
committed
Arguments
=========
1. config_file: ``str``
Jeremy Auclair
committed
json configuration file
Jeremy Auclair
committed
2. ndvi_path: ``str``
Jeremy Auclair
committed
path to ndvi cube, used for weather data reprojection
Jeremy Auclair
committed
3. raw_S2_image_ref: ``str`` ``default = None``
unmodified sentinel-2 image at correct resolution for
weather data reprojection in pixel mode
4. shapefile: ``str`` ``default = None``
path to shapefile for extraction in parcel mode
5. mode: ``str`` ``default = 'pixel'``
choose between ``'pixel'`` and ``'parcel'`` mode
Jeremy Auclair
committed
Returns
=======
1. weather_file: ``str``
path to netCDF4 file containing weather data
Jeremy Auclair
committed
"""
Jeremy Auclair
committed
# Get config file
Jeremy Auclair
committed
config_params = config(config_file)

Jérémy AUCLAIR
committed
# Get parameters
data_path = config_params.data_path
shapefile_path = config_params.shapefile_path
run_name = config_params.run_name
mode = config_params.mode
start_date = config_params.start_date
end_date = config_params.end_date
weather_overwrite = config_params.weather_overwrite
# Set output paths

Jérémy AUCLAIR
committed
outpath = os.path.join(data_path, 'WEATHER', run_name)

Jérémy AUCLAIR
committed
# Create daily wheather ncfile name

Jérémy AUCLAIR
committed
weather_daily_ncFile = os.path.join(outpath, start_date + '_' + end_date + '_' + run_name + '_era5-land-daily-meteo')

Jérémy AUCLAIR
committed
# Check if file exists and ndvi overwrite is false
if os.path.exists(weather_daily_ncFile + '.nc' * (mode == 'pixel') + '_parcel.nc' * (mode == 'parcel' )) and not weather_overwrite:
return weather_daily_ncFile
Jeremy Auclair
committed
# Geometry configuration
wgs84_epsg = 'epsg:4326' # WGS84 is the ERA5 epsg
# ERA5 product parameters
grid_size = 0.25
Jeremy Auclair
committed
wind_height = 10 # height of ERA5 wind measurements in meters
print('REQUEST CONFIGURATION INFORMATIONS:')

Jérémy AUCLAIR
committed
if shapefile_path:
if os.path.exists(shapefile_path):
print('shapeFile: ', shapefile_path)
Jeremy Auclair
committed
else:
print('shapeFile not found')
else:
# print('specify either shapeFile, boxbound or point coordinate in json file')
print('specify shapeFile in json file')
sys.exit(-1)

Jérémy AUCLAIR
committed
print('period: ', start_date, ' - ', end_date)
print('experiment name:', run_name)
Jeremy Auclair
committed
if os.path.exists(outpath):
print('path for nc files: ', outpath)
else:
os.mkdir(outpath)
print('mkdir path for nc files: ', outpath)
print('----------')
Jeremy Auclair
committed
# Request ERA5-land BoxBound Determination

Jérémy AUCLAIR
committed
if shapefile_path:
Jeremy Auclair
committed
# Load shapefile to access geometrics informations for ERA5-Land request

Jérémy AUCLAIR
committed
gdf_expe_polygons = gpd.read_file(shapefile_path)
Jeremy Auclair
committed
print('Input polygons CRS :', gdf_expe_polygons.crs)
expe_epsg = gdf_expe_polygons.crs
# verification que les polygones sont tous fermés
liste_polygons_validity = gdf_expe_polygons.geometry.is_valid
if list(liste_polygons_validity).count(False) > 0:
Jeremy Auclair
committed
Jeremy Auclair
committed
print('some polygons of Shapefile are not valid')
polygons_invalid = liste_polygons_validity.loc[liste_polygons_validity == False]
print('invalid polygons:', polygons_invalid)
Jeremy Auclair
committed
Jeremy Auclair
committed
for i in polygons_invalid.index:
gdf_expe_polygons.geometry[i]
# Application d'un buffer de zero m
gdf_expe_polygons_clean = gdf_expe_polygons.geometry.buffer(0)
gdf_expe_polygons = gdf_expe_polygons_clean
# search for the total extent of the whole polygons in lat/lon [xlo/ylo/xhi/yhi] [W S E N]
expe_polygons_boxbound = gdf_expe_polygons.total_bounds
Jeremy Auclair
committed
expe_polygons_boxbound = list(expe_polygons_boxbound)
print('shape extend in ', expe_epsg.srs, ':', expe_polygons_boxbound)
if expe_epsg.srs != wgs84_epsg:
print('--- convert extend in wgs84 coordinates ---')
Jeremy Auclair
committed
Jeremy Auclair
committed
# idem en wgs84 pour des lat/lon en degree (format utilisé par google earth engine)
expe_polygons_boxbound_wgs84 = gdf_expe_polygons.to_crs(wgs84_epsg).total_bounds
Jeremy Auclair
committed
else:
expe_polygons_boxbound_wgs84 = expe_polygons_boxbound
print('boxbound [N W S E] extend in ', wgs84_epsg)
print(expe_polygons_boxbound_wgs84)
Jeremy Auclair
committed
print('--start request--')
# Get number of available CPUs
nb_processes = 2 * min([cpu_count(logical = False), len(os.sched_getaffinity(0)), config_params.max_cpu]) # downloading data demands very little computing power, each processor core can manage multiple downloads
Jeremy Auclair
committed
# Define variables
variables = ['2m_temperature', '2m_dewpoint_temperature', 'surface_solar_radiation_downwards', '10m_u_component_of_wind', '10m_v_component_of_wind', 'total_precipitation']
# Split download per year if necessary
dates = era5land.split_dates_by_year(start_date, end_date)
if len(dates) == 1:
# Generate arguments
args = [(var, outpath, start_date, end_date, expe_polygons_boxbound_wgs84, grid_size) for var in variables]
else:
args = []
for d in dates:
date1, date2 = d
temp = [(var, outpath, date1, date2, expe_polygons_boxbound_wgs84, grid_size) for var in variables]
args.extend(temp)
# Download in multiprocess
weather_files = p_map(era5land.call_era5landhourly, args, **{"num_cpus": nb_processes})
Jeremy Auclair
committed
print('\nPath to daily weather netCDF file: ', outpath)
Jeremy Auclair
committed
Jeremy Auclair
committed
# Generate pandas dataframe for parcel mode
if mode == 'parcel':
# Create save path
weather_datframe = weather_daily_ncFile + '_df.csv'
weather_dataset = weather_daily_ncFile + '_parcel.nc'
Jeremy Auclair
committed
# Check if weather dataset already exists

Jérémy AUCLAIR
committed
if os.path.exists(weather_dataset) and not weather_overwrite:
Jeremy Auclair
committed
return weather_dataset
# Generate daily weather datasets as Geotiffs for each variable
weather_daily_rain, weather_daily_ET0 = era5land.era5Land_daily_to_yearly_parcel(weather_files, variables, weather_daily_ncFile, start_date, end_date, h = wind_height)
Jeremy Auclair
committed
Jeremy Auclair
committed
# Generate and save weather dataframe

Jérémy AUCLAIR
committed
era5land.extract_weather_dataframe(weather_daily_rain, weather_daily_ET0, shapefile, config_file, weather_datframe, max_cpu = config_params.max_cpu)
Jeremy Auclair
committed
# Convert dataframe to xarray dataset
convert_dataframe_to_xarray(weather_datframe, weather_dataset, variables = ['Rain', 'ET0'], data_types = ['u2', 'u2'])
Jeremy Auclair
committed
print('\nWeather dataset:', weather_daily_ncFile)
Jeremy Auclair
committed
return weather_dataset
Jeremy Auclair
committed
# Calculate ET0 over the whole time period
weather_daily_ncFile = era5land.era5Land_daily_to_yearly_pixel(weather_files, variables, weather_daily_ncFile, raw_S2_image_ref, ndvi_path, start_date, end_date, h = wind_height, max_ram = config_params.max_ram, weather_overwrite = weather_overwrite)
Jeremy Auclair
committed
print('\nWeather dataset:', weather_daily_ncFile)
Jeremy Auclair
committed
Jeremy Auclair
committed
return weather_daily_ncFile