anoved / phstl

Convert GDAL rasters (like GeoTIFF heightmaps) to 3D STL surfaces.
MIT License
73 stars 21 forks source link

Issue with using Degree Units: Breaks on Certain Tiffs #10

Closed Skylion007 closed 8 years ago

Skylion007 commented 8 years ago

Hello,

Excellent script, it proved really helpful. I did notice one bug that caused me about three hours headache yesterday though. Mainly that while this repo seems to works perfectly if you are using a UTM coordinate system in your rasters, if use one that relies on degrees in the Geotransform (such as WGS84 if I am not mistaken) then it fails horrendously. The z-scale gets divided by the xyres variable which with my particular Geotiff was in the order of magnitude of 10^-6. This meant that the values in my geotiff went from having about a hundred meters of differences to 10^6 meters of difference, which really screwed up the mesh.

I fixed it in my version of the script by removing xyres variable entirely, but that only works for my particular geotif since the height and x y resolution are the exact same. I am not sure how to fix this in a more generic fashion, but I thought I would bring the issue to your attention and anyone else who finds this otherwise fabulous repo.

Happy bug hunting and thanks for the helpful comments!

anoved commented 8 years ago

Thanks for the feedback! I'm happy to hear you find the script useful.

It's been a while since I've looked at the code, but I believe it does assume the same units for elevation as for xy pixel resolution. Works fine with UTM/meter data but would indeed go haywire if xy coordinates were degrees. Good catch.

I may be able to better utilize metadata in the Geotiff to recognize if the units differ. Anyway, thanks again for bringing this to my attention.

anoved commented 8 years ago

Any chance you could provide a sample of your GeoTIFFs or link to similar source? Sorry for the delay but I still intend to fix this better. Thanks.

anoved commented 8 years ago

After reviewing some samples and discussion of similar datasets I recommend projecting the original dataset first, since with geographic (lat/lon) coordinates scale inherently varies with latitude, precluding a trivial conversion between xy and z units.

If that's not an option, you can fudge it. Notes in the gdaldem documentation provide a starting point: for equatorial datasets, you can use a scale of 111120 vertical units per horizontal unit to bring meters and lat/lon degrees into the same ballpark. Use the new -s 111120 option or use the reciprocal with the existing z scale option: -z 0.000008999.