fix wrong reprojction to wrf domain

[matthiasdemuzere]

Addressing Issue #56

Issue: Reprojection from the high-res LCZ grid to the lower-res WRF grid resulted in errors for larger domains.
Reason: WRF has an irregular latlong grid, and this was not taken into account when reprojecting.
Solution: Inspired by the blog information provided by Fabien Maussion here, the procedure is not completely changed:

• the create_wrf_gridinfo() routine that wrote a temporary .tif file, containing WRF's grid info, is now removed.
• instead a new routine _get_wrf_grid_info() is introduced to get a proper description of the WRF grid, in its native Lambert Conformal Conic projection.
• this routine provide the proper crs and transform values needed to Reproject the lcz info to the wrf grid.

[matthiasdemuzere]

Ok, as far as I can see this all works now:

• new reprojection procedure seems to work well.
• tests have been updated and pass
• new test has been added to check the _get_wrf_grid_info() routine
• relevant sample_data/ and testing/ files have been updated

@dargueso @theendlessriver13: perhaps you can review?

If accepted, I guess we should bump the version number, and update pypi.
Especially for the larger WRF domains this change will have some impact ....

[dargueso]

@matthiasdemuzere, sure, I can check it out, but won't be able to do so before Monday.
Hope that works.

[matthiasdemuzere]

I realized there is another level of complication that we did not yet take care of.

Geo_em files can be on different projections. This is information is stored in MAP_PROJ.

for sample data:

>>> ds = xr.open_dataset(geo_em.d04.nc)
>>> print(f'{case}: MAP_PROJ = {ds.MAP_PROJ}')
sample_data: MAP_PROJ = 6

for Italy data from pippo2013:

>>> ds = xr.open_dataset(geo_em.d04.nc)
>>> print(f'{case}: MAP_PROJ = {ds.MAP_PROJ}')
sample_data: MAP_PROJ = 1

MAP_PROJ can have the following values:

• 1: Lambert,
• 2: polar stereographic,
• 3: mercator,
• 6: lat-lon

More info:

• https://www2.mmm.ucar.edu/wrf/users/namelist_best_prac_wps.html#map_proj
• http://140.112.69.65/research/coawst/COAWST_TUTORIAL/training_2019/monday/werner_wps.pdf

So far this new code caters for 1-Lambert only.

So a check is needed within _get_wrf_grid_info(), to:

• get the MAP_PROJ value
• assign the corresponding wrf_proj string

We can probably get a long way borrowing some code from salem.
And maybe this issue as well.

I can only look into this later, no time at the moment ...

[dargueso]

@matthiasdemuzere, sorry, I probably wasn't clear enough in my last messages. I thought you were aware of that. When I discussed the issue (#56 (comment)), I had this multiple projection option in mind. You have to add the rotated-pole option to the above.

I think we can indeed do most of the work with salem or wrf-python. WRF-python reads all possible wrf projections and makes it easy to reproject after.

[matthiasdemuzere]

So far I developed a solution without additional packages. Could be good to keep it this way, as e.g. salem is somewhat maintained but not further developed ...

Let me add these other ones.
For MAP_PROJ = 6, I'll have to check how this should be done.

[matthiasdemuzere]

Ok, I now added the different type of projections.

MAP_PROJ = 1 was tested with the Italy data. That looks fine, let's see if Pippo2013 can confirm this.

MAP_PROJ = 6 was tested with input from
costa_del_sol_input.zip

Unfortunately, the projection does not seem to be correct yet. So this needs more work.
MAP_PROJ = 2 or 3 also have not been tested.

So all of this still needs a bit of work ...

[matthiasdemuzere]

I believe this is done now.

I tested the reprojection implementation for MAP_PROJ 1, 2 and 6. All seems good.

@theendlessriver13 do you want to review this one last time?

If accepted, I guess we should bump the version number, and update pypi.
Especially for the larger WRF domains this change will have some impact ....

[jkittner]

Looking very good, nice that this is done, thank you!
A few final questions:

• can you elaborate a little what has changed in the netcdf files - sample and testing data?
• Why are so many hard coded values in the tests changed, all due to the different projection and slightly moved pixels? What does this mean for data generated with the 0.2.0 version are they all slightly off?
• Would you classify this patch as a "bug fix" or a "new feature" so should this become version 0.2.1 or 0.3.0?

[matthiasdemuzere]

• can you elaborate a little what has changed in the netcdf files - sample and testing data?

The netcdf files in sample_data/ are updated with the correct re-projected for the LCZ classes and urban canopy parameters, using this new implementation. As a result, the LCZ information in LCZ_extent.nc and LCZ_params.nc is slightly different.
The same is true for some of the files in `testing/

• Why are so many hard coded values in the tests changed, all due to the different projection and slightly moved pixels? What does this mean for data generated with the 0.2.0 version are they all slightly off?

The previous reprojection was slightly off, as I was treating the WRF grid as a regular grid, but it is not.
The WRF grid can have 4 projection types (MAP_RPOJ 1, 2, 3 and 6) which are now all properly used when reprojecting the higher-res LCZ data onto the lower-res WRF grid.
That indeed means that the data generated with the 0.2.0 version is slightly off. Impact will be minimal for smaller domains, but more pronounced for bigger domains.

• Would you classify this patch as a "bug fix" or a "new feature" so should this become version 0.2.1 or 0.3.0?

This patch fixes a bug in the reprojection from high-res LCZ to low-res WRF.
Yet since outputs are changing between these two versions, it is probably better to assign this as a new feature, to bump the version from 0.2.1 to 0.3.0?

[jkittner]

Nice, looking very good! Great that you were able to fix this. I'll make a new release to PyPi!