fmaussion / salem

Add geolocalised subsetting, masking, and plotting operations to xarray
http://salem.readthedocs.io
Other
186 stars 44 forks source link

Reproject WRF out file to WGS84 to match ERA5 data #232

Closed krober10nd closed 1 year ago

krober10nd commented 1 year ago

Hi,

Thanks for the great package. I have a seemingly easy task that has become quite tricky. I want to compute a correlation map between a WRF run and the corresponding ERA5 data file I used to force the model.

To compute the correlation both datasets need to be in the same CRS which is why I'm doing what I'm doing.

Here's a minimum failing example reprojecting the wrf file to WGS84 to match ERA5 data:

import Salem 

era5_file = "ERA5_2018070800.nc"
wrf_file = "wrfout_d01_2018-07-08_00/00/00"
era5_data = salem.open_xr_dataset(era5_file)
# make into -180 to 180 deg
era5_data = era5_data.assign_coords(longitude=(((era5_data.longitude + 180) % 360) - 180)).sortby('longitude')
wrf_data = salem.open_wrf_dataset(wrf_files[0])
# Reproject WRF data from Mercator to WGS84
wrf_data_reproj = wrf_data.salem.transform(wrf_data)

The last line of the above example fails with

RuntimeError: dataset Grid not understood

despite both datasets containing a Salem.Grid

Data is attached DataForExample.zip

Thanks for any help!

fmaussion commented 1 year ago

This line here does not make sense?

wrf_data_reproj = wrf_data.salem.transform(wrf_data)

This attempts to reproject wrf onto itself?

krober10nd commented 1 year ago

Just a typo on my part. It should read 'era5_data'. Regardless it doesn't work either way.

fmaussion commented 1 year ago

Replace the last line with

# Reproject WRF data from Mercator to WGS84
T2_data_reproj = era5_data.salem.transform(wrf_data.T2, grid=wrf_data.salem.grid)
T2_data_reproj.salem.quick_map();

The problem occurs because you are trying to reproject all variables in the wrf file, some of them cannot be reprojected because they don't have a coordinate system (e.g. the pressure coordinates). salem could be a bit cleverer and not try to reproject this, but I would recommend to only reproject the data you really need anyway for performance reasons.

krober10nd commented 1 year ago

Thanks @fmaussion. I see, restrict to only one variable at a time.

Interestingly, it plots with latitude and longitude ticks with the quick_map but the pyproj_srs in the supposedly reprojected DataArray remains in Mercator and not in EPSG:4326.

The same remains for the inverse transform Mercator to EPSG:4326.

Edit: I also don't see any change to the salem.grid object nor the coordinates .coords in the reprojected dataset.

fmaussion commented 1 year ago

restrict to only one variable at a time.

Or select the ones you really need to reproject, e.g.:

ds_sel = ds[['var1', 'var2', ...]]
# reproject ds_sel

For the rest I'm not sure to understand what you mean, sorry! Please provide a MWE and clearly state what the problem is.

krober10nd commented 1 year ago

era5_file = "ERA5_2018070800.nc"
wrf_file = "wrfout_d01_2018-07-08_00/00/00"
era5_data = salem.open_xr_dataset(era5_file)
# make into -180 to 180 deg
era5_data = era5_data.assign_coords(longitude=(((era5_data.longitude + 180) % 360) - 180)).sortby('longitude')
era5_u10 = era5_data['u10']

wrf_data = salem.open_wrf_dataset(wrf_files[0])
wrf_u10 = wrf_data['U10']
# Reproject WRF data from Mercator to WGS84
wrf_u10_reproj = wrf_u10.salem.transform(era5_u10)

wrf_u10_reproj.salem.quick_map() # works appears in lat/lon 

print(wrf_u10_reproj.salem.pyproj_srs) # appears in mercator and coordinate vectors are in mercator. 
fmaussion commented 1 year ago
wrf_u10.salem.transform(wrf_u10)

here you are reprojecting unto oneself...

krober10nd commented 1 year ago

typo. fixed above. same problem. I have to say this is remarkably difficult to manipulate these files.

fmaussion commented 1 year ago

This code here:

import salem
import xarray as xr

era5_file = "ERA5_2018070800.nc"
wrf_file = "WRFout/wrfout_d01_2018-07-08_00:00:00"
era5_data = salem.open_xr_dataset(era5_file)

# make into -180 to 180 deg
era5_data = era5_data.assign_coords(longitude=(((era5_data.longitude + 180) % 360) - 180)).sortby('longitude')

wrf_data = salem.open_wrf_dataset(wrf_file)

# Reproject WRF data from Mercator to WGS84
T2_data_reproj = era5_data.salem.transform(wrf_data.T2, grid=wrf_data.salem.grid)
T2_data_reproj.salem.grid

Returns

<salem.Grid>
  proj: +datum=WGS84 +no_defs+proj=longlat
  pixel_ref: center
  origin: upper-left
  (nx, ny): (75, 69)
  (dx, dy): (0.5, -0.5)
  (x0, y0): (-102.0, 47.0)

Which is working as expected.

I have to say this is remarkably difficult to manipulate these files.

Yes, I think you should have a look at another library than salem. I recommend a few on the doc: https://salem.readthedocs.io/en/stable/faq.html#what-others-tools-should-i-know-about

krober10nd commented 1 year ago

Thanks Fabien. I've tried all of them (ncl, cdo, wrf-python, xwrf). My recommendation was to use a different model than WRF.

I wanted to evaluate its performance of a simulation I performed but it was almost impossible to extract a point time series and then there's the headache with installation.