From ff9751395dec3c8412a967e22363395de907d664 Mon Sep 17 00:00:00 2001 From: "Timothy W. Hilton" Date: Thu, 5 Sep 2024 16:08:13 +1200 Subject: [PATCH 1/4] - branch to change Munich example to NZ split NZ ERA5 annual netCDF files to monthly --- vprm_predictions/era5_to_monthly.py | 30 +++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 vprm_predictions/era5_to_monthly.py diff --git a/vprm_predictions/era5_to_monthly.py b/vprm_predictions/era5_to_monthly.py new file mode 100644 index 0000000..8a79f78 --- /dev/null +++ b/vprm_predictions/era5_to_monthly.py @@ -0,0 +1,30 @@ +"""separate annual ERA5 netCDF files into monthly""" + +import socket +from pathlib import Path + +import pandas as pd +import xarray as xr + +from loguru import logger + +ERA5_DATA_DIR = Path("/mnt/amp/CarbonWatchUrban/urbanVPRM/ERA5/") +if socket.gethostname() == "hutl21264.gns.cri.nz": + # prepend /home/timh to path for laptop + ERA5_DATA_DIR = Path("/home/timh", *ERA5_DATA_DIR.parts[1:]) + +if __name__ == "__main__": + this_year = 2023 + for this_year in range(2020, 2024): + for this_month in range(1, 13): + fname = f"ERA5_2mT_msdswrf_ssrd_NZ_{this_year}.nc" + xds = xr.open_dataset(Path(ERA5_DATA_DIR, fname)) + month = pd.to_datetime(xds["time"]).month + + fname = Path( + ERA5_DATA_DIR + / "monthly" + / f"ERA5_2mT_msdswrf_ssrd_NZ_{this_year}_{this_month:02d}.nc" + ) + logger.info(f"writing {fname}") + # xds.sel(time=(month == this_month)).to_netcdf(fname, mode="w") From 153f95e0b6cd460f64e1d661d2b8e18e4f5d3c1f Mon Sep 17 00:00:00 2001 From: "Timothy W. Hilton" Date: Fri, 6 Sep 2024 11:39:28 +1200 Subject: [PATCH 2/4] conda environment spec --- conda_environment.yaml | 47 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 conda_environment.yaml diff --git a/conda_environment.yaml b/conda_environment.yaml new file mode 100644 index 0000000..54d7f95 --- /dev/null +++ b/conda_environment.yaml @@ -0,0 +1,47 @@ +name: pyVPRM +dependencies: + - dask + - netCDF4 + - esmf + - esmpy + + - tqdm + - gitsummary + - loguru + + # plotting + - hvplot + - geoviews + - nzgeom + - matplotlib + + # development + - conda-build # for `mamba develop` + - isort + - jupyterlab + - jupyter_console + - jupyterlab_code_formatter + - tox + - pytest + - black + - pylint + - mypy + - pyright + - ipdb + - sphinx + - sphinx_rtd_theme + - sphinx-autodoc-typehints + - m2r2 + + # # github dependencies + # - pip + # - pip: + # - "--editable=git+https://github.com/Timothy-W-Hilton/pyVPRM" + + # - pip + # - pip: + # - "git+https://github.com/Timothy-W-Hilton/pyVPRM" + + - pip + - pip: + - "--editable=/home/timh/Code/pyVPRM" From d2f7945f11957ab743348ab176127f6346dfd089 Mon Sep 17 00:00:00 2001 From: "Timothy W. Hilton" Date: Mon, 23 Sep 2024 13:59:14 +1200 Subject: [PATCH 3/4] blacken python source files apply black formatter (https://black.readthedocs.io/en/stable/index.html) to Python source files. Ran black with default options in top-level pyVPRM directory (`black --verbose .`). versions used: black, 24.8.0 (compiled: no) Python (CPython) 3.12.5 --- fit_vprm_parameters/fit_params_draft.py | 68 ++++-- fit_vprm_parameters/get_fit_files_for_site.py | 141 +++++++---- .../download_satellite_images.py | 58 ++--- vprm_predictions/make_vprm_predictions.py | 168 +++++++------ wrf_preprocessor/vprm_preprocessor_new.py | 226 +++++++++++------- 5 files changed, 399 insertions(+), 262 deletions(-) diff --git a/fit_vprm_parameters/fit_params_draft.py b/fit_vprm_parameters/fit_params_draft.py index 0353931..730142c 100644 --- a/fit_vprm_parameters/fit_params_draft.py +++ b/fit_vprm_parameters/fit_params_draft.py @@ -1,9 +1,10 @@ import sys import os import pathlib -sys.path.append(os.path.join(pathlib.Path(__file__).parent.resolve(), '..', 'lib')) + +sys.path.append(os.path.join(pathlib.Path(__file__).parent.resolve(), "..", "lib")) from sat_manager import VIIRS, modis, copernicus_land_cover_map -from VPRM import vprm +from VPRM import vprm import yaml import glob import time @@ -13,8 +14,8 @@ import argparse p = argparse.ArgumentParser( - description = "Commend Line Arguments", - formatter_class = argparse.RawTextHelpFormatter) + description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter +) p.add_argument("--h", type=int) p.add_argument("--v", type=int) p.add_argument("--config", type=str) @@ -25,17 +26,17 @@ with open(args.config, "r") as stream: try: - cfg = yaml.safe_load(stream) + cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: print(exc) h = args.h v = args.v -if not os.path.exists(cfg['out_path']): - os.makedirs(cfg['out_path']) +if not os.path.exists(cfg["out_path"]): + os.makedirs(cfg["out_path"]) -outfile = os.path.join(cfg['out_path'], 'h{}v{}_{}.pickle'.format(h, v, cfg['year'])) +outfile = os.path.join(cfg["out_path"], "h{}v{}_{}.pickle".format(h, v, cfg["year"])) print(outfile) @@ -51,25 +52,37 @@ vprm_inst = vprm(n_cpus=args.n_cpus) -for c, i in enumerate(glob.glob(os.path.join(cfg['sat_image_path'], '*h{:02d}v{:02d}*.h*'.format(h, v)))): +for c, i in enumerate( + glob.glob(os.path.join(cfg["sat_image_path"], "*h{:02d}v{:02d}*.h*".format(h, v))) +): print(i) - if cfg['satellite'] == 'modis': + if cfg["satellite"] == "modis": handler = modis(sat_image_path=i) handler.load() - vprm_inst.add_sat_img(handler, b_nir='B02', b_red='B01', - b_blue='B03', b_swir='B06', - which_evi='evi', - drop_bands=True) + vprm_inst.add_sat_img( + handler, + b_nir="B02", + b_red="B01", + b_blue="B03", + b_swir="B06", + which_evi="evi", + drop_bands=True, + ) else: handler = VIIRS(sat_image_path=i) handler.load() - vprm_inst.add_sat_img(handler, b_nir='SurfReflect_I2', b_red='SurfReflect_I1', - b_blue='no_blue_sensor', b_swir='SurfReflect_I3', - which_evi='evi2', - drop_bands=True) + vprm_inst.add_sat_img( + handler, + b_nir="SurfReflect_I2", + b_red="SurfReflect_I1", + b_blue="no_blue_sensor", + b_swir="SurfReflect_I3", + which_evi="evi2", + drop_bands=True, + ) # Uncomment if you want to do a spatial smearing with size=(3,3). Otherwise give size manually. -#vprm_inst.smearing(lonlats=lonlats) +# vprm_inst.smearing(lonlats=lonlats) # Sort the satellite images by time and merge internally for easier computations vprm_inst.sort_and_merge_by_timestamp() @@ -84,15 +97,18 @@ # Loop over the year and get the vprm_variables -time_range = pd.date_range(start="{}-01-01".format(cfg['year']), - end="{}-01-01".format(cfg['year'] + 1), - freq='H') +time_range = pd.date_range( + start="{}-01-01".format(cfg["year"]), + end="{}-01-01".format(cfg["year"] + 1), + freq="H", +) -for c,t in enumerate(time_range): +for c, t in enumerate(time_range): t0 = time.time() # Returns dict with VPRM input variables for the different lats and lons # Additional ERA5 variables can be added trough the add_era_variables key - vrbls = vprm_inst.get_vprm_variables(t, lat=lats, lon=lons, - add_era_variables=['svwl1']) - + vrbls = vprm_inst.get_vprm_variables( + t, lat=lats, lon=lons, add_era_variables=["svwl1"] + ) + # ToDo: write to your own format and store for further analysis diff --git a/fit_vprm_parameters/get_fit_files_for_site.py b/fit_vprm_parameters/get_fit_files_for_site.py index bcbb16f..e44d314 100644 --- a/fit_vprm_parameters/get_fit_files_for_site.py +++ b/fit_vprm_parameters/get_fit_files_for_site.py @@ -2,7 +2,7 @@ from pyVPRM.lib.flux_tower_class import fluxnet, icos from pyVPRM.sat_managers.viirs import VIIRS from pyVPRM.sat_managers.modis import modis -from pyVPRM.VPRM import vprm +from pyVPRM.VPRM import vprm import xarray as xr import pickle import yaml @@ -16,6 +16,7 @@ from astropy.convolution import Gaussian2DKernel from functions import lat_lon_to_modis + def all_files_exist(item): for key in item.assets.keys(): path = str(item.assets[key].href[7:]) @@ -23,9 +24,10 @@ def all_files_exist(item): return False return True + p = argparse.ArgumentParser( - description = "Commend Line Arguments", - formatter_class = argparse.RawTextHelpFormatter) + description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter +) p.add_argument("--site", type=str) p.add_argument("--veg_type", type=str) p.add_argument("--this_year", type=str) @@ -34,100 +36,137 @@ def all_files_exist(item): with open(args.cfg_path, "r") as stream: try: - cfg = yaml.safe_load(stream) + cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: print(exc) s = args.site veg_type = args.veg_type -this_year = int(args.this_year) +this_year = int(args.this_year) -veg_type_id = {'GRA': 7, 'MF': 3 , 'CRO':6, - 'EF': 1, 'DF': 2, 'SH':4, 'WET': 9} +veg_type_id = {"GRA": 7, "MF": 3, "CRO": 6, "EF": 1, "DF": 2, "SH": 4, "WET": 9} # GRA = Grassland | MF = Mixed Forest | CRO = Cropland | EF = Evergreen Forest | DF = Deciduous Forest | SH = Shrubland all_data = [] vprm_insts = [] -tmin = parser.parse('{}0101'.format(this_year)) -tmax = parser.parse('{}1231'.format(this_year)) +tmin = parser.parse("{}0101".format(this_year)) +tmax = parser.parse("{}1231".format(this_year)) tower_data_list = [] print(s) # if len(tower_data_list) >1: # continue if this_year == 2012: - data_files = glob.glob(os.path.join(cfg['fluxnet_path'], '*{}*/*_FULLSET_H*'.format(s))) - if len(data_files) == 0: exit() - flux_tower_inst = fluxnet(data_files[0],'SW_IN_F', 'TA_F', - t_start= tmin, t_stop=tmax) + data_files = glob.glob( + os.path.join(cfg["fluxnet_path"], "*{}*/*_FULLSET_H*".format(s)) + ) + if len(data_files) == 0: + exit() + flux_tower_inst = fluxnet( + data_files[0], "SW_IN_F", "TA_F", t_start=tmin, t_stop=tmax + ) elif this_year in [2021, 2022]: - data_files = glob.glob(os.path.join(cfg['icos_path'], '*/*_{}_FLUXNET_HH_L2.csv'.format(s))) - if len(data_files) == 0: exit() - flux_tower_inst = icos(data_files[0],'SW_IN_F', 'TA_F', t_start= tmin, t_stop=tmax) + data_files = glob.glob( + os.path.join(cfg["icos_path"], "*/*_{}_FLUXNET_HH_L2.csv".format(s)) + ) + if len(data_files) == 0: + exit() + flux_tower_inst = icos(data_files[0], "SW_IN_F", "TA_F", t_start=tmin, t_stop=tmax) lon = flux_tower_inst.get_lonlat()[0] -lat= flux_tower_inst.get_lonlat()[1] +lat = flux_tower_inst.get_lonlat()[1] okay = flux_tower_inst.add_tower_data() -if not okay: exit() +if not okay: + exit() print(flux_tower_inst.get_site_name()) flux_tower_inst.set_land_type(veg_type_id[veg_type]) lat = flux_tower_inst.get_lonlat()[1] lon = flux_tower_inst.get_lonlat()[0] h, v = lat_lon_to_modis(lat, lon) -inp_files1 = sorted(glob.glob(os.path.join(cfg['sat_image_path'], str(this_year-1), - '*h{:02d}v{:02d}*.h*'.format(h, v))))[-8:] -inp_files2 = sorted(glob.glob(os.path.join(cfg['sat_image_path'], str(this_year), - '*h{:02d}v{:02d}*.h*'.format(h, v)))) -inp_files3 = sorted(glob.glob(os.path.join(cfg['sat_image_path'], str(this_year+1), - '*h{:02d}v{:02d}*.h*'.format(h, v))))[:8] +inp_files1 = sorted( + glob.glob( + os.path.join( + cfg["sat_image_path"], + str(this_year - 1), + "*h{:02d}v{:02d}*.h*".format(h, v), + ) + ) +)[-8:] +inp_files2 = sorted( + glob.glob( + os.path.join( + cfg["sat_image_path"], str(this_year), "*h{:02d}v{:02d}*.h*".format(h, v) + ) + ) +) +inp_files3 = sorted( + glob.glob( + os.path.join( + cfg["sat_image_path"], + str(this_year + 1), + "*h{:02d}v{:02d}*.h*".format(h, v), + ) + ) +)[:8] inp_files = np.concatenate([inp_files1, inp_files2, inp_files3]) -inp_files = np.array([i for i in inp_files if '.xml' not in i]) -_, inds = np.unique([os.path.basename(f) for f in inp_files], - return_index=True) +inp_files = np.array([i for i in inp_files if ".xml" not in i]) +_, inds = np.unique([os.path.basename(f) for f in inp_files], return_index=True) inp_files = inp_files[inds] vprm_inst = vprm(n_cpus=1, sites=[flux_tower_inst]) for c, i in enumerate(inp_files): print(i) - if cfg['satellite'] == 'modis': + if cfg["satellite"] == "modis": handler = modis(sat_image_path=i) handler.load() handler.crop(lonlat=(lon, lat), radius=4) - vprm_inst.add_sat_img(handler, b_nir='sur_refl_b02', - b_red='sur_refl_b01', - b_blue='sur_refl_b03',b_swir='sur_refl_b06', - which_evi='evi', - timestamp_key='sur_refl_day_of_year', - mask_bad_pixels=True, - mask_clouds=True, - mask_snow=True, - drop_bands=True) + vprm_inst.add_sat_img( + handler, + b_nir="sur_refl_b02", + b_red="sur_refl_b01", + b_blue="sur_refl_b03", + b_swir="sur_refl_b06", + which_evi="evi", + timestamp_key="sur_refl_day_of_year", + mask_bad_pixels=True, + mask_clouds=True, + mask_snow=True, + drop_bands=True, + ) else: handler = VIIRS(sat_image_path=i) handler.load() handler.crop(lonlat=(lon, lat), radius=4) - vprm_inst.add_sat_img(handler, b_nir='SurfReflect_I2', b_red='SurfReflect_I1', - b_blue='no_blue_sensor', b_swir='SurfReflect_I3', - which_evi='evi2', - drop_bands=True) + vprm_inst.add_sat_img( + handler, + b_nir="SurfReflect_I2", + b_red="SurfReflect_I1", + b_blue="no_blue_sensor", + b_swir="SurfReflect_I3", + which_evi="evi2", + drop_bands=True, + ) # Uncomment if you want to do a spatial smearing with size=(3,3). Otherwise give size manually. # Sort the satellite images by time and merge internally for easier computations vprm_inst.sort_and_merge_by_timestamp() -vprm_inst.smearing(lonlats=[(lon, lat)], - keys=['evi', 'lswi'], - kernel=Gaussian2DKernel(1)) #(21, 21)) +vprm_inst.smearing( + lonlats=[(lon, lat)], keys=["evi", "lswi"], kernel=Gaussian2DKernel(1) +) # (21, 21)) vprm_inst.reduce_along_lat_lon() # Run lowess smoothing -vprm_inst.lowess(times=[datetime(this_year, 1, 1)+ timedelta(days=i) for i in np.arange(365.)], - frac=0.15, it=3, - keys=['evi', 'lswi']) +vprm_inst.lowess( + times=[datetime(this_year, 1, 1) + timedelta(days=i) for i in np.arange(365.0)], + frac=0.15, + it=3, + keys=["evi", "lswi"], +) # Calculate necessary parameters for the vprm calculation vprm_inst.calc_min_max_evi_lswi() @@ -137,10 +176,12 @@ def all_files_exist(item): # Loop over the year and get the vprm_variables data_list = vprm_inst.data_for_fitting() -opath = '/home/b/b309233/software/VPRM_preprocessor/analysis_scripts/site_data_for_fit_modis/{}_{}.pickle'.format(this_year, s) +opath = "/home/b/b309233/software/VPRM_preprocessor/analysis_scripts/site_data_for_fit_modis/{}_{}.pickle".format( + this_year, s +) if os.path.exists(opath): os.remove(opath) - -with open(opath, 'wb+') as ofile: + +with open(opath, "wb+") as ofile: pickle.dump(data_list, ofile) diff --git a/sat_data_download/download_satellite_images.py b/sat_data_download/download_satellite_images.py index b70dfd2..9caf37e 100644 --- a/sat_data_download/download_satellite_images.py +++ b/sat_data_download/download_satellite_images.py @@ -1,14 +1,14 @@ import os from pyVPRM.sat_managers.viirs import VIIRS from pyVPRM.sat_managers.modis import modis -import yaml +import yaml from datetime import date import argparse import shutil p = argparse.ArgumentParser( - description = "Commend Line Arguments", - formatter_class = argparse.RawTextHelpFormatter) + description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter +) p.add_argument("--config", type=str) p.add_argument("--login_data", type=str) p.add_argument("--year", type=int, default=None) @@ -18,56 +18,60 @@ with open(args.config, "r") as stream: try: - cfg = yaml.safe_load(stream) + cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: print(exc) with open(args.login_data, "r") as stream: try: - logins = yaml.safe_load(stream) + logins = yaml.safe_load(stream) except yaml.YAMLError as exc: print(exc) if args.year is not None: years = [args.year] else: - years = cfg['years'] + years = cfg["years"] -hvs=cfg['hvs'] +hvs = cfg["hvs"] if (args.h is not None) & (args.v is not None): hvs = [(args.h, args.v)] for year in years: - savepath = os.path.join(cfg['sat_image_path'], str(year)) - if cfg['satellite'] == 'modis': + savepath = os.path.join(cfg["sat_image_path"], str(year)) + if cfg["satellite"] == "modis": for i in hvs: - print('Tile {}'.format(i)) + print("Tile {}".format(i)) handler = modis() try: - handler.download(date(year, 1, 1), - savepath = savepath, - username = logins['modis'][0], - pwd = logins['modis'][1], - hv=i, - delta = 1, - enddate=date(year + 1 , 1, 1)) + handler.download( + date(year, 1, 1), + savepath=savepath, + username=logins["modis"][0], + pwd=logins["modis"][1], + hv=i, + delta=1, + enddate=date(year + 1, 1, 1), + ) except Exception as e: print(e) - elif cfg['satellite'] == 'viirs': + elif cfg["satellite"] == "viirs": for i in hvs: - print('Tile {}'.format(i)) + print("Tile {}".format(i)) handler = VIIRS() try: - handler.download(date(year, 1, 1), - savepath = savepath, - username = logins['modis'][0], - pwd = logins['modis'][1], - hv = i, - delta = 1, - enddate=date(year + 1, 1, 1)) + handler.download( + date(year, 1, 1), + savepath=savepath, + username=logins["modis"][0], + pwd=logins["modis"][1], + hv=i, + delta=1, + enddate=date(year + 1, 1, 1), + ) except Exception as e: print(e) else: - print('No download function for this satellite implemented') + print("No download function for this satellite implemented") diff --git a/vprm_predictions/make_vprm_predictions.py b/vprm_predictions/make_vprm_predictions.py index 5dfe645..e2f29aa 100644 --- a/vprm_predictions/make_vprm_predictions.py +++ b/vprm_predictions/make_vprm_predictions.py @@ -4,7 +4,7 @@ from pyVPRM.sat_managers.viirs import VIIRS from pyVPRM.sat_managers.modis import modis from pyVPRM.sat_managers.copernicus import copernicus_land_cover_map -from pyVPRM.VPRM import vprm +from pyVPRM.VPRM import vprm from pyVPRM.meteorologies import era5_monthly_xr, era5_class_dkrz from pyVPRM.lib.functions import lat_lon_to_modis from pyVPRM.vprm_models import vprm_modified, vprm_base @@ -24,9 +24,13 @@ def get_hourly_time_range(year, day_of_year): - - start_time = datetime(year, 1, 1) + timedelta(days=int(day_of_year)-1) # Set the starting time based on the day of the year - end_time = start_time + timedelta(hours=1) # Add 1 hour to get the end time of the first hour + + start_time = datetime(year, 1, 1) + timedelta( + days=int(day_of_year) - 1 + ) # Set the starting time based on the day of the year + end_time = start_time + timedelta( + hours=1 + ) # Add 1 hour to get the end time of the first hour hourly_range = [] while start_time.timetuple().tm_yday == day_of_year: @@ -35,72 +39,88 @@ def get_hourly_time_range(year, day_of_year): end_time = start_time + timedelta(hours=1) return hourly_range + # Read command line arguments p = argparse.ArgumentParser( - description = "Commend Line Arguments", - formatter_class = argparse.RawTextHelpFormatter) + description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter +) p.add_argument("--h", type=int) p.add_argument("--v", type=int) p.add_argument("--config", type=str) p.add_argument("--n_cpus", type=int, default=1) p.add_argument("--year", type=int) args = p.parse_args() -print('Run with args', args) +print("Run with args", args) h = args.h v = args.v -#Read config +# Read config with open(args.config, "r") as stream: try: - cfg = yaml.safe_load(stream) + cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: print(exc) -if not os.path.exists(cfg['predictions_path']): - os.makedirs(cfg['predictions_path']) +if not os.path.exists(cfg["predictions_path"]): + os.makedirs(cfg["predictions_path"]) # Initialize VPRM instance with the copernicus land cover config -vprm_inst = vprm(vprm_config_path=os.path.join(pyVPRM.__path__[0], 'vprm_configs/copernicus_land_cover.yaml'), - n_cpus=args.n_cpus) - -#Note: There is no need to convert HDF4 into Netcdf files. You can also use HDF4 files directly. -files = glob.glob(os.path.join(cfg['sat_image_path'],# str(args.year), - '*h{:02d}v{:02d}*.nc'.format(h, v))) +vprm_inst = vprm( + vprm_config_path=os.path.join( + pyVPRM.__path__[0], "vprm_configs/copernicus_land_cover.yaml" + ), + n_cpus=args.n_cpus, +) + +# Note: There is no need to convert HDF4 into Netcdf files. You can also use HDF4 files directly. +files = glob.glob( + os.path.join( + cfg["sat_image_path"], "*h{:02d}v{:02d}*.nc".format(h, v) # str(args.year), + ) +) # Add satellite images to the VPRM instance for c, i in enumerate(sorted(files)): - if '.xml' in i: + if ".xml" in i: continue print(i) - if cfg['satellite'] == 'modis': + if cfg["satellite"] == "modis": handler = modis(sat_image_path=i) handler.load() - vprm_inst.add_sat_img(handler, b_nir='sur_refl_b02', b_red='sur_refl_b01', - b_blue='sur_refl_b03', b_swir='sur_refl_b06', - which_evi='evi', - drop_bands=True, - timestamp_key='sur_refl_day_of_year', - mask_bad_pixels=True, - mask_clouds=True) + vprm_inst.add_sat_img( + handler, + b_nir="sur_refl_b02", + b_red="sur_refl_b01", + b_blue="sur_refl_b03", + b_swir="sur_refl_b06", + which_evi="evi", + drop_bands=True, + timestamp_key="sur_refl_day_of_year", + mask_bad_pixels=True, + mask_clouds=True, + ) else: handler = VIIRS(sat_image_path=i) handler.load() - vprm_inst.add_sat_img(handler, b_nir='SurfReflect_I2', b_red='SurfReflect_I1', - b_blue='no_blue_sensor', b_swir='SurfReflect_I3', - which_evi='evi2', - drop_bands=True) + vprm_inst.add_sat_img( + handler, + b_nir="SurfReflect_I2", + b_red="SurfReflect_I1", + b_blue="no_blue_sensor", + b_swir="SurfReflect_I3", + which_evi="evi2", + drop_bands=True, + ) # Sort the satellite data by time and run the lowess smoothing vprm_inst.sort_and_merge_by_timestamp() -vprm_inst.lowess(keys=['evi', 'lswi'], - times='daily', - frac=0.2, it=3) +vprm_inst.lowess(keys=["evi", "lswi"], times="daily", frac=0.2, it=3) # Clip EVI and LSWI values to allows ranges -vprm_inst.clip_values('evi', 0, 1) -vprm_inst.clip_values('lswi',-1, 1) +vprm_inst.clip_values("evi", 0, 1) +vprm_inst.clip_values("lswi", -1, 1) # Calculate the minimum and maximum EVI/LSWI vprm_inst.calc_min_max_evi_lswi() @@ -108,7 +128,7 @@ def get_hourly_time_range(year, day_of_year): # Add land covery map(s) by iterating over all maps in the `copernicus path` and picking those that overlap with our satellite images lcm = None -for c in glob.glob(os.path.join(cfg['copernicus_path'], '*')): +for c in glob.glob(os.path.join(cfg["copernicus_path"], "*")): # Generate a copernicus_land_cover_map instance thandler = copernicus_land_cover_map(c) thandler.load() @@ -117,75 +137,85 @@ def get_hourly_time_range(year, day_of_year): # Check overlap with our satellite images dj = rasterio.coords.disjoint_bounds(bounds, thandler.sat_img.rio.bounds()) if dj: - print('Do not add {}'.format(c)) + print("Do not add {}".format(c)) continue - print('Add {}'.format(c)) + print("Add {}".format(c)) if lcm is None: - lcm=copernicus_land_cover_map(c) + lcm = copernicus_land_cover_map(c) lcm.load() else: lcm.add_tile(thandler, reproject=False) # Crop land cover map to the extend of our satellite images (to speed up computations and save memory) geom = box(*vprm_inst.sat_imgs.sat_img.rio.bounds()) -df = gpd.GeoDataFrame({"id":1,"geometry":[geom]}) +df = gpd.GeoDataFrame({"id": 1, "geometry": [geom]}) df = df.set_crs(vprm_inst.sat_imgs.sat_img.rio.crs) df = df.scale(1.3, 1.3) lcm.crop_to_polygon(df) # Add land cover map to the VPRM instance. This wil regrid the land cover map to the satellite grid -vprm_inst.add_land_cover_map(lcm, regridder_save_path=os.path.join(cfg['predictions_path'], - 'regridder.nc'), mpi=False) +vprm_inst.add_land_cover_map( + lcm, + regridder_save_path=os.path.join(cfg["predictions_path"], "regridder.nc"), + mpi=False, +) # Set meteorology -era5_inst = era5_monthly_xr.met_data_handler(args.year, 1, 1, 0, - './data/era5', - keys=['t2m', 'ssrd']) +era5_inst = era5_monthly_xr.met_data_handler( + args.year, 1, 1, 0, "./data/era5", keys=["t2m", "ssrd"] +) -# Load VPRM parameters from a dictionary -with open(cfg['vprm_params_dict'], 'rb') as ifile: +# Load VPRM parameters from a dictionary +with open(cfg["vprm_params_dict"], "rb") as ifile: res_dict = pickle.load(ifile) -vprm_model = vprm_base.vprm_base(vprm_pre=vprm_inst, - met=era5_inst, - fit_params_dict= res_dict) - +vprm_model = vprm_base.vprm_base( + vprm_pre=vprm_inst, met=era5_inst, fit_params_dict=res_dict +) + # Make NEE/GPP flux predictions and save them days_in_year = 365 + calendar.isleap(args.year) -met_regridder_weights = os.path.join(cfg['predictions_path'], - 'met_regridder_weights.nc') +met_regridder_weights = os.path.join( + cfg["predictions_path"], "met_regridder_weights.nc" +) -for i in np.arange(160,161, 1): - time_range=get_hourly_time_range(int(args.year), i) +for i in np.arange(160, 161, 1): + time_range = get_hourly_time_range(int(args.year), i) preds_gpp = [] preds_nee = [] ts = [] for t in time_range[:]: - t0=time.time() + t0 = time.time() print(t) - pred = vprm_model.make_vprm_predictions(t, met_regridder_weights=met_regridder_weights) + pred = vprm_model.make_vprm_predictions( + t, met_regridder_weights=met_regridder_weights + ) if pred is None: continue - preds_gpp.append(pred['gpp']) - preds_nee.append(pred['nee']) + preds_gpp.append(pred["gpp"]) + preds_nee.append(pred["nee"]) ts.append(t) - preds_gpp = xr.concat(preds_gpp, 'time') - preds_gpp = preds_gpp.assign_coords({'time': ts}) - outpath = os.path.join(cfg['predictions_path'], - 'gpp_h{:02d}v{:02d}_{}_{:03d}.h5'.format(h, v, args.year, i)) + preds_gpp = xr.concat(preds_gpp, "time") + preds_gpp = preds_gpp.assign_coords({"time": ts}) + outpath = os.path.join( + cfg["predictions_path"], + "gpp_h{:02d}v{:02d}_{}_{:03d}.h5".format(h, v, args.year, i), + ) if os.path.exists(outpath): os.remove(outpath) preds_gpp.to_netcdf(outpath) preds_gpp.close() - - preds_nee = xr.concat(preds_nee, 'time') - preds_nee = preds_nee.assign_coords({'time': ts}) - outpath = os.path.join(cfg['predictions_path'], - 'nee_h{:02d}v{:02d}_{}_{:03d}.h5'.format(h, v,args.year, i)) + + preds_nee = xr.concat(preds_nee, "time") + preds_nee = preds_nee.assign_coords({"time": ts}) + outpath = os.path.join( + cfg["predictions_path"], + "nee_h{:02d}v{:02d}_{}_{:03d}.h5".format(h, v, args.year, i), + ) if os.path.exists(outpath): os.remove(outpath) preds_nee.to_netcdf(outpath) preds_nee.close() -print('Done. In order to inspect the output use evaluate_output.ipynb') +print("Done. In order to inspect the output use evaluate_output.ipynb") diff --git a/wrf_preprocessor/vprm_preprocessor_new.py b/wrf_preprocessor/vprm_preprocessor_new.py index 21131cb..0d26be2 100644 --- a/wrf_preprocessor/vprm_preprocessor_new.py +++ b/wrf_preprocessor/vprm_preprocessor_new.py @@ -3,7 +3,11 @@ from pyVPRM.sat_managers.viirs import VIIRS from pyVPRM.sat_managers.modis import modis from pyVPRM.sat_managers.copernicus import copernicus_land_cover_map -from pyVPRM.lib.functions import lat_lon_to_modis, add_corners_to_1d_grid, parse_wrf_grid_file +from pyVPRM.lib.functions import ( + lat_lon_to_modis, + add_corners_to_1d_grid, + parse_wrf_grid_file, +) from pyVPRM.VPRM import vprm import yaml import glob @@ -18,8 +22,8 @@ # Read command line arguments p = argparse.ArgumentParser( - description = "Commend Line Arguments", - formatter_class = argparse.RawTextHelpFormatter) + description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter +) p.add_argument("--config", type=str) p.add_argument("--year", type=int) p.add_argument("--n_cpus", type=int, default=1) @@ -29,103 +33,129 @@ args = p.parse_args() print(args) -#Load config file +# Load config file this_year = int(args.year) with open(args.config, "r") as stream: try: - cfg = yaml.safe_load(stream) + cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: print(exc) -if not os.path.exists(cfg['out_path']): - os.makedirs(cfg['out_path']) +if not os.path.exists(cfg["out_path"]): + os.makedirs(cfg["out_path"]) -# Use WRF netcdf to generate an xarray instance. With chunk_x and chunk_y +# Use WRF netcdf to generate an xarray instance. With chunk_x and chunk_y # the output grid can be calculated in disjoint pieces -out_grid = parse_wrf_grid_file(cfg['geo_em_file'], n_chunks=cfg['n_chunks'], - chunk_x=args.chunk_x, chunk_y=args.chunk_y) +out_grid = parse_wrf_grid_file( + cfg["geo_em_file"], + n_chunks=cfg["n_chunks"], + chunk_x=args.chunk_x, + chunk_y=args.chunk_y, +) # Calculate the required MODIS tiles for our out grid -hvs = np.unique([lat_lon_to_modis(out_grid['lat_b'].values.flatten()[i], - out_grid['lon_b'].values.flatten()[i]) - for i in range(len(out_grid['lat_b'].values.flatten()))], - axis=0) +hvs = np.unique( + [ + lat_lon_to_modis( + out_grid["lat_b"].values.flatten()[i], out_grid["lon_b"].values.flatten()[i] + ) + for i in range(len(out_grid["lat_b"].values.flatten())) + ], + axis=0, +) print(hvs) insts = [] -#Load the data +# Load the data # days of interest -days = [datetime(this_year, 1, 1)+ timedelta(days=i) for i in np.arange(365.)] +days = [datetime(this_year, 1, 1) + timedelta(days=i) for i in np.arange(365.0)] # read the data for c, i in enumerate(hvs): - + print(i) - - #Note: There is no need to convert MODIS HDF4 into Netcdf files. You can also use HDF4 files directly. - file_collections = glob.glob(os.path.join(cfg['sat_image_path'], - '*h{:02d}v{:02d}*.nc'.format(i[0], i[1]))) - + + # Note: There is no need to convert MODIS HDF4 into Netcdf files. You can also use HDF4 files directly. + file_collections = glob.glob( + os.path.join(cfg["sat_image_path"], "*h{:02d}v{:02d}*.nc".format(i[0], i[1])) + ) + if len(file_collections) == 0: continue - new_inst = vprm(vprm_config_path = os.path.join(pyVPRM.__path__[0], - 'vprm_configs/copernicus_land_cover.yaml') , - n_cpus=args.n_cpus) + new_inst = vprm( + vprm_config_path=os.path.join( + pyVPRM.__path__[0], "vprm_configs/copernicus_land_cover.yaml" + ), + n_cpus=args.n_cpus, + ) for c0, fpath in enumerate(file_collections): print(fpath) - if cfg['satellite'] == 'modis': + if cfg["satellite"] == "modis": handler = modis(sat_image_path=fpath) - handler.load() - elif cfg['satellite'] == 'viirs': + handler.load() + elif cfg["satellite"] == "viirs": handler = VIIRS(sat_image_path=fpath) handler.load() else: - print('Set the satellite in the cfg either to modis or viirs.') + print("Set the satellite in the cfg either to modis or viirs.") - # In order to save memory crop the satellite images to our WRF out grid + # In order to save memory crop the satellite images to our WRF out grid if c0 == 0: - trans = Transformer.from_crs('+proj=longlat +datum=WGS84', - handler.sat_img.rio.crs) + trans = Transformer.from_crs( + "+proj=longlat +datum=WGS84", handler.sat_img.rio.crs + ) # To save memory crop satellite images to WRF grid - x_a, y_a = trans.transform(out_grid['lon_b'], out_grid['lat_b']) - b = box(float(np.min(x_a)), float(np.min(y_a)), - float(np.max(x_a)), float(np.max(y_a))) + x_a, y_a = trans.transform(out_grid["lon_b"], out_grid["lat_b"]) + b = box( + float(np.min(x_a)), + float(np.min(y_a)), + float(np.max(x_a)), + float(np.max(y_a)), + ) b = gpd.GeoSeries(Polygon(b), crs=handler.sat_img.rio.crs) b = b.scale(1.2, 1.2) handler.crop_box(b) # Add satellite image to VPRM instance - if cfg['satellite'] == 'modis': - new_inst.add_sat_img(handler, b_nir='sur_refl_b02', b_red='sur_refl_b01', - b_blue='sur_refl_b03', b_swir='sur_refl_b06', - which_evi='evi', - drop_bands=True, - timestamp_key='sur_refl_day_of_year', - mask_bad_pixels=True, - mask_clouds=True) - elif cfg['satellite'] == 'viirs': - new_inst.add_sat_img(handler, b_nir='SurfReflect_I2', b_red='SurfReflect_I1', - b_blue='no_blue_sensor', b_swir='SurfReflect_I3', - which_evi='evi2', - drop_bands=True) - + if cfg["satellite"] == "modis": + new_inst.add_sat_img( + handler, + b_nir="sur_refl_b02", + b_red="sur_refl_b01", + b_blue="sur_refl_b03", + b_swir="sur_refl_b06", + which_evi="evi", + drop_bands=True, + timestamp_key="sur_refl_day_of_year", + mask_bad_pixels=True, + mask_clouds=True, + ) + elif cfg["satellite"] == "viirs": + new_inst.add_sat_img( + handler, + b_nir="SurfReflect_I2", + b_red="SurfReflect_I1", + b_blue="no_blue_sensor", + b_swir="SurfReflect_I3", + which_evi="evi2", + drop_bands=True, + ) + # Sort and merge satellite images new_inst.sort_and_merge_by_timestamp() - + # Apply lowess smoothing in time - new_inst.lowess(keys=['evi', 'lswi'], - times=days, - frac=0.25, it=3) #0.2 + new_inst.lowess(keys=["evi", "lswi"], times=days, frac=0.25, it=3) # 0.2 # Clip EVI and LSWI to the allowed range - new_inst.clip_values('evi', 0, 1) - new_inst.clip_values('lswi',-1, 1) - new_inst.sat_imgs.sat_img = new_inst.sat_imgs.sat_img[['evi', 'lswi']] + new_inst.clip_values("evi", 0, 1) + new_inst.clip_values("lswi", -1, 1) + new_inst.sat_imgs.sat_img = new_inst.sat_imgs.sat_img[["evi", "lswi"]] insts.append(new_inst) @@ -134,76 +164,92 @@ vprm_inst = insts[0] if len(insts) > 1: vprm_inst.add_vprm_insts(insts[1:]) - + print(vprm_inst.sat_imgs.sat_img) # Add the land cover map -if not os.path.exists(cfg['out_path']): - os.makedirs(cfg['out_path']) +if not os.path.exists(cfg["out_path"]): + os.makedirs(cfg["out_path"]) # Add the land cover map and perform regridding to the satellite grid -print('Generate land cover map') -veg_file = os.path.join(cfg['out_path'], 'veg_map_on_modis_grid_{}_{}.nc'.format(args.chunk_x, - args.chunk_y)) +print("Generate land cover map") +veg_file = os.path.join( + cfg["out_path"], "veg_map_on_modis_grid_{}_{}.nc".format(args.chunk_x, args.chunk_y) +) if os.path.exists(veg_file): vprm_inst.add_land_cover_map(veg_file) else: - regridder_path = os.path.join(cfg['out_path'], 'regridder_{}_{}_lcm.nc'.format(args.chunk_x, - args.chunk_y)) + regridder_path = os.path.join( + cfg["out_path"], "regridder_{}_{}_lcm.nc".format(args.chunk_x, args.chunk_y) + ) handler_lt = None - copernicus_data_path = cfg['copernicus_path'] + copernicus_data_path = cfg["copernicus_path"] if copernicus_data_path is not None: tiles_to_add = [] - for i, c in enumerate(glob.glob(os.path.join(copernicus_data_path, '*'))): + for i, c in enumerate(glob.glob(os.path.join(copernicus_data_path, "*"))): print(c) temp_map = copernicus_land_cover_map(c) temp_map.load() dj = vprm_inst.is_disjoint(temp_map) if dj: continue - print('Add {}'.format(c)) + print("Add {}".format(c)) if handler_lt is None: handler_lt = temp_map else: tiles_to_add.append(temp_map) - + handler_lt.add_tile(tiles_to_add, reproject=False) - + # Crop the land cover map to the extend of our (cropped) satellite images to save memory geom = box(*vprm_inst.sat_imgs.sat_img.rio.bounds()) - df = gpd.GeoDataFrame({"id":1,"geometry":[geom]}) + df = gpd.GeoDataFrame({"id": 1, "geometry": [geom]}) df = df.set_crs(vprm_inst.sat_imgs.sat_img.rio.crs) df = df.scale(1.2, 1.2) - handler_lt.crop_to_polygon(df) - vprm_inst.add_land_cover_map(handler_lt, save_path=veg_file, - regridder_save_path=regridder_path) + handler_lt.crop_to_polygon(df) + vprm_inst.add_land_cover_map( + handler_lt, save_path=veg_file, regridder_save_path=regridder_path + ) -regridder_path = os.path.join(cfg['out_path'], 'regridder_{}_{}.nc'.format(args.chunk_x, - args.chunk_y)) +regridder_path = os.path.join( + cfg["out_path"], "regridder_{}_{}.nc".format(args.chunk_x, args.chunk_y) +) # Use all the information in the VPRM instance to generate the WRF input files -print('Create regridder') -wrf_op = vprm_inst.to_wrf_output(out_grid, driver = 'xEMSF', # currently only xESMF works here when the WRF grid is 2D - regridder_save_path=regridder_path, - mpi=False) +print("Create regridder") +wrf_op = vprm_inst.to_wrf_output( + out_grid, + driver="xEMSF", # currently only xESMF works here when the WRF grid is 2D + regridder_save_path=regridder_path, + mpi=False, +) # Save to NetCDF files -file_base = 'VPRM_input_' -filename_dict = {'lswi': 'LSWI', 'evi': 'EVI', 'veg_fraction': 'VEG_FRA', - 'lswi_max': 'LSWI_MAX', 'lswi_min': 'LSWI_MIN', - 'evi_max': 'EVI_MAX', 'evi_min': 'EVI_MIN'} +file_base = "VPRM_input_" +filename_dict = { + "lswi": "LSWI", + "evi": "EVI", + "veg_fraction": "VEG_FRA", + "lswi_max": "LSWI_MAX", + "lswi_min": "LSWI_MIN", + "evi_max": "EVI_MAX", + "evi_min": "EVI_MIN", +} for key in wrf_op.keys(): - ofile = os.path.join(cfg['out_path'],file_base + filename_dict[key] +'_{}_part_{}_{}.nc'.format(this_year, - args.chunk_x, - args.chunk_y)) + ofile = os.path.join( + cfg["out_path"], + file_base + + filename_dict[key] + + "_{}_part_{}_{}.nc".format(this_year, args.chunk_x, args.chunk_y), + ) if os.path.exists(ofile): os.remove(ofile) - if ('lswi' in key) | ('evi' in key): - t = wrf_op[key][key].loc[{'vprm_classes': 8}].values + if ("lswi" in key) | ("evi" in key): + t = wrf_op[key][key].loc[{"vprm_classes": 8}].values t[~np.isfinite(t)] = 0 - wrf_op[key][key].loc[{'vprm_classes': 8}] = t + wrf_op[key][key].loc[{"vprm_classes": 8}] = t wrf_op[key].to_netcdf(ofile) -print('Done. In order to inspect the output use evaluate_wrf_input.ipynb') +print("Done. In order to inspect the output use evaluate_wrf_input.ipynb") From d57eb1239f3f161d019b5325ed29565489e18871 Mon Sep 17 00:00:00 2001 From: "Timothy W. Hilton" Date: Mon, 23 Sep 2024 15:45:40 +1200 Subject: [PATCH 4/4] use loguru package for logging replace "print(...)" with "loguru.logger.info(...)" provides time-stamps and source line attributions to progress updates. The time stamps are particularly useful for understanding which pieces are taking a long time on a long-running simulation. https://loguru.readthedocs.io/en/stable/index.html https://github.com/Delgan/loguru --- fit_vprm_parameters/fit_params_draft.py | 9 ++++--- fit_vprm_parameters/get_fit_files_for_site.py | 9 ++++--- .../download_satellite_images.py | 15 +++++------ vprm_predictions/make_vprm_predictions.py | 15 +++++------ wrf_preprocessor/vprm_preprocessor_new.py | 25 ++++++++++--------- 5 files changed, 39 insertions(+), 34 deletions(-) diff --git a/fit_vprm_parameters/fit_params_draft.py b/fit_vprm_parameters/fit_params_draft.py index 730142c..e4acd2a 100644 --- a/fit_vprm_parameters/fit_params_draft.py +++ b/fit_vprm_parameters/fit_params_draft.py @@ -12,6 +12,7 @@ import pandas as pd import pickle import argparse +from loguru import logger p = argparse.ArgumentParser( description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter @@ -21,14 +22,14 @@ p.add_argument("--config", type=str) p.add_argument("--n_cpus", type=int, default=1) args = p.parse_args() -print(args) +logger.info(str(args)) with open(args.config, "r") as stream: try: cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: - print(exc) + logger.info(exc) h = args.h v = args.v @@ -37,7 +38,7 @@ os.makedirs(cfg["out_path"]) outfile = os.path.join(cfg["out_path"], "h{}v{}_{}.pickle".format(h, v, cfg["year"])) -print(outfile) +logger.info(outfile) ### ToDo: Provide a list of lats and lons for the data extraction, i.e. the flux tower positions @@ -55,7 +56,7 @@ for c, i in enumerate( glob.glob(os.path.join(cfg["sat_image_path"], "*h{:02d}v{:02d}*.h*".format(h, v))) ): - print(i) + logger.info(i) if cfg["satellite"] == "modis": handler = modis(sat_image_path=i) handler.load() diff --git a/fit_vprm_parameters/get_fit_files_for_site.py b/fit_vprm_parameters/get_fit_files_for_site.py index e44d314..8e84c9d 100644 --- a/fit_vprm_parameters/get_fit_files_for_site.py +++ b/fit_vprm_parameters/get_fit_files_for_site.py @@ -15,6 +15,7 @@ import datetime from astropy.convolution import Gaussian2DKernel from functions import lat_lon_to_modis +from loguru import logger def all_files_exist(item): @@ -38,7 +39,7 @@ def all_files_exist(item): try: cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: - print(exc) + logger.info(exc) s = args.site veg_type = args.veg_type @@ -53,7 +54,7 @@ def all_files_exist(item): tmin = parser.parse("{}0101".format(this_year)) tmax = parser.parse("{}1231".format(this_year)) tower_data_list = [] -print(s) +logger.info(s) # if len(tower_data_list) >1: # continue if this_year == 2012: @@ -77,7 +78,7 @@ def all_files_exist(item): okay = flux_tower_inst.add_tower_data() if not okay: exit() -print(flux_tower_inst.get_site_name()) +logger.info(flux_tower_inst.get_site_name()) flux_tower_inst.set_land_type(veg_type_id[veg_type]) lat = flux_tower_inst.get_lonlat()[1] @@ -117,7 +118,7 @@ def all_files_exist(item): vprm_inst = vprm(n_cpus=1, sites=[flux_tower_inst]) for c, i in enumerate(inp_files): - print(i) + logger.info(i) if cfg["satellite"] == "modis": handler = modis(sat_image_path=i) handler.load() diff --git a/sat_data_download/download_satellite_images.py b/sat_data_download/download_satellite_images.py index 9caf37e..72a82b4 100644 --- a/sat_data_download/download_satellite_images.py +++ b/sat_data_download/download_satellite_images.py @@ -5,6 +5,7 @@ from datetime import date import argparse import shutil +from loguru import logger p = argparse.ArgumentParser( description="Commend Line Arguments", formatter_class=argparse.RawTextHelpFormatter @@ -20,13 +21,13 @@ try: cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: - print(exc) + logger.info(exc) with open(args.login_data, "r") as stream: try: logins = yaml.safe_load(stream) except yaml.YAMLError as exc: - print(exc) + logger.info(exc) if args.year is not None: years = [args.year] @@ -41,7 +42,7 @@ savepath = os.path.join(cfg["sat_image_path"], str(year)) if cfg["satellite"] == "modis": for i in hvs: - print("Tile {}".format(i)) + logger.info("Tile {}".format(i)) handler = modis() try: handler.download( @@ -54,11 +55,11 @@ enddate=date(year + 1, 1, 1), ) except Exception as e: - print(e) + logger.info(e) elif cfg["satellite"] == "viirs": for i in hvs: - print("Tile {}".format(i)) + logger.info("Tile {}".format(i)) handler = VIIRS() try: handler.download( @@ -71,7 +72,7 @@ enddate=date(year + 1, 1, 1), ) except Exception as e: - print(e) + logger.info(e) else: - print("No download function for this satellite implemented") + logger.info("No download function for this satellite implemented") diff --git a/vprm_predictions/make_vprm_predictions.py b/vprm_predictions/make_vprm_predictions.py index e2f29aa..15fdf59 100644 --- a/vprm_predictions/make_vprm_predictions.py +++ b/vprm_predictions/make_vprm_predictions.py @@ -21,6 +21,7 @@ from datetime import datetime, timedelta from shapely.geometry import box import geopandas as gpd +from loguru import logger def get_hourly_time_range(year, day_of_year): @@ -50,7 +51,7 @@ def get_hourly_time_range(year, day_of_year): p.add_argument("--n_cpus", type=int, default=1) p.add_argument("--year", type=int) args = p.parse_args() -print("Run with args", args) +logger.info("Run with args: " + str(args)) h = args.h v = args.v @@ -60,7 +61,7 @@ def get_hourly_time_range(year, day_of_year): try: cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: - print(exc) + logger.info(exc) if not os.path.exists(cfg["predictions_path"]): os.makedirs(cfg["predictions_path"]) @@ -85,7 +86,7 @@ def get_hourly_time_range(year, day_of_year): for c, i in enumerate(sorted(files)): if ".xml" in i: continue - print(i) + logger.info(i) if cfg["satellite"] == "modis": handler = modis(sat_image_path=i) handler.load() @@ -137,9 +138,9 @@ def get_hourly_time_range(year, day_of_year): # Check overlap with our satellite images dj = rasterio.coords.disjoint_bounds(bounds, thandler.sat_img.rio.bounds()) if dj: - print("Do not add {}".format(c)) + logger.info("Do not add {}".format(c)) continue - print("Add {}".format(c)) + logger.info("Add {}".format(c)) if lcm is None: lcm = copernicus_land_cover_map(c) lcm.load() @@ -186,7 +187,7 @@ def get_hourly_time_range(year, day_of_year): ts = [] for t in time_range[:]: t0 = time.time() - print(t) + logger.info(t) pred = vprm_model.make_vprm_predictions( t, met_regridder_weights=met_regridder_weights ) @@ -218,4 +219,4 @@ def get_hourly_time_range(year, day_of_year): preds_nee.to_netcdf(outpath) preds_nee.close() -print("Done. In order to inspect the output use evaluate_output.ipynb") +logger.info("Done. In order to inspect the output use evaluate_output.ipynb") diff --git a/wrf_preprocessor/vprm_preprocessor_new.py b/wrf_preprocessor/vprm_preprocessor_new.py index 0d26be2..fe94b88 100644 --- a/wrf_preprocessor/vprm_preprocessor_new.py +++ b/wrf_preprocessor/vprm_preprocessor_new.py @@ -19,6 +19,7 @@ import geopandas as gpd from datetime import datetime, timedelta from pyproj import Transformer +from loguru import logger # Read command line arguments p = argparse.ArgumentParser( @@ -31,7 +32,7 @@ p.add_argument("--chunk_y", type=int, default=1) args = p.parse_args() -print(args) +logger.info(str(args)) # Load config file this_year = int(args.year) @@ -39,7 +40,7 @@ try: cfg = yaml.safe_load(stream) except yaml.YAMLError as exc: - print(exc) + logger.info(exc) if not os.path.exists(cfg["out_path"]): os.makedirs(cfg["out_path"]) @@ -64,7 +65,7 @@ ], axis=0, ) -print(hvs) +logger.info(hvs) insts = [] # Load the data @@ -75,7 +76,7 @@ # read the data for c, i in enumerate(hvs): - print(i) + logger.info(i) # Note: There is no need to convert MODIS HDF4 into Netcdf files. You can also use HDF4 files directly. file_collections = glob.glob( @@ -92,7 +93,7 @@ n_cpus=args.n_cpus, ) for c0, fpath in enumerate(file_collections): - print(fpath) + logger.info(fpath) if cfg["satellite"] == "modis": handler = modis(sat_image_path=fpath) handler.load() @@ -100,7 +101,7 @@ handler = VIIRS(sat_image_path=fpath) handler.load() else: - print("Set the satellite in the cfg either to modis or viirs.") + logger.info("Set the satellite in the cfg either to modis or viirs.") # In order to save memory crop the satellite images to our WRF out grid if c0 == 0: @@ -165,14 +166,14 @@ if len(insts) > 1: vprm_inst.add_vprm_insts(insts[1:]) -print(vprm_inst.sat_imgs.sat_img) +logger.info(vprm_inst.sat_imgs.sat_img) # Add the land cover map if not os.path.exists(cfg["out_path"]): os.makedirs(cfg["out_path"]) # Add the land cover map and perform regridding to the satellite grid -print("Generate land cover map") +logger.info("Generate land cover map") veg_file = os.path.join( cfg["out_path"], "veg_map_on_modis_grid_{}_{}.nc".format(args.chunk_x, args.chunk_y) ) @@ -188,13 +189,13 @@ if copernicus_data_path is not None: tiles_to_add = [] for i, c in enumerate(glob.glob(os.path.join(copernicus_data_path, "*"))): - print(c) + logger.info(c) temp_map = copernicus_land_cover_map(c) temp_map.load() dj = vprm_inst.is_disjoint(temp_map) if dj: continue - print("Add {}".format(c)) + logger.info("Add {}".format(c)) if handler_lt is None: handler_lt = temp_map else: @@ -217,7 +218,7 @@ ) # Use all the information in the VPRM instance to generate the WRF input files -print("Create regridder") +logger.info("Create regridder") wrf_op = vprm_inst.to_wrf_output( out_grid, driver="xEMSF", # currently only xESMF works here when the WRF grid is 2D @@ -252,4 +253,4 @@ wrf_op[key].to_netcdf(ofile) -print("Done. In order to inspect the output use evaluate_wrf_input.ipynb") +logger.info("Done. In order to inspect the output use evaluate_wrf_input.ipynb")