Re: [gdal-dev] USGS 3DEP (3D Elevation Program) - feet should be metres, how can I fix this?

Fri, 22 Mar 2024 01:57:13 -0700 

Hi,
having looked at the file, I believe GDAL behaves correctly given the metadata in the file. This is clearly a unit type attached to the Band, and it should be corrected by the data producer to reflect what is really used in it 

One way to correct it on your side is for example this little Python snippet

from osgeo import gdal
ds = gdal.Open('test.tif')
ds.GetRasterBand(1).SetUnitType('metre')


As the file is used in read-only mode, this will generate an auxiliary 
.aux.xml side car file overriding the unit.


You can also directly override the wrong metadata inside the TIFF with

ds = gdal.OpenEx('test.tif', open_options=['IGNORE_COG_LAYOUT_BREAK=YES'])
ds.GetRasterBand(1).SetUnitType('metre')


The IGNORE_COG_LAYOUT_BREAK=YES open option is because modifying TIFF 
metadata causes it to be rewritten at the end of the file, which 
"breaks" the COG convention. But that's still a perfectly valid GeoTIFF 
file, and will not change anything for local use of such file. You could 
reprocess that to a fully compliant COG with gdal_translate -of COG



Regarding how that is used by gdalwarp, the doc enhancement at 
https://github.com/OSGeo/gdal/pull/9527 captures my initial attempt at 
answering your questions.


Even


Le 22/03/2024 à 01:16, Conrad Bielski a écrit :

Just realised I should have put you also in CC

Here is my response message again.

Thanks Even,

you are correct about the tag number (I sent the incorrrect one).


I have provided a sample tif and some other descriptions of the issue 
here:
bad_geotiff - Google Drive 
<https://drive.google.com/drive/folders/1LVJqoKrnLRf6_5FpE7a5SrUGIWBPQmqa?usp=sharing>


        


    bad_geotiff - Google Drive


<https://drive.google.com/drive/folders/1LVJqoKrnLRf6_5FpE7a5SrUGIWBPQmqa?usp=sharing> 



The three files there are the following:
1.  A bad tiff (this is slightly smaller; they are all pretty big).
2.  A GDAL info dump (survey feet is almost at the bottom)
3. A PDF with an explanation of the issue

thanks for your help.
Conrad



On Thursday, March 21, 2024 at 04:49:57 PM EDT, Even Rouault 
<even.roua...@spatialys.com> wrote:




Le 21/03/2024 à 21:45, Conrad Bielski via gdal-dev a écrit :
Hello GDALers,


I have a question about reading USGS 3DEP (3D Elevation Program) data. 
Inside of this data, a GEOTIFF tag 42114 is provided which is causing 
problems with datum shifts.


There's no such thing as a GEOTIFF tag 42114.

The closest tags are:

#define TIFFTAG_GDAL_METADATA 42112
#define TIFFTAG_GDAL_NODATA 42113


It would be helpuful if you could provide the file, or the result of a 
"tiffdump -m 1000 the.tif" and "listgeo the.tif" on it.



And was the .tif.aux.xml provided with the TIFF, or computed by GDAL. 
Because the information might come from it too.






So when I use GDAL to compute the datum shifts, the tag is read and 
interprets that the DEM is showing elevation in 'feet' while the DEM 
is actually in metres. The DEMs are in fact in meter elevations and 
meter UTM horizontal coordinates. This is obviously erroneously 
integrated into the tag.



So my question is whether anyone has come across this issue and found 
a solution? Is there a way to edit the offending tag so that it is 
correctly interpreted as metres instead of feet? Other potential 
solutions so that GDAL interprets the elevation correctly?



Thanks in advance for your help. I provide below the output of 
gdalinfo which I hope could help with the offending info in bold.


Have a great day,
Conrad




Driver: GTiff/GeoTIFF

Files: 
I:\v3_us_3dep\asheville\ashevilleUSGS_1M_17_x34y392_NC_Phase5_2018_A18.tif

 
I:\v3_us_3dep\asheville\ashevilleUSGS_1M_17_x34y392_NC_Phase5_2018_A18.tif.aux.xml
Size is 10012, 10012
Coordinate System is:
PROJCRS["NAD83 / UTM zone 17N",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["UTM zone 17N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-81,
ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],

        AREA["North America - between 84°W and 78°W - onshore and 
offshore. Canada - Nunavut; Ontario; Quebec. United States (USA) - 
Florida; Georgia; Kentucky; Maryland; Michigan; New York; North 
Carolina; Ohio; Pennsylvania; South Carolina; Tennessee; Virginia; 
West Virginia."],

        BBOX[23.81,-84,84,-78]],
    ID["EPSG",26917]]
Data axis to CRS axis mapping: 1,2
Origin = (339993.999958705157042,3920005.999981591943651)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
  LAYOUT=COG
  PREDICTOR=3
Corner Coordinates:
Upper Left  (  339994.000, 3920006.000) ( 82d45'43.98"W, 35d24'38.21"N)
Lower Left  (  339994.000, 3909994.000) ( 82d45'36.92"W, 35d19'13.36"N)
Upper Right (  350006.000, 3920006.000) ( 82d39' 7.17"W, 35d24'43.82"N)
Lower Right (  350006.000, 3909994.000) ( 82d39' 0.55"W, 35d19'18.95"N)
Center      (  345000.000, 3915000.000) ( 82d42'22.15"W, 35d21'58.63"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  Min=669.509 Max=1604.949
  Minimum=669.509, Maximum=1604.949, Mean=988.295, StdDev=159.804
  NoData Value=-999999
  Overviews: 5006x5006, 2503x2503, 1251x1251, 625x625, 312x312
*Unit Type: US survey foot*
  Metadata:
    STATISTICS_APPROXIMATE=YES
    STATISTICS_MAXIMUM=1604.9490966797
    STATISTICS_MEAN=988.29463693706
    STATISTICS_MINIMUM=669.50915527344
    STATISTICS_STDDEV=159.80361228109
    STATISTICS_VALID_PERCENT=100


_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org  <mailto:gdal-dev@lists.osgeo.org>
https://lists.osgeo.org/mailman/listinfo/gdal-dev  
<https://lists.osgeo.org/mailman/listinfo/gdal-dev>

--
http://www.spatialys.com  <http://www.spatialys.com>
My software is free, but my time generally not.


--
http://www.spatialys.com
My software is free, but my time generally not.

_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev






















					Reply via email to