(RESOLVED) New high-resolution topography in WPS / Geogrid

[Arty]

RESOLUTION NOTE :
The issue revolved around both the binary format output from GDAL and the data order. Notably, the origin of the conversion didn't align with the lower-left corner of the raster. GDAL lacks a built-in option for determining the desired endianness.
Furthermore, it's vital to explicitly specify the Int16 binary format, which isn't mentioned in the commonly suggested commands like:

gdal_translate -of ENVI in.tif out.bin

Which can be done this way :

gdal_translate -ot Int16 -of ENVI in.tif out.bin

These commands fall short of addressing the problem, which is why I've provided tools at the end of this thread to overcome the encountered issues.

---

Hello,

I tried to implement a new high resolution topography within WPS/geogrid following the tutorial here, but I can't find the reason why it didn't work. I followed the protocol from this post using the "convert_geotiff" tool, and I obtained the desired binary and index files (see attached). I edited GEOGRID.TBL and WPS namelist.input accordingly.
Nevertheless, the whole HGT_M field in selected domains are filled with zeroes. The geogrid.log indicates at some points (see attached) :

Processing HGT_M
WARNING: In GEOGRID.TBL, no specification for truelat1 in entry 1.
WARNING: In GEOGRID.TBL, no specification for truelat2 in entry 1.
WARNING: In GEOGRID.TBL, no specification for stdlon in entry 1.

I also attached the .tiff file I used.

Also, I checked for the binary files with a small Python script plotting elevation data. It gives fine results ; so the binary file seems OK (see Python script and screenshot attached).

Then I also tried the method exposed in this post. The GDAL tool gives the same binary file as the previous method (compared with the .py script I wrote).

Being new to some programming languages, I didn't yet succeed to use the provided WPS/geogrid tool write_geogrid.c. It seems I have to compile the script first.

Anyway, any help is welcome concerning the warning message from geogrid.log as the problem doesn't seem to come from the binary file.

Thank you

 

[Arty]

Maybe those files can help :

- GEOGRID.TBL
- namelist.wps

 

[Arty]

I figured out I had misnamed the binary file confunding " . " and " - " separators. I will try again after having it renamed properly.

Nevertheless, I'm still wandering about this Mercator projection I have in my index file and if ti correct given the original GeoTIFF file doesn't seem to be projected as such :

gdalinfo tahi.tiff
Driver: GTiff/GeoTIFF
Files: tahi.tiff
Size is 625, 601
Origin = (-150.261250000000217,-16.798750000000002)
Pixel Size = (0.002500000000000,-0.002500000000000)
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-150.2612500, -16.7987500)
Lower Left  (-150.2612500, -18.3012500)
Upper Right (-148.6987500, -16.7987500)
Lower Right (-148.6987500, -18.3012500)
Center      (-149.4800000, -17.5500000)
Band 1 Block=625x3 Type=Float32, ColorInterp=Gray
  Min=0.000 Max=1980.548
  Minimum=0.000, Maximum=1980.548, Mean=17.268, StdDev=106.653
  NoData Value=nan
  Metadata:
    STATISTICS_MAXIMUM=1980.5482177734
    STATISTICS_MEAN=17.268496548744
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=106.65328874517
    STATISTICS_VALID_PERCENT=100

It might just be regular_ll so I should modify my index accordingly. Any help or confirmation on these matters are welcome.

Edit : it seems there are no use of projection when it comes to elevation data, usually gridded under regular_ll format ; which seems quite obvious now I've been confirmed by a friend.

 

[Arty]

Hello all,

I corrected several dumb mistakes mentioned above, but am still stuck - even though I feel really closer to the goal (see screenshot).
I could really use your help, so please leave a comment if you have any idea how I can overcome the issue.
All needed files are attached.

I have tested several index combinations (signed/unsigned | wordsize=2,4,8 | scale_factor=0.000001 | endianness ...) but nothing seems to do the trick ; and the current index attached below is the closer I can get from the desired result despite elevation value are far off (and kind of "flat").

Thank you so much for your help.

 

[Arty]

After dedicating several days to resolving this issue, I resorted to "reverse engineering" the WRF/WPS topography files (topo_30s). Using GDAL and Python, I've developed a compact script that can convert a GeoTIFF topography file into a binary format that's compatible with WRF/WPS. This script also generates the required index file. Please don't hesitate to reach out if you encounter any issues while using it. While the script works in my specific case, it may require some adjustments to ensure broader compatibility.
Note: the script cannot manage the creation of multiple tiles or the use of borders; it is therefore limited according to the WRF/WPS format to a single tile of maximum size 99999x99999.

Step 1: Open GeoTIFF input file
Step 2: Retrieve information from the GeoTIFF file
- Extract the lower-left corner coordinates directly from the dataset
- Round the values to a maximum of 5 decimals
- Create the index file with the specified variables
Step 3: Convert GeoTIFF elevation data to binary Int16
Step 4: Convert big-endian binary Int16 data to little-endian and flip x-axis
- Save the flipped data to the binary file
Step 5: Export the final binary file with the desired name format

