Opened 12 years ago

Closed 9 years ago

#4677 closed defect (fixed)

unspecified projection parameters get nan in geotiff output

Reported by: kyngchaos Owned by: warmerdam
Priority: highest Milestone:
Component: GDAL_Raster Version: 1.9.1
Severity: major Keywords: GTIFF
Cc:

Description

I'm warping a raster to:

"+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +ellps=GRS80 +datum=NAD83"

The unspecified false easting and northing are now getting set to nan in the geotiff output, instead of the default 0. Also, the default units are wrong.

PROJCS["unnamed",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",23],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting","nan"],
    PARAMETER["false_northing","nan"],
    UNIT["degree minute second hemisphere",0,
        AUTHORITY["EPSG","9108"]]]

ogr2ogr works in a direct projection change, but gdal_contour does propogate the nans from a warped geotiff to the shapefile output.

This is new in 1.9.1. 1.9.0 and previous never did this.

Change History (9)

comment:1 by warmerdam, 12 years ago

Keywords: GTIFF added
Priority: normalhighest
Status: newassigned

I try:

gdalwarp -t_srs '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +ellps=GRS80 +datum=NAD83' utm.tif out.tif

The listgeo report on the output file looks like:

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                 0                 0                
         -1976601.56760975 1429332.91862047  0                
      ModelPixelScaleTag (1,3):
         60.022242806635   60.022242806635   0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,8): "unnamed"
      GeographicTypeGeoKey (Short,1): GCS_NAD83
      GeogCitationGeoKey (Ascii,6): "NAD83"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogInvFlatteningGeoKey (Double,1): 298.257222101    
      GeogTOWGS84GeoKey (Double,3): 0                0                0                
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_AlbersEqualArea
      ProjLinearUnitsGeoKey (Short,1): Angular_DMS_Hemisphere
      ProjStdParallel1GeoKey (Double,1): 29.5             
      ProjStdParallel2GeoKey (Double,1): 45.5             
      ProjNatOriginLongGeoKey (Double,1): -96              
      ProjNatOriginLatGeoKey (Double,1): 23               
      ProjFalseEastingGeoKey (Double,1): 0                
      ProjFalseNorthingGeoKey (Double,1): 0                
      End_Of_Keys.
   End_Of_Geotiff.

which has reasonable seeming false easting/northing (but an odd choice of linear units). But the gdalinfo report is:

PROJCS["unnamed",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",23],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting","-nan"],
    PARAMETER["false_northing","-nan"],
    UNIT["degree minute second hemisphere",0,
        AUTHORITY["EPSG","9108"]]]

Hmm....

Investigating.

comment:2 by warmerdam, 12 years ago

Ah, it is actually the linear units that is screwing up the false easting/northing when the values are transformed to/from meters.

comment:3 by warmerdam, 12 years ago

Also in trunk...

comment:4 by warmerdam, 12 years ago

The problem is this code in GTIFSetfromOGCDefn():

    if( poSRS->GetAuthorityName("PROJCS|UNIT") != NULL 
        && EQUAL(poSRS->GetAuthorityName("PROJCS|UNIT"),"EPSG") )
        nUOMLengthCode = atoi(poSRS->GetAuthorityCode("PROJCS|UNIT"));

If there is no PROJCS UNIT object, the search continues down into the GEOGCS and finds it's angular units. This didn't use to happen because we didn't save WKT defined units in GeoTIFF. The issue it tickled in this case because the PROJ.4 translator does not set linear units in the WKT if they aren't set in the PROJ.4 string.

comment:5 by warmerdam, 12 years ago

I have applied a minimal fix in 1.9 branch to avoid this error (r24496).

I would like to make things substantially more robust in trunk.

comment:6 by kyngchaos, 12 years ago

That works applied to the 1.9.1 release - I now have the correct units and 0,0 easting/northing.

comment:7 by Even Rouault, 12 years ago

I've just added a regression test in trunk (r24502) and branches/1.9 (r24503). Letting Frank design his solution for trunk.

comment:8 by Even Rouault, 12 years ago

I've forward-ported the workaround in trunk, because I've just hit that issue when doing a gdalwarp from a BSB to a GTIFF, and that might also annoy other people.

r24790 /trunk/gdal/frmts/gtiff/gt_wkt_srs.cpp: GTiff: minimal change to avoid picking up geog units as projcs units (#4677, forward port of r24496)

Still leaving the ticket open if Frank has a better fix.

comment:9 by Jukka Rahkonen, 9 years ago

Resolution: fixed
Status: assignedclosed

Two years is longer than the reaction speed of Frank. Closing because I believe no better fix is in a queue.

Note: See TracTickets for help on using tickets.