Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

gdalwarp: "Unknown units=US survey foot" warning [GDAL 3.5.3] #6839

Closed
pinno opened this issue Dec 2, 2022 · 3 comments
Closed

gdalwarp: "Unknown units=US survey foot" warning [GDAL 3.5.3] #6839

pinno opened this issue Dec 2, 2022 · 3 comments
Assignees
@rouault

Comments

@pinno
Copy link

pinno commented Dec 2, 2022

Expected behavior and actual behavior.

gdalwarp (GDAL 3.5.3) is not working properly when resampling a GeoTiff raster in US survey foot units.

The same operation with the same input file works as expected with previous versions of GDAL, e.g. 3.0.4 and 3.1.4.

Steps to reproduce the problem.

gdalwarp input_raster.tif output_raster.tif -tr 1.0 1.0 -co compress=LZW -r max
Creating output file that is 11739P x 11705L.
Processing input_raster.tif [1/1] : 0Warning 1: Unknown units=US survey foot
...10...20...30...40...50...60...70...80...90...100 - done.

The output file is produced, but it is not as expected.

gdalinfo input_raster.tif gives the following output:

gdalinfo input_raster.tif 
Driver: GTiff/GeoTIFF
Files: input_raster.tif
Size is 17770, 17718
Coordinate System is:
COMPOUNDCRS["NAD83(2011) / Washington North (ftUS) + NAVD88 height (ftUS) - Geoid18 (ftUS)",
    PROJCRS["NAD83(2011) / Washington North (ftUS)",
        BASEGEOGCRS["NAD83(2011)",
            DATUM["NAD83 (National Spatial Reference System 2011)",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",6318]],
        CONVERSION["unnamed",
            METHOD["Lambert Conic Conformal (2SP)",
                ID["EPSG",9802]],
            PARAMETER["Latitude of false origin",47,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8821]],
            PARAMETER["Longitude of false origin",-120.833333333333,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8822]],
            PARAMETER["Latitude of 1st standard parallel",48.7333333333333,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8823]],
            PARAMETER["Latitude of 2nd standard parallel",47.5,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8824]],
            PARAMETER["Easting at false origin",1640416.667,
                LENGTHUNIT["US survey foot",0.304800609601219],
                ID["EPSG",8826]],
            PARAMETER["Northing at false origin",0,
                LENGTHUNIT["US survey foot",0.304800609601219],
                ID["EPSG",8827]]],
        CS[Cartesian,2],
            AXIS["easting",east,
                ORDER[1],
                LENGTHUNIT["US survey foot",0.304800609601219]],
            AXIS["northing",north,
                ORDER[2],
                LENGTHUNIT["US survey foot",0.304800609601219]],
        ID["EPSG",6597]],
    VERTCRS["NAVD88 height (ftUS)",
        VDATUM["North American Vertical Datum 1988"],
        CS[vertical,1],
            AXIS["gravity-related height",up,
                LENGTHUNIT["US survey foot",0.304800609601219]],
        ID["EPSG",6360]]]
