Monday 29 May 2017

hillshade - Losing resolution of my DEM/TIFF when converting from HGT


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:


enter image description here


enter image description here


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


enter image description here


... 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:



enter image description here


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


enter image description here


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

enter image description here



Answer




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:


enter image description here enter image description here


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 (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):












name
Feature

5000


grid









multiply






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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...