Geotiff DEM problem: NULL not recognized

[Markus Neteler]

Frank,

currently I am teaching GRASS and GDAL/OGR at university, and came across a 
NULL problem with US DEM data for Spearfish converted from ArcBinGRID to
GeoTIFF:

gdalinfo -mm d44103f7.tif
Driver: GTiff/GeoTIFF
Size is 337, 467
Coordinate System is:
PROJCS["NAD83 / UTM zone 13N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-105],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26913"]]
Origin = (589024.160521,4955929.485738)
Pixel Size = (30.000000,-30.000000)
Corner Coordinates:
Upper Left  (  589024.161, 4955929.486) (103d52'31.18"W, 44d45'4.38"N)
Lower Left  (  589024.161, 4941919.486) (103d52'39.96"W, 44d37'30.42"N)
Upper Right (  599134.161, 4955929.486) (103d44'51.48"W, 44d44'59.59"N)
Lower Right (  599134.161, 4941919.486) (103d45'1.25"W, 44d37'25.66"N)
Center      (  594079.161, 4948924.486) (103d48'45.97"W, 44d41'15.08"N)
Band 1 Block=337x12 Type=Int16, ColorInterp=Gray
    Computed Min/Max=-32768.000,1041.000

It seems that NULL is not recognized properly, subsequently r.in.gdal
fails and also gdal_merge.py does not produce useful results.

The original data was an ARC grid coverage:

 gdalinfo d44103f7/hdr.adf
Driver: AIG/Arc/Info Binary Grid
Size is 337, 467
Coordinate System is:
PROJCS["UTM Zone 13, Northern Hemisphere",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AXIS["Lat",NORTH],
        AXIS["Long",EAST],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-105],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["METERS",1]]
Origin = (589024.160521,4955929.485738)
Pixel Size = (30.000000,-30.000000)
Corner Coordinates:
Upper Left  (  589024.161, 4955929.486) (103d52'31.18"W, 44d45'4.38"N)
Lower Left  (  589024.161, 4941919.486) (103d52'39.96"W, 44d37'30.42"N)
Upper Right (  599134.161, 4955929.486) (103d44'51.48"W, 44d44'59.59"N)
Lower Right (  599134.161, 4941919.486) (103d45'1.25"W, 44d37'25.66"N)
Center      (  594079.161, 4948924.486) (103d48'45.97"W, 44d41'15.08"N)
Band 1 Block=256x4 Type=Int16, ColorInterp=Undefined
  Min=896.000 Max=1051.000
  NoData Value=-32768

So it seems that GDAL is "forgetting" about the NULL value when
doing:
gdal_translate -of GTiff -a_srs '+init=epsg:26913' d44103f7/hdr.adf d44103f7.tif

However, when I do:
gdal_translate -of GTiff  d44103f7/hdr.adf d44103f7.tif
everything is fine. Yes, ok, I just realized that -a_srs '+init=epsg:26913' 
isn't necessary. But it should not eat the NULL value :-)

The Arc coverage file comes from
http://www.sdgs.usd.edu/register/index.html
"Digital Elevation Model (DEM)"

I don't know yet if -a_srs always deletes the NULL value and just want
to notify you.

Best regards,

 Markus

PS: I'll attach the file to the report.

[dron]

Frank,

I have fixed this in gdal_translate.cpp. NoData value now successfully copied
over the virtual dataset.

[warmerdam]

Cool, thanks Andrey.