Enable Copernicus DEM #374
Conversation
…s the numerical accuracy for e.g. computing the ASTER grid spacing, it produces erroneous results when the boundary coordinates of a tile does not exactly match integer coordinates (e.g. the case for MERIT). Furthermore, the signs for adding ‘half_gridp’ to latitudinal coordinates is erroneous because the decreasing nature of latitude values was already accounted for while reading them in.
|
launch jenkins |
|
Your PR changes many result by a lot: [------------] TEST clm/ecoclimap_sg: Ecoclimap-SG landuse, Bolivia
Starting test
Executing: ./run_extpar_icon.sh
Test finished
[ MATCH ] existence_extpar_out.sh
PARAM ABSDIFF DIFF/TOTAL
HORIZON 0.41544 20878/42496
HORIZON 0.38630 20414/42496
HORIZON 0.23093 20748/42496
HORIZON 0.24857 21610/42496
HORIZON 0.31576 22193/42496
HORIZON 0.20184 23824/42496
HORIZON 0.19192 23283/42496
HORIZON 0.28927 23742/42496
HORIZON 0.18990 24689/42496
HORIZON 0.26653 24841/42496
HORIZON 0.31748 23344/42496
HORIZON 0.34628 23445/42496
SKYVIEW 0.00012463 14034/42496
SSO_OROMAX 486.35 11120/42496
SSO_OROMIN 132.07 12049/42496
SSO_STDH 13.052 38087/42496
T_CL 0.18784 41564/42496
Z0 0.23920 37700/42496
topography_c 28.898 38076/42496
Differences from reference found
[ FAIL ] tolerance_check_cdo.py
[ FAIL ] RESULT clm/ecoclimap_sg: Ecoclimap-SG landuse, Bolivia
Creating directory for icon_ecci
Fetching executable /workspace/test/testsuite/run_extpar_icon.sh
Modify namelists (according to XML specification)
[------------] TEST dwd/icon_ecci: ECCI landuse, Bolivia
Starting test
Executing: ./run_extpar_icon.sh
Test finished
[ MATCH ] existence_extpar_out.sh
PARAM ABSDIFF DIFF/TOTAL
SSO_OROMAX 486.35 11120/42496
SSO_OROMIN 132.07 12049/42496
SSO_STDH 13.052 38087/42496
Z0 0.23920 37669/42496
topography_c 28.898 38076/42496
Differences from reference found
[ FAIL ] tolerance_check_cdo.py
[ FAIL ] RESULT dwd/icon_ecci: ECCI landuse, Bolivia
Creating directory for corine_icon
Fetching executable /workspace/test/testsuite/run_extpar_icon.sh
Modify namelists (according to XML specification) |
| @@ -294,10 +294,10 @@ SUBROUTINE fill_topo_data(raw_data_orography_path,topo_files, & | |||
| ! reads in the last latitude value of tile i | |||
| CALL check_netcdf(nf90_close(ncid)) | |||
| ! the netcdf file is closed again | |||
| tiles_lon_min(i) = REAL(NINT(tiles_lon_min(i) - half_gridp)) !< half of a grid point must be | |||
| tiles_lon_max(i) = REAL(NINT(tiles_lon_max(i) + half_gridp)) !< added, as the ASTER/GLOBE/MERIT/COPERNICUS data | |||
| tiles_lat_min(i) = REAL(NINT(tiles_lat_min(i) + half_gridp)) !< is located at the pixel center | |||
There was a problem hiding this comment.
I think the NINT could be necessary for ASTER and GLOBE to ensure that we get a natural number for lat_min and lat_max, since the bounds of the tiles are always at natural numbers for these two dataset (by not using NINT we could have cases where the float value is almost equal to the correct natural number, but not exactly). However, this is not the case for MERIT, where the cells are shifted by half a grid point and thus their edges are not at integer latitude values. Therefore I think we should remove NINT only for the MERIT (and Copernicus) dataset, which is what you propose in "Option B" (see email conversation).
Regarding the signs, I think these are correct as they are and don't need to be swapped for ASTER, GLOBE and MERIT, since the latitude values in the NetCDF files are laid out "from North to South". Note that tiles_lat_min and tiles_lat_max are not the actual minimum and maximum latitude values in the tiles, but simply the first and last element of the lat field, respectively. Not sure how the lat field looks in the Copernicus dataset though.
There was a problem hiding this comment.
Correct, I think the NINT makes sense for ASTER and GLOBE because it yields integer values that are more precise. I removed this suggestion from the current pull request and will add a new branch / pull request to handle MERIT/Copernicus DEM correctly. I'm still confused about the signs - I will look into that once more in the same new branch.
There was a problem hiding this comment.
Ok, perfect!
Actually, you are right, the signs also need to be swapped. I was fooled by the fact that these are also swapped compared to the reading of longitude min and max:
CALL check_netcdf(nf90_get_var(ncid, varID_lat, tiles_lat_max(i), start = (/1/)))
! reads in the first latitude value of tile i
CALL check_netcdf(nf90_get_var(ncid, varID_lat, tiles_lat_min(i), start = (/tiles_nrows(i)/)))
! reads in the last latitude value of tile i
So tiles_lat_max is the first element of the lat field and tiles_lat_min is the last one. In the previous message I wrote the opposite.
Sadly, I think this means more work, because this is not the only place in the code where these signs are swapped. E.g., see lines 327-328 in mo_topo_routines.f90:
topo_grid%start_lat_reg = aster_lat_max + 0.5_wp * dlat ! latitude from north to south, note the negative increment!
topo_grid%end_lat_reg = aster_lat_min - 0.5_wp * dlat ! latitude from north to south, note the negative increment!
Judging from the corresponding code for lon these signs need to be swapped as well. BTW, this computations seem unnecessary, since they are basically undoing the one we do for tiles_lat_min and tiles_lat_max (I think you mentioned that when we met in person as well). Maybe I should look a bit more into this and clean up the code to remove/simplify such things.
There was a problem hiding this comment.
I think in the above lines, the computation is actually correct. Because here, dlat itself is negative, and you want to compute the centre coordinates of the first and last grid cell. So you subtract half a cell from the upper bound (max) and add half a cell from the lower bound (min).
Generally, I think one has to be extremely careful about making chances here because certain terms are not used consistently. And I agree - occasionally there is a lot of forth and back computing, which seems superfluous and makes following the code extra hard...
There was a problem hiding this comment.
You are right, dlat is negative, then this is fine. Yes, I think I should look at this code and make it a bit more intuitive and less prone to errors.
For now let's proceed as discussed, then I will look into the topography code in a separate PR to not slow you down any further.
…increases the numerical accuracy for e.g. computing the ASTER grid spacing, it produces erroneous results when the boundary coordinates of a tile does not exactly match integer coordinates (e.g. the case for MERIT). Furthermore, the signs for adding ‘half_gridp’ to latitudinal coordinates is erroneous because the decreasing nature of latitude values was already accounted for while reading them in." This reverts commit 238a2e7.
Enables the usage of Copernicus DEM (30 m resolution) as input topography data for EXTPAR.
https://dataspace.copernicus.eu/explore-data/data-collections/copernicus-contributing-missions/collections-description/COP-DEM
Copernicus DEM is a high-resolution (30 m) topography data set, which is based on recent data. Its superior accuracy compared to other DEMs was demonstrated in various publications (e.g., Guth (2021)). It is a particular useful as a replacement for ASTER data for high-resolution simulations (~1km grid spacing and smaller) and for regions with complex terrain (ASTER data reveals, for instance, severe artefacts in the Alps - particularly at steep north-facing slopes).
References:
Guth, P. L. and Geoffroy, T. M. (2021). LiDAR point cloud and ICESat-2 evaluation of 1 second global digital elevation models: Copernicus wins. Transactions in GIS, 25, 2245–2261.
Python script to download and prepare Copernicus DEM data:
process_copernicus_dem.txt