Losing resolution of my DEM/TIFF when converting from HGT

Question

I need to get a nice hillshade from this area. My steps:

1) Get HGT files from the area: 
 -- phyghtmap --download-only --a -43.7544:-23.2363:-42.0378:-22.3183 --earthdata-user=myuser --earthdata-password=mypass --viewfinder-mask=1

2) Create a VRT file to join all downloaded HGT:
  -- gdalbuildvrt ./test.vrt hgt/SRTM1v3.0/S23W043.hgt  hgt/SRTM1v3.0/S23W044.hgt  hgt/SRTM1v3.0/S24W043.hgt  hgt/SRTM1v3.0/S24W044.hgt

3) Create the TIF:
  -- gdaldem hillshade test.vrt test.tif -z 10 -s 90000

The result is very ugly:

What I'm doing wrong?

EDIT

This is my GDALINFO result:

Driver: GTiff/GeoTIFF
Files: test.tif
Size is 7201, 7201
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-44.000138888888891,-21.999861111111109)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -44.0001389, -21.9998611) ( 44d 0' 0.50"W, 21d59'59.50"S)
Lower Left  ( -44.0001389, -24.0001389) ( 44d 0' 0.50"W, 24d 0' 0.50"S)
Upper Right ( -41.9998611, -21.9998611) ( 41d59'59.50"W, 21d59'59.50"S)
Lower Right ( -41.9998611, -24.0001389) ( 41d59'59.50"W, 24d 0' 0.50"S)
Center      ( -43.0000000, -23.0000000) ( 43d 0' 0.00"W, 23d 0' 0.00"S)
Band 1 Block=7201x1 Type=Int16, ColorInterp=Gray
    Computed Min/Max=-38.000,2277.000
  NoData Value=-32768

EDIT 2: Using gdaldem hillshade test.vrt test.tif -z 1 -s 80000

... and this is my best result so far. After load the tif as a Geoserver layer. You can see the very poor "pixel as big sqares" result:

...same image from https://www.opencyclemap.org/ ( Pão de Açúcar, Rio de Janeiro, RJ, Brasil):

EDIT 3: Making progress:

gdal_translate -tr 0.00009 -0.00009 -r cubicspline -of GTiff test.vrt test.tif

gdaldem hillshade -co TILED=YES -co compress=lzw -s 111120 -z 1 -az 315 -alt 60 -compute_edges test.tif final.tif

It looks like your DEM is bumpy and needs to be smoothed, though looking nice isn't a quantitative term, you could try resample to 5 times the cell size using bilinear then back to the original cell size with the same resample method to smooth the raster prior to creating the hillshade. Just remember the DEM is now inaccurate, you will loose the peaks and dips, so don't do anything with it except for hillshade. — Michael Stimson, Jun 02 '17 at 02:45
I'm a complete newbie on this. I'll appreciate some code examples. I'll consider as nice something like this:http://blog.mastermaps.com/2012/06/creating-hillshades-with-gdaldem.html or this https://i.stack.imgur.com/b6i3y.jpg — Magno C, Jun 02 '17 at 02:48
After you've created the VRT use GDAL_TRANSLATE http://www.gdal.org/gdal_translate.html with -tr XRes YRes -r bilinear to convert the VRT into a single image (I would recommend -of HFA or GTIFF) to resample the raster. It's not clear what your original cell size is but you can get it with GDALINFO http://www.gdal.org/gdalinfo.html — Michael Stimson, Jun 02 '17 at 02:52
Your pixel size is 0.000277777777778 so a -tr 0.001389 0.001389 will resample to 5 times the cell size, you might need to go higher than 5 so it might be a good idea to try only one of your downloaded images resampling coarse and back to fine until you're happy with the output then run it on the whole VRT. — Michael Stimson, Jun 02 '17 at 03:03
Can't find -tr option in my gdal_translate. Didn't you mean http://www.gdal.org/gdalwarp.html ? — Magno C, Jun 02 '17 at 03:09
-tr (starting from gdal 2.0) set target resolution. The values must be expressed in georeferenced units. Both must be positive values. This is mutually exclusive with -outsize and -a_ullr. What version is your GDAL? You might need to get a more recent version; multiple GDAL installs can happily coexist, I use 2.1 x64 and am quite happy with it but have older versions installed for dependent software. In the older versions you could still use -outsize. — Michael Stimson, Jun 02 '17 at 03:13
Same poor results. Trying now with gdaldem hillshade test.vrt test.tif -z 1 -s 80000 and got better result but far from what I need. See my edit. Seems I need to play with z and s parameters but have no idea what I'm doing... — Magno C, Jun 02 '17 at 13:59
This may help : https://gis.stackexchange.com/questions/144535/creating-transparent-hillshade — Magno C, Jun 02 '17 at 16:03
Or this: https://gis.stackexchange.com/questions/12833/smoothing-dem-using-grass — Magno C, Jun 02 '17 at 18:00
@Michael Miles-Stimson please tell me more about the pixel size. I'm playing with the -tr XRes YRes and discovered lower is better. changing to a small value like 0.00005 seems to improve the resolution but the file size become huge. — Magno C, Jun 03 '17 at 23:35

Magno C · Accepted Answer · 2017-06-04T03:01:25.690

Well, after my EDIT 3 I found this post: How to antialiase tiles when seeding a layer from an GeoTiff in GeoServer

In GeoServer under the point WMS you can activate antialiasing. This was already checked but the raster rendering option was nearest neighbor. I switched it to bilinear or bicubic and now the resulting tiles are nice and smooth looking. (User: strangeoptics)

I changed it to bilinear and now I can see this nice result. I've noticed some horizontal lines and found no way do get rid of them:

The final setting is:

1) Get the HGT files and create contour lines to import to OSM:
    phyghtmap --pbf --srtm=1 --a -43.7544:-23.2363:-42.0378:-22.3183 --earthdata-user=myuser --earthdata-password=mypassword

2) Create VRT from the HGT files:
    gdalbuildvrt ./teste.vrt hgt/SRTM1v3.0/S23W043.hgt  hgt/SRTM1v3.0/S23W044.hgt  hgt/SRTM1v3.0/S24W043.hgt  hgt/SRTM1v3.0/S24W044.hgt

3) Do the abracadabra: 
    gdal_translate -tr 0.000050 0.000050 -r cubicspline -of GTiff test.vrt test.tif

4) And now do some kung-fu-code:
    gdaldem hillshade -co TILED=YES -co compress=lzw -s 111120 -z 5 -combined -compute_edges test.tif final.tif

5) Import contour lines to OSM. I prefer to create a separate database and give a small `.style` to import just the needed columns.
    osm2pgsql --verbose --create --style ./srtm.style --database contour --username postgres -W --host 127.0.0.1 lon-43.23_-43.05lat-23.00_-22.90_srtm1v3.0.osm.pbf
    osm2pgsql --verbose --append <all other pbf files> (be careful with --append parameter)

