usgs-r/drb_gridmet_tools

Repository with functions to aggregate raster data from pixel grid to hru polygon grid. This repository heavily relies on grd2shp_xagg library.

drb_gridmet_tools's Issues

Add note about OS requirements

I just tried to build on Windows and it seems like the esmpy package is not available for Windows using the channels in environment.yml. Recommending adding a note to the README about OS requirements, or finding the Windows channel for esmpy.

1 Dimensional error when running regridding

Getting the following error when running g2s.run_weights() on with the polygon file nhru_02.shp

ValueError: shape mismatch: value array of shape (1, 15342) could not be broadcast to indexing result with 1 dimensions

Overarching function found gridmet_regridding() - error occuring at L:152

Extract Gridmet data to points of interest

Precip aggregation to PRMS scale: Area-weighted average, area-weighted sum, vs sum

Recap:

8 gridmet variables were extracted for the drb. These gridmet metrics are daily timesteps since 1979 and available for each hru.

Since we are working with the PRMS segment network, we choose to match the PRMS segments with adjacent HRU polygons and aggregate the values of the Hrus to PRMS.

For temp ('tmmx', 'tmmn'), radiation (srad), windspeed (vs), and humidity measurements (rmin, rmax, sph), the aggregation is done as a area-weighted average such that the size of the hru polygon influences the mean value tied to that segment.

Our final gridmet variable pr, is accumulated precipitaiton (inch). Should the area of the hru affect the precip aggregation? We could compute:

• area-weighted mean like the other variables - but this would not give a strong indication of the total amount of potential run off into the segment on a given day.
• sum - but possible bias present for prms segments tied to many hrus.

Example:

Take for example, PRMS segment 886_1 (blue) which has two HRUs on each side of the segment:

Day: 1885-09-09

The area-weighted average gives us a value of 0.1814 inch non-weighted mean = 0.185 inch)

(Note - i merged the hru geometry information to the aggregated table for this image).

Meanwhile the basic sum of accumulated precip for 886_1 = ~0.37 inch.

Please provide your preference on this @jds485 @jzwart

Error writing g2s to netcdf on selected data variables

Function in gridmet.py extracted the gridmet data for selected variables and run the regridding process from the grd2shp_xagg package

Function inputs for inland salinity are:

polygon_file_path = ' nhru_02.shp',
var_short = ['tmmx', 'tmmn', 'pr', 'srad', 'vs','rmax','rmin','sph'] ,
start_date='1979-01-01', 
end_date='2021-01-01',  # set right now to 1989-01-01 to test script and speed up processing
lat_max = lat_max,  lat_min = lat_min , lon_max = lon_max, lon_min = lon_min

Where var_short are defined as:

{'tmmx': 'daily_maximum_temperature',
 'tmmn': 'daily_minimum_temperature',
 'pr': 'precipitation_amount',
 'srad': 'daily_mean_shortwave_radiation_at_surface',
 'vs': 'daily_mean_wind_speed',
 'rmax': 'daily_maximum_relative_humidity',
 'rmin': 'daily_minimum_relative_humidity',
 'sph': 'daily_mean_specific_humidity'}

and lat lon vars are defined by the ~ bbox of the drb

• Comments for review:
  • Getting Various errors: timeout error when running weights on grd2shp_xagg obj, dimension error the last step when writing the g2s object. g2s.write_gm_file(opath=output_data_folder, prefix= g2s_file_prefix)
  • For the initializing of grid2shpxagg object (grd2shp_xagg.Grd2ShpXagg().initialize), the time_var = 'day'. Do we want to keep our datasets at daily time intervals?
  • Process takes quite a long time. Any idea how to speed up process?
  • Currently one function for entire regridding. This can be split into 3 different functions if we prefer.

Lmk your thoughts.

Issue reviewers:
Code review: @jsadler2
Comment review: @jds485

Convert gridmet processing into snakemake pipeline

currently, the processing occurs in python scripts. It would be helpful to have it as a single snakemake pipeline.

Processing time for `np.agg()` long

Processing time for running group by on gridmet data, already converted to xarray ds is very long when using np.agg() such as:

gridmet_drb_gdf..groupby(['PRMS_segid',"time"]).agg(
    area = ("hru_area_m2", "sum"),
    pr = ("pr", 'sum'),
    tmmx = ("tmmx", area_weighted_avg),
    tmmn = ("tmmn", area_weighted_avg),
    srad = ('srad', area_weighted_avg),
    vs = ('vs', area_weighted_avg),
    rmax = ('rmax',area_weighted_avg),
    rmin = ('rmin', area_weighted_avg),
    sph = ('sph',  area_weighted_avg)
)

Error regridding for precipitation amount

Getting an error when regridding the precipitation data.

