Shift in LCZ substituted values?

[pippo2013]

@matthiasdemuzere thanks a lot for sharing w2w! I am trying to use it in a domain in Italy but I obtain a shift in urban settlements. You can see in the left plot original LU_Index (from geo_em.d03.nc) and in the left one LU_INDEX from geo_em.d03_LCZ_extent.nc New classes are added but they are shifted with respect to original urban data, I get urban settlements in the sea! This is my geo file https://drive.google.com/file/d/1BAvlREj9bJo_sgVQ4EMqs96QAvlhZPug/view?usp=sharing I cannot understand what the problem might be. The LCZ file used for w2w is EU_LCZ_map_epsg4326.tif (and was downloaded from https://figshare.com/articles/dataset/European_LCZ_map/). I also tried with your Shanghai files in the testing folder and it seems that not all the LU_Index urban cells are substituted.

Is there something I am missing?

Thanks a lot for your help!

[dargueso]

Dear @pippo2013, The w2w package gets info from LCZ maps (tiff file) and creates a new geo_em file with that information, replacing the original one from WRF (generally MODIS). Those two datasets do not coincide in all grid points for various reasons. One of them is the interpolation WPS does from MODIS to your WRF grid. Therefore, the w2w will redefine some rural grid points to urban, but it could also happen the opposite.

In your case, the urban grid points over the ocean may be due to the fact that there are small land areas there (islands, land strips, etc). You can confirm this with any mapping/satellite product such as Google Maps.

WRF will likely ignore those grid points because they will be masked out by the landmask variable, but I would need to confirm this. @matthiasdemuzere shall we use the landmask to mask them out before writing the files out? How would you like to proceed with this?

[matthiasdemuzere]

@dargueso: for the Italy case, I have the impression that there is truly a shift in the while grid? Looking at the shape of the pixels that are shifted, it is as if the the LCZ urban information is shifted to the South (South of Napoli), yet also to the East (North of Pescara) ...

I am not sure what is going on here .... but I don't think using the landmask is the right approach here, as this would still lead to a bad urban representation (assuming that the LCZ layer is shifted as a whole).

I should probably check the source LCZ EU data first, to make sure this re-projected version is ok?

[dargueso]

@matthiasdemuzere, yes, you're right. I was looking only at the south, but looking to the east of the domain, it is quite clear that it is shifted... Overall, the main city (Napoli?) is shifted south.

[pippo2013]

Dear @dargueso, thanks a lot for your reply. I understand the possible mismach between original LCS map due to the upscaling to WPS grid. However, at least in my data, there is a shift between the two geo_em files (the original WPS file and the one obtained with w2w, namely geo_em.d03_LCZ_extent.nc and geo_em.d03_LCZ_parameter.nc). Rome is clearly shifted towards South-East and the red spots in the sea draw the Eastern Italian coastal line.

[pippo2013]

Thank you @dargueso and @matthiasdemuzere , I wrote my last message before refreshing the page! I double checked LCZ EU dataset and is correct. But I cannot find the way out... Thank you for checking.

[matthiasdemuzere]

@pippo2013

The _epsg4326 version was created afterwards, to easen some of the processing.

But this new version of w2w is able to detect different projection system, and reproject if needed. Could you therefore re-try your analysis using the original EU LCZ tif?

I just want to see whether perhaps the reprojected EPSG:4326 version is wrong ...

I actually tried to quickly process this myself, see output here.

Perhaps you can first have a look at this, to see if it is any better?

[pippo2013]

Thanks @matthiasdemuzere for your prompt help! I plotted your w2w data but the result is just the same

I cannot do it immediately, but I will try to use the original EU LCZ tif asap and let you know.

[pippo2013]

Dear @matthiasdemuzere, I have tried with the orginal EU LCZ tif, as you suggested. But unfortunately I still experience the artifact. Hoping to identify a possible mistake in the procedure I am using, I also did a test on US and implemented a grid on Miami, but this time everything worked like a charm. I attach in the following the two plots with the overlaid coasts.

I also tried different parts of Europe but obtained the same artifact I got for Italy. Thanks in advance for your help!

[matthiasdemuzere]

Hi @pippo2013,

I am actually a bit puzzled by the behaviour ... not sure where to look for a solution.

As a start, could you report here your w2w version: > w2w --version

and the versions of the packages installed in your environment? > pip list

[matthiasdemuzere]

@pippo2013:

I double-checked the projection definition of the European LCZ map, and I believe something was wrong in that file. I have fixed this, and uploaded a new version of the file to Zenodo (here).

May I ask you to:

1. Download the EU map again, this time only available in the European LAEA projection
2. Clip your area of interest, that is larger in all directions than your geo_em file
3. Run w2w again -> the routine should also reproject the LCZ data
4. Analyse the results
5. Report back here.

I hope this will solve your issue.

[EDIT] I actually tried this myself, and I still get the same strange behavior ... so bear with me until I find whats is going on here ...

[matthiasdemuzere]

A small update.

So, I double checked the new European LCZ map in QGIS. I do think the projection definition is now ok.

