SRTM data into Geogrid -

L Alvarez · Aug 25, 2019

Hello there,

There is information in the old forum, a bit ambiguous, and I can't get results by applying that procedure to correctly process SRTM data in geogrid. I would like that together we can provide a standard methodology to do this.
I tell you my procedure and at the point where I am, where the result has failed.

1) I have downloaded the data from ->
http://dwtkns.com/srtm/
http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_33_07.zip

2) I have compiled the "convert_geotiff" package from ->
https://github.com/openwfm/convert_geotiff

3) And these are the instructions that I have used ->
./convert_geotiff -w 4 -t 6000 -u “meters MSL” -d “Topography height” sourceSRTM.tif
3.1) When I execute the previous command, it tells me that "Too many positional arguments", then i do ->
./convert_geotiff -w 4 -t 6000 sourceSRTM.tif

It seems to end correctly and generates an index file (attached) and another files 0001xxxx. I modify by hand the fields of "units" and "description" of index file.

4) I continue to modify ->
/geogrid/GEOGRID.TBL.ARW
interp_option= 3s: average_gcell(4.0)+four_pt+average_4pt
rel_path= 3s: topo_3s
namelist.wps -> geog_data_res = “3s”

5) And I run geogrid.exe, it ends without failure, but when I open the domain (to verify)->
(min and max both 32768 for variable HGT_M)

Have I done something wrong, any suggestions?..

index file ->

Code:

projection = regular_ll
known_x = 1
known_y = 6000
known_lat = 30.000000
known_lon = -20.000000
dx = 0.00833333
dy = 0.00833333
type = continuous
signed = yes
units = "meters MSL"
description = "Topography height"
wordsize = 4
tile_x = 6000
tile_y = 6000
tile_z = 1
tile_bdr = 3
missing_value = 0.000000
scale_factor = 1.000000
row_order = bottom_top
endian = little

L Alvarez · Aug 25, 2019

Hello there,

This is another procedure that I have performed ->

1) Go to http://dwtkns.com/srtm/
Click on a tile to download it. Download all available GeoTiff tiles for the area of your nest domain
http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_33_07.zip

Each tile contains 4 files: a .txt, .hdr, .tfw and .tif file. This latest SRTM version 4 has a resolution of 90m at the equator and is provided in 5 deg x 5 deg tiles. Their size is 6001x6001 pixels: in fact 5999x5999 pixels with an overlap of 1 pixel on each side (not 6000x6000 pixels as indicated in the text file).

2) The WRF preprocessing system (WPS) requires the static geographic data to be in binary format. This file is read by the geogrid program used to set up the WRF domain. Each SRTM tile must be converted to the WPS binary format via the command:

Code:

$ gdal_translate -of ENVI tilename.tif tilename.bil

The .hdr file is overwritten by a slightly different one and 2 new .bil and .bil.aux.xml files are created.

3) To create the needed index file, use a text editor to copy and paste the following in a new file named "index".

Code:

type = continuous
signed = yes
projection = regular_ll
dx = 0.000833333333333
dy = 0.000833333333333
known_x = 1.0
known_y = 5999
known_lat = 20
known_lon = -30
wordsize = 2
endian = little
tile_x = 5999
tile_y = 5999
tile_z = 1
tile_bdr=1
row_order = top_bottom
missing_value = 32768
units = "meters MSL"
description = "SRTM 3-sec Topography height"

- Next change the known_lat and known_lon. There are two options here:

a. If you were able to download your most southwest (lower left) tile, open the .hdr file for this tile. Locate the known_lat and known_lon values and copy into your index file. For example you might find "map info = {Geographic Lat/Lon, 1, 1, -30, 20, 0.000833333333333333, 0.000833333333333333,WGS-84}" The first value is always known_lon and the second is known_lat.

b. If you were not able to download your most southwest or lower left tile, you obviously can't open its .hdr file to copy the values of known_lat and known_lon into your index file!

c. Rename all .bil files into "00001-ncols.00001-nrows" without an extension. Ncols and nrows are the number of columns and rows in the datafile. Start with the tile in the most southwest or lower left corner. The name of this tile will be "00001-05999.00001-05999".

4) I continue to modify ->

/geogrid/GEOGRID.TBL.ARW
interp_option= 3s: average_gcell(4.0)+four_pt+average_4pt
rel_path= 3s: topo_3s
namelist.wps -> geog_data_res = “3s”