1979-01-01 2021-01-01
numtimesteps: 15342 and Start date: 1979-01-01
generating mapped vales for daily_maximum_temperature ...
aggregating daily_maximum_temperature...
all variables aggregated to polygons!
finished mapped values for daily_maximum_temperature
generating mapped vales for daily_minimum_temperature ...
aggregating daily_minimum_temperature...
all variables aggregated to polygons!
finished mapped values for daily_minimum_temperature
generating mapped vales for precipitation_amount ...
Warning:curl error: Transferred a partial file; ignored
Error:DAP DATADDS packet is apparently too short
aggregating precipitation_amount...
/Users/msleckman/miniconda3/envs/gridmet_tools_IS_1/lib/python3.10/site-packages/xagg/core.py:454: UserWarning: One or more of the pixels in variable precipitation_amount have nans in them in the dimensions day. The code can currently only deal with pixels for which the *entire* pixel has nan values in all dimensions, however there is currently no  support for data in which pixels have only some nan values. The aggregation calculation is likely incorrect.
  warnings.warn('One or more of the pixels in variable '+var+' have nans in them in the dimensions '+
/Users/msleckman/miniconda3/envs/gridmet_tools_IS_1/lib/python3.10/site-packages/xagg/core.py:454: UserWarning: One or more of the pixels in variable precipitation_amount have nans in them in the dimensions day. The code can currently only deal with pixels for which the *entire* pixel has nan values in all dimensions, however there is currently no  support for data in which pixels have only some nan values. The aggregation calculation is likely incorrect.
  warnings.warn('One or more of the pixels in variable '+var+' have nans in them in the dimensions '+
/Users/msleckman/miniconda3/envs/gridmet_tools_IS_1/lib/python3.10/site-packages/xagg/core.py:454: UserWarning: One or more of the pixels in variable precipitation_amount have nans in them in the dimensions day. The code can currently only deal with pixels for which the *entire* pixel has nan values in all dimensions, however there is currently no  support for data in which pixels have only some nan values. The aggregation calculation is likely incorrect.
  warnings.warn('One or more of the pixels in variable '+var+' have nans in them in the dimensions '+
/Users/msleckman/miniconda3/envs/gridmet_tools_IS_1/lib/python3.10/site-packages/xagg/core.py:454: UserWarning: One or more of the pixels in variable precipitation_amount have nans in them in the dimensions day. The code can currently only deal with pixels for which the *entire* pixel has nan values in all dimensions, however there is currently no  support for data in which pixels have only some nan values. The aggregation calculation is likely incorrect.
  warnings.warn('One or more of the pixels in variable '+var+' have nans in them in the dimensions '+
/Users/msleckman/miniconda3/envs/gridmet_tools_IS_1/lib/python3.10/site-packages/xagg/core.py:454: UserWarning: One or more of the pixels in variable precipitation_amount have nans in them in the dimensions day. The code can currently only deal with pixels for which the *entire* pixel has nan values in all dimensions, however there is currently no  support for data in which pixels have only some nan values. The aggregation calculation is likely incorrect.
  warnings.warn('One or more of the pixels in variable '+var+' have nans in them in the dimensions '+
/Users/msleckman/miniconda3/envs/gridmet_tools_IS_1/lib/python3.10/site-packages/xagg/core.py:454: UserWarning: One or more of the pixels in variable precipitation_amount have nans in them in the dimensions day. The code can currently only deal with pixels for which the *entire* pixel has nan values in all dimensions, however there is currently no  support for data in which pixels have only some nan values. The aggregation calculation is likely incorrect.
  warnings.warn('One or more of the pixels in variable '+var+' have nans in them in the dimensions '+
all variables aggregated to polygons!
finished mapped values for precipitation_amount

No nan directly visible in xarray output:

xarray_dict['pr']
Out[10]: 
<xarray.Dataset>
Dimensions:               (lat: 210, crs: 1, lon: 208, day: 15342)
Coordinates:
  * lat                   (lat) float64 45.36 45.32 45.28 ... 36.73 36.69 36.65
  * crs                   (crs) float32 3.0
  * lon                   (lon) float64 -80.52 -80.47 -80.43 ... -71.93 -71.89
  * day                   (day) datetime64[ns] 1979-01-01 ... 2021-01-01
Data variables:
    precipitation_amount  (day, lat, lon) float32 ...
Attributes: (12/19)
    geospatial_bounds_crs:      EPSG:4326
    Conventions:                CF-1.0
    geospatial_bounds:          POLYGON((-124.7666666333333 49.40000000000000...
    geospatial_lat_min:         25.066666666666666
    geospatial_lat_max:         49.40000000000000
    geospatial_lon_min:         -124.7666666333333
    ...                         ...
    date:                       25 January 2022
    note1:                      The projection information for this file is: ...
    note2:                      Citation: Abatzoglou, J.T., 2013, Development...
    note3:                      Data in slices after last_permanent_slice (1-...
    note4:                      Data in slices after last_provisional_slice (...
    note5:                      Days correspond approximately to calendar day...

Graphically, we see No data after a specific 6000 days since 01-01-1979.
green: geomid = 2, orange : geomid = 3, blue : geomid =1

This is the only var with missing data atm.

• Will test is new env

Build Container for gridmet data extraction process

Extract gridmet data along segments polylines

• Approach 1:

Create small buffer around polyline and run regridding as a polygon.

• Approach 2:

aggregate to NHD catchments

For the PGDL water temperature sub-quest to account for GW influence, the gridmet data will need to be aggregated to the nhd network. The DO project might also use this more fine-scale data.

I think the only things we'll need to change are:

1. make available a new geopackage with the nhd catchment polygons for the DRB. I believe it makes most sense for this to come from "https://github.com/USGS-R/drb-network-prep/"
2. Revise the Snakemake code to accept the layer name as variable instead of hard coded

   drb_gridmet_tools/Snakefile

   Line 34 in d2a4daf

   gdf = gpd.read_file(input[0], layer="GFv1_catchments_edited")

3. Probably change the "prefix" value to say something about "nhd" so we can distinguish the output

   drb_gridmet_tools/Snakefile

   Line 73 in d2a4daf

   g2s_file_prefix='drb_',

4. Change the value on this line to the nhd catchment geopackage

   drb_gridmet_tools/Snakefile

   Line 9 in d2a4daf

   drb_hru_file = "data/GFv1_catchments_edited.gpkg"

cc @msleckman @lekoenig