Shifted topography for SRTM1 data

[userwrfch]

Dear WRF support team,

I am using WRF version 4.1.3. over several mountain regions with domains of 25 km, 5 km, 1 km, 200 m, 40 m. I downloaded 1) SRTM1 data (~30 m) via the USGS EarthExplorer, and 2) ASTGTM (~30 m) via EarthDataSearch from Nasa, merged all tif files to one big tif file for each SRTM1 and ASTGTM, and then followed steps 4-8 from this post: https://forum.mmm.ucar.edu/phpBB3/viewtopic.php?t=9628

When running geogrid and plotting the output together with shapefiles of the region, I see large shifts for both 1) SRTM1 and 2) ASTGTM. I have attached the plots, where you can see that the shape aligns well when plotted on the original SRTM1/ASTGTM tif file, but not for WRF geo_em output. However, I did a test with the default USGS GMTED2010 30 arc-sec data from WPS (for 30 km, 10 km, 3.33 km, 1.11 km domains), and there is NO shift here.

My input SRTM1/ASTGTM tif files have the coordinate reference system (crs): "+proj=longlat +datum=WGS84 +no_defs".

My WRF-model uses Lambert conformal conic (LCC) projection.

I have read this paper (https://www.researchgate.net/publication/258805934_Overlapping_Interests_The_Impact_of_Geographic_Coordinate_Assumptions_on_Limited-Area_Atmospheric_Model_Simulations), where a huge shift could be seen because the data was transformed from WGS84 to the sphere BEFORE they were ingested and remapped to the WRF domains by geogrid.exe - but this does not seem to be the case for me, if I understand correctly, because I am using WGS84 for my input data.

Can you please help me? Thanks a lot!

My index file for SRTM1 data () looks as follows:

======================================================
projection = regular_ll
known_x = 1
known_y = 21601
known_lat = 55.000137
known_lon = -123.000275
dx = 5.524358e-04
dy = 2.762179e-04
type = continuous
signed = yes
units = "meters MSL"
description = "Topography height"
wordsize = 4
tile_x = 2291
tile_y = 2291
tile_z = 1
tile_bdr = 3
missing_value = 32768
scale_factor = 1.000000
row_order = bottom_top
endian = little


==================================
I have also tried to actually implement different coordinate reference system (crs) transformations to the original SRTM1 file (with "+proj=longlat +datum=WGS84 +no_defs") BEFORE using convert_geotiff:
1) I transformed it to regular lat long, because I had the suspicion that convert_geotiff assumes that --> didn't help, I can still see the shift
2) I transformed it to LCC, but was unsuccessful so far (HGT_M shows 0 in the geo_em files)

 

[mgduda]

I downloaded a couple of tiles of ASTGTM v003 data and manually processed them for use with the WPS through a combination of gdal_translate and a small Python script to remove redundant rows and columns. I'm not seeing any positional shift of the ASTGTM data relative to the GMTED2010 data in my geo_em.d01.nc files:

(The domain using ASTGTM on the left is missing topography in the northern part of the domain because I didn't download enough data to cover the entire domain.)

I've never tried the workflow from the forum post you referenced, but perhaps there's some issue there that may be introducing a shift?

 

[userwrfch]

Thanks a lot for your reply!

Could you please let me know which gdal_translate commands you used (did you transform the reference system to LCC already here, or left it in lat-long)? Did you use convert_geotiff (https://github.com/openwfm/convert_geotiff) at all, or how did you get the data in a format suitable for WPS?

Thank you very much for your time!

 

[mgduda]

Taking the tile 'ASTGTMV003_N60W139_dem.tif' as an example, I first converted to an ENVI raster, which is essentially what the geogrid program uses:

gdal_translate -of ENVI ASTGTMV003_N60W139_dem.tif temp.dat

Then, I removed the last row and column, which overlap (i.e., are redundant with) the first row and column of neighboring tiles, using Python:

import numpy as np
a = np.fromfile('temp.dat',dtype='i2')
b = a.reshape((3601,3601))
a = b[0:3600,0:3600]
a.tofile('00001-03600.00001-03600')

I didn't perform any sort of datum transform or re-projection of the original data.

 

[userwrfch]

Thank you very much for your reply!

However, I still struggle to reproduce your results. Can you please send me the index file that you used for your example?

Thanks a lot!

 

[mgduda]

Here's my 'index' file:

type = continuous
row_order = top_bottom
endian = little
signed = no
projection = regular_ll
dx = 0.000277777777777778
dy = 0.000277777777777778
known_x = 1.0
known_y = 3600.0
known_lat = 61.0001388888889
known_lon = -140.000138888889
wordsize = 2
tile_x = 3600
tile_y = 3600
tile_z = 1
tile_bdr = 0
units = "meters MSL"
description = "ASTGTM 1-arc-second"

I downloaded two tiles of ASTGTM data, ASTGTMV003_N60W140_dem.tif and ASTGTMV003_N60W139_dem.tif, which resulted in two tiles of data in the geogrid dataset named 00001-03600.00001-03600 and 03601-07200.00001-03600, respectively.

 

[userwrfch]

Thank you very much!