5) And I run geogrid.exe, it ends without failure, but when I open the domain with ncview (to verify)->
(min and max both 32768 for variable HGT_M)

Have I done something wrong, any suggestions?..

L Alvarez

mgduda · Aug 26, 2019

Could you try using the attached script to plot the tile of SRTM data that you've processed with either or both of the methods you've outlined? You'll need to make some changes to the parameters in the script, but these parameters should be explained in the comments. If the resulting plot shows what looks like correct topography, we can look further into the geogrid steps to see what might be going wrong there; otherwise, the problem may be in the processing of the GeoTIFF data.

L Alvarez · Aug 26, 2019

Hello Mduda,

It seems that this time it worked for me. I repeated the first method again.
The image from below to the left, in theory at 90 m (3s). resolution and to the right to 1 km resolution (30s).

In any case, how could I verify that it is really 3s and has not made a resample to another size?..
I'm not familiar with python, I'm embarrassed, I don't know how to run your script.

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

1) I have downloaded the data from ->
http://dwtkns.com/srtm/
http://srtm.csi.cgiar.org/wp-content/up ... _33_07.zip

2) I have compiled the "convert_geotiff" package from ->
https://github.com/openwfm/convert_geotiff

3) When I execute the previous command, it tells me that "Too many positional arguments", then i do ->
./convert_geotiff -w 4 -t 1200 sourceSRTM.tif

I think, although my tile has 6000 points, I have put 1200 points and I think what it has done is to break it down.

4) The application automatically created an index file ->

Code:

projection = regular_ll
known_x = 1
known_y = 6000
known_lat = 30.000000
known_lon = -20.000000
dx = 8.333334e-04
dy = 8.333334e-04
type = continuous
signed = yes
units = "meters MSL"
description = "Topography height"
wordsize = 4
tile_x = 1200
tile_y = 1200
tile_z = 1
tile_bdr = 3
missing_value = 0.000000
scale_factor = 1.000000
row_order = bottom_top
endian = little

I had to make this little change ->

Code:

missing_value = 32768

And my GEOGRID.TBL file ->

Code:

===============================
name = HGT_M
        priority = 1
        dest_type = continuous
        smooth_option = smth-desmth_special; smooth_passes=1
        fill_missing=0.
        interp_option = gmted2010_30s:average_gcell(4.0)+four_pt+average_4pt
        interp_option =     gtopo_30s:average_gcell(4.0)+four_pt+average_4pt
        interp_option =      gtopo_2m:four_pt
        interp_option =      gtopo_5m:four_pt
        interp_option =     gtopo_10m:four_pt
        interp_option =        lowres:average_gcell(4.0)+four_pt
        interp_option =       default:average_gcell(4.0)+four_pt+average_4pt
        interp_option = 3s:average_gcell(4.0)+four_pt+average_4pt
        rel_path = gmted2010_30s:topo_gmted2010_30s/
        rel_path =     gtopo_30s:topo_30s/
        rel_path =      gtopo_2m:topo_2m/
        rel_path =      gtopo_5m:topo_5m/
        rel_path =     gtopo_10m:topo_10m/
        rel_path =        lowres:topo_gmted2010_5m/
        rel_path =       default:topo_gmted2010_30s/
        rel_path = 3s:topo_3s/

Rules!!

L Alvarez · Aug 27, 2019

Hello Mduda,

According to the procedure above, how could I verify that it is really 3s and has not made a resample to another size?..

On the other hand, I found an MDT at 25 meters resolution that I tried to implement with the same method above.

the index file is ->

Code:

projection = mercator
truelat1 = 0.000000
stdlon = 0.000000
known_x = 1
known_y = 2747
known_lat = 28.605396
known_lon = -16.940907
dx = 2.594085e-04
dy = 2.220986e-04
type = continuous
signed = yes
units = "NO UNITS"
description = "NO DESCRIPTION"
wordsize = 4
tile_x = 1200
tile_y = 1200
tile_z = 1
tile_bdr = 3
missing_value = 0.000000
scale_factor = 1.000000
row_order = bottom_top
endian = little[code]

I have modified it to ->

Code:

