update grid in nc format gives xarray error

Devaraj-G commented 3 years ago

Running the following command: remove-restore --base ./Work/Data_Code/IAAC2/2021_07_06_gebco_2020/gebco_2020_n60.0_s37.0_w-140.0_e-119.0.nc --output ./Work/Data_Code/IAAC2/temp/merged_dem.nc --spacing 0.00099 --region_of_interest -137.250000000 -121.250000000 39.750000000 56.250000000 ./Work/Data_Code/IAAC2/2021_07_05_noaa_Strait_of_Juan_de_Fuca_1_by_3_arcsecond/strait_of_juan_de_fuca_13_navd88_2015.nc

The nc file is available at Strait of Juan de Fuca 1/3 arc-second NAVD 88 Coastal Digital Elevation Model

gives: _Loading update grid Loading ./Work/Data_Code/IAAC2/2021_07_05_noaa_Strait_of_Juan_de_Fuca_1_by_3_arcsecond/strait_of_juan_de_fuca_13_navd88_2015.nc Traceback (most recent call last): File "/home/devaraj/anaconda3/envs/cascadia/bin/remove-restore", line 8, in sys.exit(main()) File "/home/devaraj/anaconda3/envs/cascadia/lib/python3.9/site-packages/pycascadia/remove_restore.py", line 119, in main update_grid = Grid(fname, convert_to_xyz=True) File "/home/devaraj/anaconda3/envs/cascadia/lib/python3.9/site-packages/pycascadia/grid.py", line 13, in init self.load(fname) File "/home/devaraj/anaconda3/envs/cascadia/lib/python3.9/site-packages/pycascadia/grid.py", line 20, in load self.grid, self.region, self.spacing = load_source(fname) File "/home/devaraj/anaconda3/envs/cascadia/lib/python3.9/site-packages/pycascadia/loaders.py", line 17, in load_source xr_data = load_netcdf(filepath) File "/home/devaraj/anaconda3/envs/cascadia/lib/python3.9/site-packages/pycascadia/loaders.py", line 70, in load_netcdf xr_data = xr.open_dataarray(filepath).astype('float32') File "/home/devaraj/anaconda3/envs/cascadia/lib/python3.9/site-packages/xarray/backends/api.py", line 671, in open_dataarray raise ValueError( ValueError: Given file dataset contains more than one data variable. Please read with xarray.opendataset and then select the variable you want.

Devaraj-G commented 3 years ago

Maybe because _NETCDFVARNAME=Band1 instead of _NETCDFVARNAME=elevation?

JamieJQuinn commented 3 years ago

So this error is because xarray (and remove-restore) expects a netcdf file with a single variable. This particular file contains two, the actual data in Band1 and an additional variable crs:

<xarray.Dataset>
Dimensions:  (lat: 7345, lon: 17281)
Coordinates:
  * lat      (lat) float64 48.04 48.04 48.04 48.04 ... 48.72 48.72 48.72 48.72
  * lon      (lon) float64 -125.1 -125.1 -125.1 -125.1 ... -123.5 -123.5 -123.5
Data variables:
    crs      |S1 b''
    Band1    (lat, lon) float32 -138.3 -138.2 -138.2 ... -68.53 -68.4 -68.29

This could be automatically handled by remove-restore but it wouldn't generalise to other files (i.e. it would be a very specific bit of logic to handle this particular file). I think a better way is to preprocess this file, removing the crs variable from the netcdf dataset. @Devaraj-G do you know if gdal can do this? Alternatively, it can be done using xarray:

import xarray
in_fname = "strait_of_juan_de_fuca_13_navd88_2015.nc"
out_fname = "strait_of_juan_de_fuca_13_navd88_2015_crs_removed.nc"
ds = xarray.load_dataset(in_fname)
del ds['crs']
ds.to_netcdf(out_fname)

Devaraj-G commented 3 years ago

gdalinfo only shows Band 1. How about gda_translate with option -b specified as 1.

JamieJQuinn commented 3 years ago

This will not be handled specifically by remove-restore. The input netcdf file should contain only one variable.

UCL / pyCascadia

update grid in nc format gives xarray error #59