PCRGlob model does not like forcings generated with eWatercycle package

[sverhoeven]

I ran

    import pandas as pd
    import ewatercycle.analysis
    import ewatercycle.forcing
    import ewatercycle.models
    import ewatercycle.observation.grdc
    import ewatercycle.parameter_sets

    parameter_set = ewatercycle.parameter_sets.get_parameter_set('pcrglobwb_example_case')

    forcing = ewatercycle.forcing.generate(
        target_model='pcrglobwb',
        dataset='ERA5',
        start_time='1992-01-01T00:00:00Z',
        end_time='1992-12-31T00:00:00Z',
        shape='Meuse/Meuse.shp',
        model_specific_options=dict(
            start_time_climatology='1991-01-01T00:00:00Z',
            end_time_climatology='1991-12-31T00:00:00Z',
        )
    )

    model = ewatercycle.models.PCRGlobWB(version="setters", parameter_set=parameter_set, forcing=forcing)

    cfg_file, cfg_dir = model.setup(max_spinups_in_years=1)

    model.initialize(cfg_file)

    observations_df, station_info = ewatercycle.observation.grdc.get_grdc_data(
        station_id=6421500,
        start_time=model.start_time_as_isostr,
        end_time=model.end_time_as_isostr,
    )
    station_lon = station_info['grdc_longitude_in_arc_degree']
    station_lat = station_info['grdc_latitude_in_arc_degree']

    simulated_discharge = []
    timestamps = []
    while (model.time < model.end_time):
        model.update()
        discharge = model.get_value_at_coords('discharge', lat=[station_lat], lon=[station_lon])
        simulated_discharge.append(discharge)
        timestamps.append(model.time_as_datetime)

But on call to model.update() it raised following exception

_InactiveRpcError: <_InactiveRpcError of RPC that terminated with:
	status = StatusCode.UNKNOWN
	details = "Exception calling application: 'precipitation'"
	debug_error_string = "{"created":"@1626168569.168831104","description":"Error received from peer ipv6:[::1]:44483","file":"src/core/lib/surface/call.cc","file_line":1066,"grpc_message":"Exception calling application: 'precipitation'","grpc_status":2}"

The log file showed:

2021-07-13 09:25:21,960 bmiPcrglobwb INFO Shape of maps is (13, 17)
2021-07-13 09:25:21,960 bmiPcrglobwb INFO PCRGlobWB Initialized
2021-07-13 09:29:29,133 pcrglobwb INFO Reading forcings for time 1992-01-01
2021-07-13 09:29:29,167 grpc._server ERROR Exception calling application: 'precipitation'
Traceback (most recent call last):
  File "/usr/local/lib/python2.7/dist-packages/grpc/_server.py", line 435, in _call_behavior
    response_or_iterator = behavior(argument, context)
  File "/usr/local/lib/python2.7/dist-packages/grpc4bmi/bmi_grpc_server.py", line 31, in update
    self.bmi_model_.update()
  File "/opt/PCR-GLOBWB_model/model/bmiPcrglobwb.py", line 211, in update
    self.model.read_forcings()
  File "/opt/PCR-GLOBWB_model/model/pcrglobwb.py", line 439, in read_forcings
    self.meteo.read_forcings(self._modelTime)
  File "/opt/PCR-GLOBWB_model/model/meteo.py", line 388, in read_forcings
    LatitudeLongitude = True)
  File "/opt/PCR-GLOBWB_model/model/virtualOS.py", line 365, in netcdf2PCRobjClone
    cropData = f.variables[varName][int(idx),:,:]       # still original data
KeyError: 'precipitation'

[sverhoeven]

The header of generated pr forcing is