projection = regular_ll
known_x = 1
known_y = 2747
known_lat = 28.605396
known_lon = -16.940907
dx = 2.594085e-04
dy = 2.220986e-04
type = continuous
signed = yes
units = "meters MSL"
description = "Topography height"
wordsize = 4
tile_x = 1200
tile_y = 1200
tile_z = 1
tile_bdr = 3
missing_value = 32768
scale_factor = 1.000000
row_order = bottom_top
endian = little

This is my GEOGRID.TBL.ARW file ->

Code:

name = HGT_M
        priority = 1
        dest_type = continuous
        smooth_option = smth-desmth_special; smooth_passes=1
        fill_missing=0.
        interp_option = gmted2010_30s:average_gcell(4.0)+four_pt+average_4pt
        interp_option =     gtopo_30s:average_gcell(4.0)+four_pt+average_4pt
        interp_option =      gtopo_2m:four_pt
        interp_option =      gtopo_5m:four_pt
        interp_option =     gtopo_10m:four_pt
        interp_option =        lowres:average_gcell(4.0)+four_pt
        interp_option =       default:average_gcell(4.0)+four_pt+average_4pt
        interp_option = 3s:average_gcell(4.0)+four_pt+average_4pt
        interp_option = 1s:average_gcell(4.0)+four_pt+average_4pt
        rel_path = gmted2010_30s:topo_gmted2010_30s/
        rel_path =     gtopo_30s:topo_30s/
        rel_path =      gtopo_2m:topo_2m/
        rel_path =      gtopo_5m:topo_5m/
        rel_path =     gtopo_10m:topo_10m/
        rel_path =        lowres:topo_gmted2010_5m/
        rel_path =       default:topo_gmted2010_30s/
        rel_path = 3s:topo_3s/
        rel_path = 1s:topo_1s/

Running this in geogrid.exe, it ends without apparent failure, but when I'm going to plot it with ncview (image below)->
does not display the contour of the relief/terrain, the values on the "scale are 2e+08 to 1.4+09", but navigating above the image show correct height.
It may be that the file is fine, but the scale is not and therefore does not paint the terrain meters, I do not know...

How could I fix this?

Best regards,

L Alvarez

mgduda · Aug 28, 2019

When you ask "how could I verify that it is really 3s and has not made a resample to another size", are you asking whether the "convert_geotiff" program has resampled the data or whether geogrid has resampled the data? The geogrid program will interpolate from the native resolution of the data to your WRF model grid resolution; in your case, you've used the "average_gcell" method in your GEOGRID.TBL file, which instructs the geogrid program to compute the average value of all pixels of source data that lie within each WRF grid cell (i.e., a grid-cell average). I've never used the "convert_geotiff" program so I couldn't say with certainty what it is doing, but I would assume that it does not perform any resampling of the data.

mgduda · Aug 28, 2019

Regarding the 25-m MDT data, there are a few points that you could verify:

1) Your index file has tile_x = 1200, tile_y = 1200, tile_z = 1, tile_bdr = 3, and wordsize = 4, so you could verify that your 00001-01200.00001-01200 file has a size of exactly (3 + 1200 + 3) * (3 + 1200 + 3) * 4 = 5817744 bytes.

2) The range of values for HGT_M in your geo_em.d04.nc suggests that there is a missing value in the dataset other than the value of 32768 that you've specified in your index file; can you determine the value for the HGT_M field in the south-east part of your domain 4, and try some missing values in your index file that are nearby to that value? For example, a dataset might reasonable use values of 2^32 - 1 as a missing value for an unsigned, 32-bit integer field; or, if the values are stored in the MDT dataset as signed 32-bit integers, a missing value of 2^31 - 1 or -2^31 might be reasonable.

3) In case the data are stored in big-endian byte order in your 00001-01200.00001-01200 file, you could try changing "endian = little" to "endian = big" in your index file to see if that helps.

L Alvarez · Aug 30, 2019

mgduda said:
When you ask "how could I verify that it is really 3s and has not made a resample to another size", are you asking whether the "convert_geotiff" program has resampled the data or whether geogrid has resampled the data? The geogrid program will interpolate from the native resolution of the data to your WRF model grid resolution; in your case, you've used the "average_gcell" method in your GEOGRID.TBL file, which instructs the geogrid program to compute the average value of all pixels of source data that lie within each WRF grid cell (i.e., a grid-cell average). I've never used the "convert_geotiff" program so I couldn't say with certainty what it is doing, but I would assume that it does not perform any resampling of the data.