Data axis to CRS axis mapping: 1,2,3
Origin = (1264367.180066652828828,231922.768387514282949)
Pixel Size = (0.660608175138232,-0.660608175138232)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 1264367.180,  231922.768) (122d21'30.26"W, 47d37'32.13"N)
Lower Left  ( 1264367.180,  220218.113) (122d21'26.87"W, 47d35'36.63"N)
Upper Right ( 1276106.187,  231922.768) (122d18'38.91"W, 47d37'34.39"N)
Lower Right ( 1276106.187,  220218.113) (122d18'35.63"W, 47d35'38.89"N)
Center      ( 1270236.684,  226070.441) (122d20' 2.92"W, 47d36'35.52"N)
Band 1 Block=17770x1 Type=Byte, ColorInterp=Gray
  Unit Type: US survey foot

Operating system

Ubuntu 22.04 64 bit

GDAL version and provenance

GDAL 3.5.3, installed with conda install -c conda-forge gdal
@rouault
Copy link
Member

rouault commented Dec 2, 2022

I can't reproduce with a file when assigning to it the CRS output by gdalinfo. I suspect there's something very specific to the GeoTIFF encoding of your exact input raster. Please provide a link where it can be downloaded, or the output of "listgeo input_raster.tif"
If you don't have the listgeo utility, it is provided with libgeotiff.
Or run docker run --rm -it -v $PWD:$PWD osgeo/gdal:alpine-small-3.5.3 sh -c "apk add libgeotiff; listgeo $PWD/input_raster.tif"

@pinno
Copy link
Author

pinno commented Dec 2, 2022

The output of listgeo input_raster.tif is:

Geotiff_Information:
   Version: 1
   Key_Revision: 1.1
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                 0                 0                
         1264367.18006665  231922.768387514  0                
      ModelPixelScaleTag (1,3):
         0.660608175138232 0.660608175138232 1                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,78): "NAD83(2011) / Washington North (ftUS) + NAVD88 height (ftUS) - Geoid18 (ftUS)"
      ProjectedCRSGeoKey (Short,1): Code-6597 (NAD83(2011) / Washington North (ftUS))
      VerticalGeoKey (Short,1): Code-6360 (NAVD88 height (ftUS))
      End_Of_Keys.
   End_Of_Geotiff.

PCS = 6597 (NAD83(2011) / Washington North (ftUS))
Projection = 15367 (SPCS83 Washington North zone (US Survey feet))
Projection Method: CT_LambertConfConic_2SP
   ProjFalseOriginLatGeoKey: 47.000000 ( 47d 0' 0.00"N)
   ProjFalseOriginLongGeoKey: -120.833333 (120d50' 0.00"W)
   ProjStdParallel1GeoKey: 48.733333 ( 48d44' 0.00"N)
   ProjStdParallel2GeoKey: 47.500000 ( 47d30' 0.00"N)
   ProjFalseEastingGeoKey: 500000.000102 m
   ProjFalseNorthingGeoKey: 0.000000 m
GCS: 6318/NAD83(2011)
Datum: 1116/NAD83 (National Spatial Reference System 2011)
Ellipsoid: 7019/GRS 1980 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
Projection Linear Units: 9003/US survey foot (0.304801m)

Corner Coordinates:
Upper Left    ( 1264367.180,  231922.768)  (122d21'30.26"W, 47d37'32.13"N)
Lower Left    ( 1264367.180,  220218.113)  (122d21'26.87"W, 47d35'36.63"N)
Upper Right   ( 1276106.187,  231922.768)  (122d18'38.91"W, 47d37'34.39"N)
Lower Right   ( 1276106.187,  220218.113)  (122d18'35.63"W, 47d35'38.89"N)
Center        ( 1270236.684,  226070.441)  (122d20' 2.92"W, 47d36'35.52"N)

@rouault
Copy link
Member

rouault commented Dec 2, 2022

ah ok, I can actually reproduce the warning when using gdalwarp (I used gdal_translate previously).

@rouault rouault self-assigned this Dec 2, 2022
rouault added a commit to rouault/gdal that referenced this issue Dec 2, 2022
@rouault
gdalwarp: fix issue with vertical shift, in particular when CRS has U…
4d5fe76 
…S survey foot as vertical unit (fixes OSGeo#6839)
rouault added a commit to rouault/gdal that referenced this issue Dec 2, 2022
@rouault
gdalwarp: fix issue with vertical shift, in particular when CRS has U…
68d7c68 
…S survey foot as vertical unit (fixes OSGeo#6839)
rouault added a commit to rouault/gdal that referenced this issue Dec 2, 2022
@rouault
gdalwarp: fix issue with vertical shift, in particular when CRS has U…
9d85202 
…S survey foot as vertical unit (fixes OSGeo#6839)
rouault added a commit that referenced this issue Dec 3, 2022
@rouault
Merge pull request #6840 from rouault/fix_gdalwarp_vertical_shift
7daa0ab 
gdalwarp: fix issue with vertical shift, in particular when CRS has US survey foot as vertical unit (fixes #6839)
rouault added a commit that referenced this issue Dec 3, 2022
@rouault @github-actions
gdalwarp: fix issue with vertical shift, in particular when CRS has U…
b1be1ab 
…S survey foot as vertical unit (fixes #6839)
rouault added a commit that referenced this issue Dec 3, 2022
@rouault
Merge pull request #6847 from OSGeo/backport-6840-to-release/3.6
17fe94c 
[Backport release/3.6] gdalwarp: fix issue with vertical shift, in particular when CRS has US survey foot as vertical unit (fixes #6839)
@rouault rouault closed this as completed Dec 3, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@rouault rouault

Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@rouault @pinno