I also checked the internal w2w reprojection, done here. This also looks fine to me.

As such, I think the problem is situated in the geo_em file itself?

In order to resample between the LCZ tif and the geo_em grid definition, w2w writes out a temporary .tif file derived from the geo_em netcdf file. This is done here, after the replacement of the urban LC with surrounding natural LC (for which the output is ok).

So I get the feeling the problem lies in the geo_em file itself?

Here you can find the temporary _gridinfo.tif that is extracted from the geo_em netcdf file. If you open this one in QGIS, you'll see that there is a shift compared to e.g. the openstreetmap layer ...

So perhaps you could first look into your geo_em netcdf file?

[pippo2013]

Dear @matthiasdemuzere, thanks a lot for your help. I made many tests and verified that for smaller domains the artifact disappears (as in the plot I posted for Miami, while for larger domains the artifact is apparent no matter if the domain is in Europe or in US.) I include here a test in Rome for a small domain (without shift) and the bigger one with a visible shift.

I double checked my geo_em file and it seems correct, I plotted results independently with Vapor and with WRF-python and they are fine, without any shift.

However, comparing your geo_em file on Spain in the test folder, I assume (from the large map_factor and from the regular lat/long grid) they were obtained through a Mercator projection, while in my WPS file I used Lambert conformal one. Hence I think that the problem may arise due to the fact that the temporary tif file just uses X_LAT and X_LONG instead of considering all the parameters (contained in the geo_em file) used to project the domain on the secant plane with the Lambert projection. May this be the reason?

Thank you so much for your time and help!

[matthiasdemuzere]

Thanks for these additional tests.

Though it is not solved yet, I do think it is getting more clear where the problem starts! That's a start. I'll look into this later on ...

[pippo2013]

Thanks a lot! :)

[pippo2013]

Hi! I read this very interesting post on WRF projection https://fabienmaussion.info/2018/01/06/wrf-projection/ it explains how to properly convert WRF Lambert conformal projection (using all the geo_em parameters and considering the earth as a sphere) and shows how shift may arise if this is not taken into account. I hope it helps.

[matthiasdemuzere]

Ha ... Fabien to the rescue. This is great, I'll have a look at this.

[matthiasdemuzere]

After testing a lot of things, I now understand where I made a mistake. I am treating the WRF coordinates as if they are on a regular grid, by taking the first row/column for lat and lon values, see here.

Yet since WRF is on an irregular grid, I should not do this.

@dargueso or @theendlessriver13 : I wonder if you have any ideas on how to do this?

So the procedure I use is:

1. I create a temporary .tif grid file from the WRF netcdf file, via create_wrf_gridinfo.
2. In sub-sequent steps, I use this .tif for all reprojections (see eg. here), to reproject the higher resolution LCZ grid to the coarser resolution WRF grid.

But since the tmp .tif grid file definition is wrong, also the resulting reprojection is wrong.

So we need to come up with a better approach for the Reprojection. Yet I am not sure how this should best be done? Any ideas?

[pippo2013]

Thanks a lot @matthiasdemuzere for your investigation on this issue. Unfortunately I do not have much experience on this topic. I am looking forward to your news!

[dargueso]

Hi @matthiasdemuzere, sorry for my late reply, I was travelling. The problem with lat/lon vs Lambert is always tricky. I would have to dig deeper to understand the whole process in order to give an opinion. I guess that the tif files must be regular lat/lon. Do you need to use a tif file for the reprojections in step 2? Because the alternative is to directly read the WRF grid and interpolate on to that grid and skip step 1, if that is possible. Does this sound like a feasible approach?

[matthiasdemuzere]

That is indeed what I have in mind as well. Mainly have to figure out how to get the right Transform tuple for the irregular wrf grid ...

If that is in place, it should work without the temporary .tif

[dargueso]

@matthiasdemuzere, I don't really know what you mean by get the right transform tuple. Do you mean the projection parameters from WRF grid? or else? Because WRF projection parameters will change from domain to domain, even if both use Lambert - that information is embedded in the geo_em file, but some care must be taken when interpreting. I can help you with that or other approach transforming from/to the WRF grid, if that would be helpful.

[matthiasdemuzere]

That's indeed what I mean. Extract the proper projection parameters from the submitted geo_em wrf file. And then use these into the reproject directly ...

Any support with this would be appreciated!

[dargueso]

@matthiasdemuzere, the projection parameters in WRF files are a nightmare. I haven't managed to fully understand how they are interpreted by WRF to generate the domain and they don't always follow a standard, as far as I can tell. Plus there are a few projections that are commonly used. For example, Lambert, Mercator, Lat/lon and rotated pole. One option is to use wrf-python, which is able to read WRF files and extract the projection information into cartopy. I have often used it to make plots in WRF projection or transform the plots into a different projection, but not sure how to use it in other packages.

[NicMan89]