Hello mgduda,

In other words, if "convert_geotiff" uses the source data at 90 meters (srtm data origen), the geogrid takes them and interpolates the grid resolution of that domain.

For example, if our grid resolution is 200 or 500 meters, the gegorid will interpolate and convert to that 200 and 500 m. resolution topographic too, is that correct?.. both are at the same resolution?..

If this is correct, is there any way to maintain the grid resolution for example at 200 m. and for the topography to be interpolated at 500 m.

L Alvarez

L Alvarez · Sep 2, 2019

mgduda said:
Regarding the 25-m MDT data, there are a few points that you could verify:

1) Your index file has tile_x = 1200, tile_y = 1200, tile_z = 1, tile_bdr = 3, and wordsize = 4, so you could verify that your 00001-01200.00001-01200 file has a size of exactly (3 + 1200 + 3) * (3 + 1200 + 3) * 4 = 5817744 bytes.

2) The range of values for HGT_M in your geo_em.d04.nc suggests that there is a missing value in the dataset other than the value of 32768 that you've specified in your index file; can you determine the value for the HGT_M field in the south-east part of your domain 4, and try some missing values in your index file that are nearby to that value? For example, a dataset might reasonable use values of 2^32 - 1 as a missing value for an unsigned, 32-bit integer field; or, if the values are stored in the MDT dataset as signed 32-bit integers, a missing value of 2^31 - 1 or -2^31 might be reasonable.

3) In case the data are stored in big-endian byte order in your 00001-01200.00001-01200 file, you could try changing "endian = little" to "endian = big" in your index file to see if that helps.

Hello mgduda,

I have not been able to make much progress.
I answer your suggestions;

1) It's OK.

5817744 ago 27 09:30 00001-01200.00001-01200
5817744 ago 27 09:30 00001-01200.01201-02400
5817744 ago 27 09:30 00001-01200.02401-03600
5817744 ago 27 09:30 01201-02400.00001-01200
5817744 ago 27 09:30 01201-02400.01201-02400
5817744 ago 27 09:30 01201-02400.02401-03600
5817744 ago 27 09:30 02401-03600.00001-01200
5817744 ago 27 09:30 02401-03600.01201-02400
5817744 ago 27 09:30 02401-03600.02401-03600

2) I did not understand this comment. In any case, look for the value of the south-east point (1.49451e+07) but the failure continues.

3) I have changed to "endian = big" (below picture)

My topo to 1seg aprox.

https://www.dropbox.com/s/0o5gbozs1j6bmym/136_MDT25_TF.zip?dl=0

My index file ->

Code:

projection = regular_ll
known_x = 1
known_y = 2747
known_lat = 28.605396
known_lon = -16.940907
dx = 2.594085e-04
dy = 2.220986e-04
type = continuous
signed = yes
units = "meters MSL"
description = "Topography height"
wordsize = 4
tile_x = 1200
tile_y = 1200
tile_z = 1
tile_bdr = 3
missing_value = 1.49451e+07
scale_factor = 1.000000
row_order = bottom_top
endian = big

L Alvarez · Sep 26, 2019

Please, some help here!

ac_soaring · Mar 24, 2021

A rather old post, but hopefully this will help someone.

I was able to load SRTM3 data yesterday (3/23/2021) using this post:
https://wrfexplorer.wordpress.com/2015/03/24/high-resolution-topography/ (very similar to your initial method at the top using convert_geotiff).

Except, I had to change the geotiff command to 2 byte words and use a tile six of 6000:

./convert_geotiff -w 2 -t 6000 -u “meters MSL” -d “Topography height” sourceSRTM.tiff

Then, I needed to change the missing data value in the resulting index file from 0 to 32768.

It seemed to work.

Good luck

Ming Chen · Mar 29, 2021

Many thanks for the kind information !

ac_soaring · Sep 10, 2021

The source of SRTM data mentioned in the web post is no longer available. USGS says the data is still available through EarthExplorer:

Thank you for contacting USGS Earth Resources Observation and Science (EROS) center.

Please navigate to: https://www.usgs.gov/centers/eros/science/usgs-eros-archive-digital-elevation-global-multi-resolution-terrain-elevation?qt-science_center_objects=0#qt-science_center_objects

You can download imagery using our EarthExplorer viewer: https://earthexplorer.usgs.gov/

