Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest

Shifted topography for SRTM1 data

This post was from a previous version of the WRF&MPAS-A Support Forum. New replies have been disabled and if you have follow up questions related to this post, then please start a new thread from the forum home page.

userwrfch

Member
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)
 

Attachments

  • ASTER.jpeg
    ASTER.jpeg
    263.6 KB · Views: 673
  • namelist.wps
    1 KB · Views: 31
  • SRTM1.jpeg
    SRTM1.jpeg
    122.5 KB · Views: 674
  • USGS.jpeg
    USGS.jpeg
    147.9 KB · Views: 672
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:
aster-gmted.png
(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?
 
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!
 
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:
Code:
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:
Code:
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.
 
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!
 
Here's my 'index' file:
Code:
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.
 
Top