@matthiasdemuzere @dargueso thanks a lot for sharing w2w tool! I have tried to extract wrf grid projection following the link posted from @pippo2013. The projection string is something like this: 'proj=lcc lat_0=41.7900199890137 lon_0=7.25 lat_1=30 lat_2=50 x_0=0 y_0=0 R=6370000 units=m no_defs' generate from: wrf_proj = pyproj.Proj(proj='lcc', lat_1=wrf.TRUELAT1, lat_2=wrf.TRUELAT2, lat_0=wrf.MOAD_CEN_LAT, lon_0=wrf.STAND_LON, a=6370000, b=6370000)

But I'm not sure this string is correct to generate the same grid in the geo_em file. After that I have tried to generate a .tif with LU_INDEX variable inside (following the same function name in the w2w.py script). I used for this rasterio package inside the attached file create_wrf_gridinfo.txt. Now the problem is that: the output .tif file is flipped horizontally and vertically. As @matthiasdemuzere says i think the mistake is in the transform tuple. I hope this may help and I apologize in advance for any inaccuracies in the code.

[pippo2013]

Hi everyone! I tried to follow the approach you suggested, using salem reprojection of WRF grid (from geo_em file) onto EU_LCZ map:

ds2 = salem.open_xr_dataset('EU_LCZ_map_epsg4326.tif') dse = salem.open_xr_dataset('geo_em.d03.nc') lu_reproj = ds2.salem.transform(dse.LU_INDEX)

The projection is fine. However, I think it would be better to skip the geotiff since WRF has an irregular grid.

[NicMan89]

Hi everyone. Following the @pippo2013 message, I have tried to use salem package to reproject native LCZ.tif in wrf grid from geo_em.

import salem
ds2 = salem.open_xr_dataset('EU_LCZ_map_epsg4326.tif')
dse = salem.open_xr_dataset('geo_em.d03.nc')

invlu_reproj = dse.salem.transform(ds2)

The result is:

In [50]: invlu_reproj
Out[50]: 
<xarray.Dataset>
Dimensions:      (west_east: 390, south_north: 399)
Coordinates:
  * west_east    (west_east) float64 2.36e+05 2.37e+05 ... 6.24e+05 6.25e+05
  * south_north  (south_north) float64 -1.66e+05 -1.65e+05 ... 2.31e+05 2.32e+05
Data variables:
    data         (south_north, west_east) int16 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
Attributes:
    pyproj_srs:  +proj=lcc +lat_0=41.7900199890137 +lon_0=7.25 +lat_1=30 +lat...

with the shape and projection of wrf grid and land use from lcz.tif as in the plot

I don't know how the salem package resample the higher resolution LCZ grid to the coarser resolution WRF grid. I hope this may help.

EDIT https://github.com/fmaussion/salem/blob/0ae5885e4f466fc333a14544df851501f11102fd/salem/gis.py#L784-L873 I think this is the function in salem package used to resample

[matthiasdemuzere]

Hi @NicMan89 @pippo2013

Thanks for looking into this.

I also reached out to Fabien, developer of Salem. He said salem could indeed be used for this, or alternatively the xesmf package.

Unfortunately I don't manage to install xesmf via pip. So we could use salem.

Yet the transformis typically used to interpolate, with the output grid of similar of finer resolution than the input grid (see here), with as default using the nearest neighbour value. So this is not what we need.

Aggregation is also possible, using lookup_transform. Yet am not sure yet how we can work with the mode, which we need to add the modal LCZ classes into the LU_INDEX table.

I'll keep looking what would be the best way to do this, ideally with the least additional packages.

[matthiasdemuzere]

@pippo2013

I believe I managed to come up with a solution. I tested this on your data, and the results look consistent now. From left to right, LU_INDEX from geo_em.d03.nc, geo_em.d03_LCZ_extent.nc, and geo_em.d03_LCZ_params.nc

It is still work in progress, as some of our tests still fail with this new implementation. But it seems we are getting there.

If you want, you can try the branch out by installing: pip install git+https://github.com/matthiasdemuzere/w2w@fix_wrf_reproject

[pippo2013]

Thank you so much @matthiasdemuzere!! I will test it immediately and let you know! :)

[pippo2013]

Dear @matthiasdemuzere, sorry for my late reply. I struggled more than expected when installing the new w2w branch (compatibilities with libraries). In the end, after beginning from scratch, I managed to install it and everything seems to work fine! I made a test also for another domain, apart from the Italy case. Everything works properly, no shift for the new LCZ domain, nor for the parent domains (the ones having label _41). Thanks a lot!! I read in the pull section that the problem still persists for stereographic projection, but my specific problem for the Lambert projection is ok, hence I guess this issue can be closed.

[matthiasdemuzere]

Hi @pippo2013 :

Thanks for testing this.

Questions / comments:

• could you provide some more details on the library issues you had? Something we should address?
• great that there is no longer a shift in the LCZ implementation
• the _41 will not show a shift anyhow, as we only change some global attributes in these files (eg. NUM_CAT_LAND).

[pippo2013]

Hi @matthiasdemuzere, my answers:

• I got problems with gdal and numpy (to upgrade), but nothing you should address. I suppose my environment got messed up
• yes, great! I look forward to using WRF with the correct LCZ!
• you are right for _41, it was a dummy test!

Thank you again for your help!!