[gdal-dev] gdalinfo coordinates problem?

Andre Joost andre+joost at nurfuerspam.de
Sat Jan 18 23:54:16 PST 2014 

Am 18.01.2014 21:52, schrieb Hermann Peifer:
>
> As far as I can see: lon_0 is the same in both cases and there is also
> no actual difference between +datum=NAD83 (gdal) and +ellps=GRS80
> +towgs84=0,0,0,0,0,0,0 (qgis).
>
> I am not sure where +x_0=600000 (qgis) comes from, but gdal's +x_0 value
> is simply 6458320.416666666 * 0.3048006096012192 = 1968500 *meters*, so
> in fact only +units=us-ft is wrong there.

No, that's the only thing right. The DEM zip file includes a html and a 
xml file which gives the intended SRS information:

Geodetic_Model:

     Horizontal_Datum_Name: North American Datum of 1983
     Ellipsoid_Name: Geodetic Reference System 80
     Semi-major_Axis: 6378137.000000
     Denominator_of_Flattening_Ratio: 298.257222

cordsysn:

     geogcsn: GCS_North_American_1983
     projcsn: NAD_1983_StatePlane_Pennsylvania_South_FIPS_3702_Feet

Planar:

     Map_Projection:

         Map_Projection_Name: Lambert Conformal Conic
         Lambert_Conformal_Conic:

             Standard_Parallel: 39.933333
             Standard_Parallel: 40.966667
             Longitude_of_Central_Meridian: -77.750000
             Latitude_of_Projection_Origin: 39.333333
             False_Easting: 1968500.000000
             False_Northing: 0.000000

     Planar_Coordinate_Information:

         Planar_Coordinate_Encoding_Method: row and column
         Planar_Distance_Units: survey feet

That should be equivalent to the ESRI definition called EPSG:102790 in 
QGIS 2.0:

+proj=lcc +lat_1=39.93333333333333 +lat_2=40.96666666666667 
+lat_0=39.33333333333334 +lon_0=-77.75 +x_0=600000.0000000001 +y_0=0 
+datum=NAD83 +units=us-ft +no_defs

with that, the DEM file will place in the extent given by the same html 
file:

     West_Bounding_Coordinate: -75.601950
     East_Bounding_Coordinate: -75.565369
     North_Bounding_Coordinate: 40.054419
     South_Bounding_Coordinate: 40.026303

and the upper left corner reflects the name of the DEM

The only fault of QGIS is that it guesses the CRS wrongly as EPSG:3219:

+proj=lcc +lat_1=40.96666666666667 +lat_2=39.93333333333333 
+lat_0=39.33333333333334 +lon_0=-77.75 +x_0=600000 +y_0=0 +ellps=GRS80 
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs

which has a different lat_1 and different units.

EPSG:2272 has the right units, the wrong lat_1, but places the file in 
the same location:

+proj=lcc +lat_1=40.96666666666667 +lat_2=39.93333333333333 
+lat_0=39.33333333333334 +lon_0=-77.75 +x_0=600000 +y_0=0 +ellps=GRS80 
+towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

x_0=600000 is the right value for the CRS. Creating a vector grid on 
that CRS, the central meridian gets the value of 1968500 feet as 
intended in the html document. You get x_0 by calculating 196850*0.3048

By the way, 600000 is the false easting in metres, as used by EPSG:3219, 
because both PA systems have their origin in the same place (which makes 
some sense).

gdalsrsinfo creates wrong x_0 for the project string, and wrong 
false_easting in OGC_WKT.

Greetings,
André Joost


>
> On 2014-01-18 20:09, Etienne Tourigny wrote:
>> I noticed that qgis-2.0 also has a problem with this file, but places it
>> around 54 degrees west longitude (if using on-the-fly reprojection to
>> WGS84).
>>
>> The CRS is different from that identified by gdal (the lon_0, units and
>> datum are different)
>>
>> qgis:
>> +proj=lcc +lat_1=40.96666666666667 +lat_2=39.93333333333333
>> +lat_0=39.33333333333334 +lon_0=-77.75 +x_0=600000 +y_0=0 +ellps=GRS80
>> +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
>>
>>
>> gdal:
>>
>> PROJ.4 : '+proj=lcc +lat_1=40.96666666666667 +lat_2=39.93333333333333
>> +lat_0=39.33333333333334 +lon_0=-77.75 +x_0=1968500 +y_0=0 +datum=NAD83
>> +units=us-ft +no_defs '
>>
>> OGC WKT :
>> PROJCS["NAD_1983_StatePlane_Pennsylvania_South_FIPS_3702_Feet",
>> 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["Lambert_Conformal_Conic_2SP"],
>> PARAMETER["standard_parallel_1",40.96666666666667],
>> PARAMETER["standard_parallel_2",39.93333333333333],
>> PARAMETER["latitude_of_origin",39.33333333333334],
>> PARAMETER["central_meridian",-77.75],
>> PARAMETER["false_easting",6458320.416666666],
>> PARAMETER["false_northing",0],
>> UNIT["us_survey_feet",0.3048006096012192],
>> AUTHORITY["EPSG","32129"]]
>>
>>
>>
>>
>> On Fri, Jan 17, 2014 at 8:37 PM, Even Rouault
>> <even.rouault at mines-paris.org <mailto:even.rouault at mines-paris.org>>
>> wrote:
>>
>> Le vendredi 17 janvier 2014 21:28:35, Reynolds, Scott a écrit :
>> > Hi,
>> >
>> > I've downloaded several DEMs (e.g 27002570PAS_dem.zip) from
>> >
>> ftp://pamap.pasda.psu.edu/pamap_lidar/cycle1/DEM/South/2008/20000000/.
>> > Gdalinfo is reporting the corner coordinates to be in the
>> vicinity of 91
>> > degrees 31 minutes west longitude which should be more like 75
>> degrees 36
>> > minutes west longitude. Can someone explain what is happening here?
>>
>> GeoTIFF support for non-meter linear units has been a major source of
>> headaches in GDAL / libgeotiff. People will appreciate looking at
>> frmts/gtiff/gt_wkt_srs.cpp
>>
>> I'm not sure if that GeoTIFF file is correct or not regarding the
>> values of the
>> false easting/northing with respect to the GeoTIFF spec (many in the
>> wild are
>> broken), but when looking at what happens, I noticed that a
>> transformation
>> from meters to feet happened twice. The attached patch fixes that
>> (and the
>> corner coordinates are around 75°), but I'm not sure if it is completely
>> correct in all cases.
>>
>> You should likely create a GDAL Trac ticket with your report and the
>> patch,
>> but I'm not sure who will dare hurting his head against this wall...
>>
>> Even
>>
>> >
>> > Thanks,
>> > Scott
>>
>> --
>> Geospatial professional services
>> http://even.rouault.free.fr/services.html
>>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org <mailto:gdal-dev at lists.osgeo.org>
>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>>
>>
>>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>

More information about the gdal-dev mailing list