ncdump -h /home/verhoes/git/eWaterCycle/ewatercycle/docs/examples/esmvaltool_output/recipe_pcrglobwb_20210713_090147/work/diagnostic_daily/script/pcrglobwb_OBS6_ERA5_reanaly_1_day_pr_1992-1992_Meuse.nc
netcdf pcrglobwb_OBS6_ERA5_reanaly_1_day_pr_1992-1992_Meuse {
dimensions:
	time = 731 ;
	lat = 16 ;
	lon = 12 ;
	bnds = 2 ;
variables:
	float pr(time, lat, lon) ;
		pr:_FillValue = 1.e+20f ;
		pr:standard_name = "precipitation_flux" ;
		pr:long_name = "Precipitation" ;
		pr:units = "m" ;

It is a bit weird that file contains 2 years, but file name only single year.

[Peter9192]

So the problem occured when it reached the second year?

[sverhoeven]

Problem occurred in first call to model.update() as timestamps and simulated_discharge vars are still empty.

[sverhoeven]

In the https://github.com/UU-Hydro/PCR-GLOBWB_input_example/blob/master/RhineMeuse30min/forcing/precipitation_2001to2010.nc the var is

	float precipitation(time, lat, lon) ;
		precipitation:_FillValue = 1.e+20f ;
		precipitation:standard_name = "precipitation" ;
		precipitation:long_name = "precipitation" ;
		precipitation:units = "undefined" ;

This is different from the nc file generated by ESMValTool, maybe the model config needs adjustment or recipe needs help?

[SarahAlidoost]

On cartesis /projects/0/wtrcycle/comparison/forcing/pcrglobwb/, there is a forcing file pcrglobwb_OBS6_ERA5_reanaly_1_day_pr_1990-1990_rhine.nc with this info:

netcdf pcrglobwb_OBS6_ERA5_reanaly_1_day_pr_1990-1990_rhine {
dimensions:
        time = 730 ;
        lat = 49 ;
        lon = 55 ;
        bnds = 2 ;
variables:
        float pr(time, lat, lon) ;
                pr:_FillValue = 1.e+20f ;
                pr:standard_name = "precipitation_flux" ;
                pr:long_name = "Precipitation" ;
                pr:units = "m" ;
        float time(time) ;
                time:axis = "T" ;
                time:bounds = "time_bnds" ;
                time:units = "days since 1850-1-1 00:00:00" ;
                time:standard_name = "time" ;
                time:long_name = "time" ;
                time:calendar = "gregorian" ;
        float time_bnds(time, bnds) ;
        float lat(lat) ;
                lat:axis = "Y" ;
                lat:bounds = "lat_bnds" ;
                lat:units = "degrees_north" ;
                lat:standard_name = "latitude" ;
                lat:long_name = "Latitude" ;
        float lat_bnds(lat, bnds) ;
        float lon(lon) ;
                lon:axis = "X" ;
                lon:bounds = "lon_bnds" ;
                lon:units = "degrees_east" ;
                lon:standard_name = "longitude" ;
                lon:long_name = "Longitude" ;
        float lon_bnds(lon, bnds) ;

// global attributes:
                :comment = "Contains modified Copernicus Climate Change Service Information 2020" ;
                :Conventions = "CF-1.7" ;
                :provenance

This file can be re-generated using the example notebook generate_forcing.ipynb. Here its info:

netcdf pcrglobwb_OBS6_ERA5_reanaly_1_day_pr_1990-1990_rhine_ewater {
dimensions:
        time = 730 ;
        lat = 23 ;
        lon = 31 ;
        bnds = 2 ;
variables:
        float pr(time, lat, lon) ;
                pr:_FillValue = 1.e+20f ;
                pr:standard_name = "precipitation_flux" ;
                pr:long_name = "Precipitation" ;
                pr:units = "m" ;
        float time(time) ;
                time:axis = "T" ;
                time:bounds = "time_bnds" ;
                time:units = "days since 1850-1-1 00:00:00" ;
                time:standard_name = "time" ;
                time:long_name = "time" ;
                time:calendar = "gregorian" ;
        float time_bnds(time, bnds) ;
        float lat(lat) ;
                lat:axis = "Y" ;
                lat:bounds = "lat_bnds" ;
                lat:units = "degrees_north" ;
                lat:standard_name = "latitude" ;
                lat:long_name = "Latitude" ;
        float lat_bnds(lat, bnds) ;
        float lon(lon) ;
                lon:axis = "X" ;
                lon:bounds = "lon_bnds" ;
                lon:units = "degrees_east" ;
                lon:standard_name = "longitude" ;
                lon:long_name = "Longitude" ;
        float lon_bnds(lon, bnds) ;

// global attributes:
                :comment = "Contains modified Copernicus Climate Change Service Information 2020" ;
                :Conventions = "CF-1.7" ;
                :provenance 

These two files are different only in lat/long dimensions. The reason is our get_extents function. So, I think this issue is related to the pcrglob config file rather than forcing data.

[SarahAlidoost]

I compared two pcrglob config files: one from Cartesius, and the other one from pcrglobwb_example_case. The one on the cartesius has these two items in the section [meteoOptions]:

precipitationVariableName = pr
temperatureVariableName = tas

I added these two items to the config file generated by ewatercycle and got a different error:

_InactiveRpcError: <_InactiveRpcError of RPC that terminated with:
	status = StatusCode.UNKNOWN
	details = "Exception calling application: Size mismatch: Number of array elements is 713, number of raster cells is 221"
	debug_error_string = "{"created":"@1626190527.745444205","description":"Error received from peer ipv4:127.0.0.1:50039","file":"src/core/lib/surface/call.cc","file_line":1067,"grpc_message":"Exception calling application: Size mismatch: Number of array elements is 713, number of raster cells is 221","grpc_status":2}"
>

[sverhoeven]

1. So the ESMValTool generates pr var while the default model value is precipitation. We could add precipitationVariableName and temperatureVariableName props to PCRGlobWBForcing and in PCRGlobWB set those fields in the config in _setup_default_config() method. Or alternatively we could add them to the PCRGlobWB.setup() method.

2. The other nc files in example parameter set have

 lat = 52.25, 51.75, 51.25, 50.75, 50.25, 49.75, 49.25, 48.75, 48.25, 47.75, 
    47.25, 46.75, 46.25 ;

 lon = 3.75, 4.25, 4.75, 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 
    9.25, 9.75, 10.25, 10.75, 11.25, 11.75 ;

While the generated forcing has

 lat = 48, 48.25, 48.5, 48.75, 49, 49.25, 49.5, 49.75, 50, 50.25, 50.5, 
    50.75, 51, 51.25, 51.5, 51.75 ;
 lon = 3.75, 4, 4.25, 4.5, 4.75, 5, 5.25, 5.5, 5.75, 6, 6.25, 6.5 ;

Now I used a random shape file, which the extents are derived from. Seems the grid positions of forcing and parameter set need to be the same? I could set the extract_region while generating forcing to same extend as parameter set.

[SarahAlidoost]

pcrglob uses a clone map. In era5-comparison, we make the clone maps for each catchment, please see cell Check the clone extents and make clone maps and Load the recipe file, change the settings and save it in this preprocessing notebook. As shown, the extents of forcing are derived from clone maps. These two steps are not done in eWaterCycle generate forcing.

[Peter9192]

So either we

• also generate the clone map (based on a shapefile) and make it part of the PCRGlobWB Forcing
• derive the bounds from an existing clonemap, which would have to be passed as an additional parameter

First option would be more in line with our philisophy. Second one would help more for the example case. Maybe we can even support both...

We might also want to make the grid resolution configurable.

[Peter9192]

pcrglob uses a clone map. In era5-comparison, we make the clone maps for each catchment, please see cell Check the clone extents and make clone maps and Load the recipe file, change the settings and save it in this preprocessing notebook. As shown, the extents of forcing are derived from clone maps. These two steps are not done in eWaterCycle generate forcing.

So in that notebook, indeed, we derive the entire clonemap. But for generating the forcing, we eventually just pass the extents.
The important note here is that

"The clonemap extent divided by the forcing resolution must yield an integer number of grid cells."

so given the forcing resolution of ERA5 (.25 deg) it seems that setting the bounds to that of the example clonemap should produce compatible forcing data.

[Peter9192]

1. add `precipitationVariableName` and `temperatureVariableName` props to PCRGlobWBForcing and in `PCRGlobWB` set those fields in the config in `_setup_default_config()` method.

Did that in eWaterCycle/ewatercycle#179.

Seems the grid positions of forcing and parameter set need to be the same? I could set the extract_region while generating forcing to same extend as parameter set.

I tried setting the extents right, but since the data resolution is higher than that of the clonemap, it still doesn't work. We need to add a regrid processor to the ESMValTool recipe, so it's not gonna have to wait for another ESMValTool release for this to work out of the box. Or we'd have to change the 'extract_region' preprocessor to the latest regional regrid from inside the ewatercycle.forcing module.

[Peter9192]

Or we'd have to change the 'extract_region' preprocessor to the latest regional regrid from inside the ewatercycle.forcing module.

Tried that but failed, see ESMValGroup/ESMValCore#1230

[Peter9192]

I was able to run the example case with forcing generated by ESMValTool, using the new regional regrid option with the bugfix branch from ESMValGroup/ESMValCore#1231.

Due to ESMValGroup/ESMValCore#1232 I couldn't use the API yet, so instead I edited the recipe obtained with esmvaltool recipes get hydrology/recipe_pcrglobwb.yml.

Specifically, I replaced the extract_region preprocessor with the following regrid spec:

  crop_basin: &crop_basin
    regrid:
      scheme: linear
      lon_offset: true
      lat_offset: true
      target_grid:
        start_longitude: 3.75
        end_longitude: 11.75
        start_latitude: 46.25
        end_latitude: 52.25
        step_longitude: 0.5
        step_latitude: 0.5

I'll push another commit to eWaterCycle/ewatercycle#179 that enables using the regridder once it's available through the API.

[Peter9192]

Even though it now doesn't complain anymore, the output doesn't look right though...

[Peter9192]

with a new parameterset at higher resolution we need to make it possible to regrid anyway.