Saturday, 25 April 2015

gdal - How to remove strange gdaldem hillshade artifacts?


I've downloaded 9 chunks of 1 arc-second NED in IMG format from the National Map. I'm trying to create a hillshade from them.


I'm using:


GDAL 1.10.0, released 2013/04/24

First, I combine them and change the projection:


gdalwarp -t_srs EPSG:900913 n40w111/imgn40w111_1.img n40w112/imgn40w112_1.img \

n40w113/imgn40w113_1.img n41w111/imgn41w111_1.img n41w112/imgn41w112_1.img \
n41w113/imgn41w113_1.img n42w111/imgn42w111_1.img n42w112/imgn42w112_1.img \
n42w113/imgn42w113_1.img uinta-projected.tif

Then I create a hillshade from uinta-projected.tif:


gdaldem hillshade -compute_edges -co compress=lzw uinta-projected.tif uinta-hillshade.tif

However, when I take a look at my shiny new hillshade, it looks like this: ugly hillshade


Does anybody have any ideas what might be causing this grid artifact? I've tried using gdal_merge.py instead of gdalwarp and I end up with the same result.


gdalinfo uinta-projected.tif:



Driver: GTiff/GeoTIFF
Files: uinta-projected.tif
Size is 9253, 12173
Coordinate System is:
PROJCS["Google Maps Global Mercator",
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"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]

Origin = (-12579287.992128284648061,5161229.105774102732539)
Pixel Size = (36.130094040215752,-36.130094040215752)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-12579287.992, 5161229.106) (113d 0' 6.00"W, 42d11'34.83"N)
Lower Left (-12579287.992, 4721417.471) (113d 0' 6.00"W, 39d11'11.36"N)
Upper Right (-12244976.232, 5161229.106) (109d59'54.57"W, 42d11'34.83"N)

Lower Right (-12244976.232, 4721417.471) (109d59'54.57"W, 39d11'11.36"N)
Center (-12412132.112, 4941323.288) (111d30' 0.29"W, 40d42'24.63"N)
Band 1 Block=9253x1 Type=Float32, ColorInterp=Gray


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