Issue with corine land cover

[wrf_alessandro]

Hi All,

I am trying to use the corine land cover in WPS, but after following the steps suggested in the user guide, the geo_em has a constant LU_INDEX and does not display any LANDUSEF. These are the steps I did to process the data:

1) Downloaded CLC from this link:
https://land.copernicus.eu/pan-european/corine-land-cover/clc2018

2) Converted the CLC to a regular projection with a resolution of 30s:
gdalwarp U2018_CLC2018_V2020_20u1.tif U2018_CLC2018_V2020_20u1_WGS84_30s.tif -s_srs EPSG:3035 -t_srs EPSG:4326 -dstnodata "255" -tr 0.008333333 0.008333333

3) I have reclassified the 44 CLC categories to 21 IGBP_MODIS categories
python Reclassify_CLC_to_IGBP_MODIS.py U2018_CLC2018_V2020_20u1_WGS84_30s.tif U2018_CLC2018_V2020_20u1_WGS84_MODIS_30s.tif

4) Converted the Geotiff to binary for WPS:
gdal_translate -of ENVI U2018_CLC2018_V2020_20u1_WGS84_MODIS_30s.tif 00001-15529.00001-5805

5) Then, I have created the index file as follow:

projection = regular_ll
known_x = 1
known_y = 5805
known_lat = 72.663269
known_lon = -56.505142
dx=0.00833333
dy=0.00833333
type = categorical
description = "CORINE LandCover"
wordsize = 1
tile_x = 15529
tile_y = 5805
tile_z = 1
category_min = 1
category_max = 21
row_order = top_bottom
units="category"
description="Noah-modified 21-category IGBP-MODIS landuse"
mminlu="MODIFIED_IGBP_MODIS_NOAH"
iswater=17
islake=21
isice=15
isurban=13

6) Finally, I modified namelist.wps and GEOGRID.TBL as follows:

geog_data_res = 'corine_30s+default',

landmask_water = corine_30s:17,21 # CORINE landcover Calculate a landmask from this field
interp_option = corine_30s:nearest_neighbor
rel_path = corine_30s:corine_30s/

My domain completely falls within the CLC bounds, so I suspect and issue with the descriptor index file; however, I tried to modify it in several ways changing the row_order and missing_value but none of the modifications worked.
I wonder if someone could help me to understand what is going wrong in my procedure; I attach the data and scripts I used in the whole process.

Any suggestion is welcome,
Alessandro

 

[kwerner]

Hi Alessandro,
Are you able to plot the tile to see if everything looks okay? If not, I'm attaching an old Python script I used to plot a tile. You will have to make some adjustments in it, but see if you're able to get it to plot, and let me know. Thanks!

 

[wrf_alessandro]

Hi,

yes, I can plot the binary file I created with "gdal_translate -of ENVI" but then when I run geogrid the LU_INDEX is empty. I tried all the combinations into index, so I do not understand why it is not working. On the other side, I have been able to overcome this issue using convert_geotiff.

Alessandro

 

[kwerner]

Alessandro.
I'm glad to hear you have been able to find a work-around. Hopefully that is enough to allow you to move forward?

 

[wrf_alessandro]

Hi Kwerner,

I have been able to move forward and, after fixing a bug into the code, the model is running fine, but now I am a bit puzzled on my set-up.

I am performing a simulation with one-way nesting, where on the parent domain I set sf_urban_physics=0, while in the nested domains I would activate one of the urban canopy schemes; thus, for consistency, I put num_land_cat=33 in both parent and nested domains, and my LU_INDEX ranges from 1 to 33.
The simulation on the parent domain starts fine and the model recognizes the 33 landcover classes:

INPUT LandUse = "MODIFIED_IGBP_MODIS_NOAH"
 LANDUSE TYPE = "MODIFIED_IGBP_MODIS_NOAH" FOUND          33  CATEGORIES           2  SEASONS WATER CATEGORY =           17  SNOW CATEGORY =           15
Climatological albedo is used instead of table values
INITIALIZE THREE Noah LSM RELATED TABLES
Skipping over LUTYPE = USGS
 LANDUSE TYPE = MODIFIED_IGBP_MODIS_NOAH FOUND          20  CATEGORIES

However, the run stops with a segmentation fault when the model initializes the dry deposition (I am running WRF-Chem). After playing with debugger I found that into module_dep_simple.F the MODIS categories are converted into USGS, but the urban categories 31-33 are not considered:

        if( .not. allocated(luse2usgs) ) then
          allocate( luse2usgs(config_flags%num_land_cat),stat=astat )
          if( astat /= 0 ) then
            CALL wrf_error_fatal( 'dep_init: failed to allocate luse2usgs array' )
          end if
          if( trim(mminlu_loc) == 'USGS' ) then
            luse2usgs(:) = (/ (iland,iland=1,config_flags%num_land_cat) /)
          elseif( trim(mminlu_loc) == 'MODIFIED_IGBP_MODIS_NOAH' ) then
            luse2usgs(:) = (/ 14,13,12,11,15,8,9,10,10,7, &
                              17,4,1,5,24,19,16,21,22,23 /)
          endif
        endif

Thus, I modified the luse2usgs as follows to fix the issue (I also noted that something similar occurs into module_ftuv_driver.F, when using and inizializing F-TUV photolysis):

        if( .not. allocated(luse2usgs) ) then
          allocate( luse2usgs(config_flags%num_land_cat),stat=astat )
          if( astat /= 0 ) then
            CALL wrf_error_fatal( 'dep_init: failed to allocate luse2usgs array' )
          end if
          if( trim(mminlu_loc) == 'USGS' ) then
            luse2usgs(:) = (/ (iland,iland=1,config_flags%num_land_cat) /)
          elseif( trim(mminlu_loc) == 'MODIFIED_IGBP_MODIS_NOAH' ) then
             if( config_flags%num_land_cat .eq. 33) then
               luse2usgs(:) = (/ 14,13,12,11,15,8,9,10,10,7, &
                              17,4,1,5,24,19,16,21,22,23, 16,-999,-999,-999,-999,-999,-999,-999,-999,-999,1,1,1/)
             else
               luse2usgs(:) = (/ 14,13,12,11,15,8,9,10,10,7, &
                              17,4,1,5,24,19,16,21,22,23 /)
             end if
          endif
        endif

Although the simulation now is running fine and I do not see any relevant differences with respect to the simulation using the default MODIS landcover, I wonder if this is correct and I can provide to the model urban categories 31-33 even if I do not use them in the coarse parent domain? If so, would it be better to change module_initialize_real.F to automatically modify grid%ivgtyp=31-31 to 13 when sf_urban_physics=0?

Thanks again for your help,
Alessandro

 

[kwerner]

Hi Alessandro,
As this issue seems to be related to the WRF Chemistry code, can you post the question to the WRF-Chem section of this forum, instead? There is a different group that works with that code and they should be able to help you better. Thanks!