You will need to download GMTED2010 data. Please view this link to read more about GMTED2010 data:

Follow the tabs in EarthExplorer to conduct your search.

1. Enter Search Criteria to find your area of interest. Alternatively, zoom into your area of interest and click the Use Map button to search the area fitted on your screen. Here is a helpful video tutorial: https://www.youtube.com/watch?v=Dz1JN7TJAvA

2. Select your Data Sets.

3. Enter Additional Criteria if necessary.

4. View your Results. From this screen, you can navigate between the original Data Sets you selected. Hover over the icons to see their functions. The following menu items appear.

Show Footprint
Show Browse Overlay <-------------- Use this option to preview the imagery
Compare Browse
Show Metadata and Browse
Download Options <-------------- Use this option to download
Add to Bulk Download
Order Scene
Exclude Scene from Results

Depending on what you are trying to download, your options may be slightly different. Clicking on one of the blue buttons will work to download the images.

Please let me know if I can further assist you.

adriana · Jan 6, 2022

Dear ac_soaring

I followed your solution (./convert_geotiff -w 2 -t 6000 -u “meters MSL” -d “Topography height” sourceSRTM.tiff
Then, I needed to change the missing data value in the resulting index file from 0 to 32768.), but it did not work for me.

When I used :./convert_geotiff -w 2 -t 6000 -u “meters MSL” -d “Topography height” sourceSRTM.tiff ; it told me there are many arguments, and I used :./convert_geotiff -w 2 -t 6000 sourceSRTM.tiff.

May I ask you how can I sort the issue out please?
Best wishes,

ac_soaring · Jan 11, 2022

Adriana:

Well it may be something simple like the odd quotes that this forum uses. The command should be:

./convert_geotiff -w 2 -t 6000 -u "meters MSL" -d "Topography height" sourceSRTM.tiff

If that's not it, I used a version called convert_geotiff-0.1.0 so perhaps there is a newer version?

I still have the source code for the version I downloaded.

I've turned on notifications for this topic so please let me know.

Best,
Alan

adriana · Jan 14, 2022

Dear Alan,
Many thanks for your help. I could follow your solution successfully. But, I face a new issue! After converting the SRTM file to binary raster,and change the following lines :
/geogrid/GEOGRID.TBL.ARW
interp_option= 3s: average_gcell(4.0)+four_pt+average_4pt
rel_path= 3s: topo_3s
namelist.wps -> geog_data_res = “3s”

I got the same result for topography as the default one!( attached)
I wonder why did happen? I expect high resolution and more details for topography using SRTM data.
May I ask you please help me in this regard?
Best wishes,

adriana · Jan 14, 2022

Dear Alan,
I forgot to mention that I downloaded the SRTM data from below link:
http://dwtkns.com/srtm/

ac_soaring · Jan 14, 2022

Adriana:

There seem to be subtle differences between the two images. They are most easily seen in the isolated dots on the lower right. Each of these areas is a slightly different shape and size.

adriana · Jan 15, 2022

Dear Alan,
So many thanks for your help. Actually, I expected more changes for topography.
I also changed the land-use of the model with Esri-30m resolution. The process is as follow:
I converted my land-use data( Esri-30m resolution) classes into MODIS classification. And then used "convert_geotiff function" to produce binary raster files. I also changed the following lines in GEOGRID.TBL.ARW :
landmask_water = 30m:17,21
interp_option = 30m:nearest_neighbor
rel_path = 30m:landuse30m
In WPS name list changed :
geog_data_res = ' 30m'
However after doing all these steps, and running geogrid.exe, I found out that land-use file produced by geogrid.exe has only one value (attached).
I wonder why did this happen?
May I ask you please advice me in this regard?
Best wishes,

ac_soaring · Jan 15, 2022

Adriana:

I'm not much help on this one other than suggesting you check the intermediate steps. Can you confirm that there is more than 1 MODIS classification after the conversion to MODIS classifications? Can you check that the binary raster file has multiple raster values (I think ARCGIS, GRASS or QGIS lets you view and check the pixel values)?

Alan

SRTM data into Geogrid -

New member

Attachments

New member

Attachments

Administrator

Attachments

New member

Attachments

New member

Attachments

Administrator

Administrator

New member

New member

Attachments

New member

New member

Moderator

New member

New member

New member

New member

Attachments

New member

New member

New member

Attachments

New member