Additionally, I'm providing a small script that enables you to visualize the topography based on the binary format used in WRF/WPS. You'll need to modify the information derived from the index file at the beginning of the script to make it work for your specific needs.

---

RESOLUTION NOTE :
The issue revolved around both the binary format output from GDAL and the data order. Notably, the origin of the conversion didn't align with the geographic lower-left corner. GDAL lacks a built-in option for determining the desired endianness. Furthermore, it's vital to explicitly specify the Int16 binary format, which isn't mentioned in the commonly suggested command :

gdal_translate -of ENVI in.tif out.bin

Which can be done this way :

gdal_translate -ot Int16 -of ENVI in.tif out.bin

These commands fall short of addressing the problem, which is why I've provided tools below to overcome the encountered issues.

 

[kwerner]

I'm glad you were able to find a solution, and thank you so much for posting it. I hope it's able to help others in the future!

 

[Josh]

After dedicating several days to resolving this issue, I resorted to "reverse engineering" the WRF/WPS topography files (topo_30s). Using GDAL and Python, I've developed a compact script that can convert a GeoTIFF topography file into a binary format that's compatible with WRF/WPS. This script also generates the required index file. Please don't hesitate to reach out if you encounter any issues while using it. While the script works in my specific case, it may require some adjustments to ensure broader compatibility.
Note: the script cannot manage the creation of multiple tiles or the use of borders; it is therefore limited according to the WRF/WPS format to a single tile of maximum size 99999x99999.



Additionally, I'm providing a small script that enables you to visualize the topography based on the binary format used in WRF/WPS. You'll need to modify the information derived from the index file at the beginning of the script to make it work for your specific needs.

---

RESOLUTION NOTE :
The issue revolved around both the binary format output from GDAL and the data order. Notably, the origin of the conversion didn't align with the geographic lower-left corner. GDAL lacks a built-in option for determining the desired endianness. Furthermore, it's vital to explicitly specify the Int16 binary format, which isn't mentioned in the commonly suggested command :

gdal_translate -of ENVI in.tif out.bin

Which can be done this way :

gdal_translate -ot Int16 -of ENVI in.tif out.bin

These commands fall short of addressing the problem, which is why I've provided tools below to overcome the encountered issues.

Thanks Arty,

This has proved very helpful to me. As a new programmer, I was also unable to figure out how to use the the read_geogrid.c and write_geogrid.c programs.

I am using SRTM data with a resolution of 1 arc-second. I originally used GDAL command line utilities but after running geogrid, the terrain heights in my dataset didn't make sense. I was expecting values between 0 and 1300 and I ended up getting values over 10,000 and they appeared not orientated correctly. With your script, the values and orientation appears to be in line with what I expect, however, I haven't done much verification or ran WRF yet.

I was wondering why you rounded the values to 5 decimal places? Could this introduce errors in the exact placement of the cells based on the geographic location? Does geogrid only have a resolution for dx/dy of out to 5 decimal places? Honestly I don't know that much about geospatial information, but this was just a concern of mine. For example, a small rounding error propagating over the domain. Since I assume my dataset resolution is exactly 1 arc-second while dx/dy rounded to 5 decimals is an approximation?

Josh

 

[Arty]

Thanks Arty,

This has proved very helpful to me. As a new programmer, I was also unable to figure out how to use the the read_geogrid.c and write_geogrid.c programs.

I am using SRTM data with a resolution of 1 arc-second. I originally used GDAL command line utilities but after running geogrid, the terrain heights in my dataset didn't make sense. I was expecting values between 0 and 1300 and I ended up getting values over 10,000 and they appeared not orientated correctly. With your script, the values and orientation appears to be in line with what I expect, however, I haven't done much verification or ran WRF yet.

I was wondering why you rounded the values to 5 decimal places? Could this introduce errors in the exact placement of the cells based on the geographic location? Does geogrid only have a resolution for dx/dy of out to 5 decimal places? Honestly I don't know that much about geospatial information, but this was just a concern of mine. For example, a small rounding error propagating over the domain. Since I assume my dataset resolution is exactly 1 arc-second while dx/dy rounded to 5 decimals is an approximation?

Josh

Hello Josh,

Thanks for your post. I'm really glad this could help.

To answer your question - and for what I remember of it - I chose to round the data to 5 digits because I had trouble with some values length in my original GeoTiff dataset or while conversion that I needed to get regular in shape (i.e : avoid -17.3864° next to -17.38657°).

Also, even at the equator, since 1° ~= 111km, a precision of 0.00001° is close to the meter ; negligible approximation while working with kilo-metric grid cells (but I understand your concerns about large domains).

I must tell I did not have any trouble working with the topographic data I created through this process and am happy with the results. Feel free to try and add 1 supplementary decimal in the Py script if needed.