This is the style file used:

srtm.style
# OsmType Tag DataType Flags
node,way contour text linear
node,way contour_ext text linear
node,way ele text linear

Applied this style to the raster layer (GeoTiff Coverage Store pointed to final.tif):

<?xml version="1.0" encoding="UTF-8"?>
<sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">
  <sld:UserLayer>
    <sld:LayerFeatureConstraints>
      <sld:FeatureTypeConstraint/>
    </sld:LayerFeatureConstraints>
    <sld:UserStyle>
      <sld:Title/>
      <sld:FeatureTypeStyle>
        <sld:Name>name</sld:Name>
        <sld:FeatureTypeName>Feature</sld:FeatureTypeName>
        <sld:Rule>
          <sld:MinScaleDenominator>5000</sld:MinScaleDenominator>
          <sld:RasterSymbolizer>
            <sld:Geometry>
              <ogc:PropertyName>grid</ogc:PropertyName>
            </sld:Geometry>
            <sld:ColorMap>
                <ColorMapEntry color="#000000" quantity="0.0" label="low" opacity="0.0"/>
                <ColorMapEntry color="#999999" quantity="1.0" label="high" opacity="0.75"/>
                <ColorMapEntry color="#FFFFFF" quantity="256" label="alpha" opacity="0.75"/> 
            </sld:ColorMap>
          </sld:RasterSymbolizer>
        </sld:Rule>
        <sld:VendorOption name="composite">multiply</sld:VendorOption>
      </sld:FeatureTypeStyle>
    </sld:UserStyle>
  </sld:UserLayer>
</sld:StyledLayerDescriptor>

Do not ask me anything because actualy I have no idea what I've done. Just getting ideas from thousands sources and adjusting to fit my needs. BTW: The result file is very huge but seems Geoserver knows how to handle it.

Finally I throw the annoying horizontal lines away with -tr 0.000150 0.000150 in gdal_translate. — Magno C, Jun 04 '17 at 02:51

Losing resolution of my DEM/TIFF when converting from HGT

1 Answers1

Linked