[gdal-dev] UTM -> LAEA

[Frank Warmerdam]

Thomas Becker wrote:
> Hello listers,
> 
> using gdalwarp to transform a Landsat image from UTM to lambert 
 > Azimuthal Equal Area, the output is showing "unknown" for the
 > Datum and "unnamed" for the Spheroid.
> 
> I was using the following command:
> gdalwarp -t_srs '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 
 >          +y_0=3210000 +ellps=GRS80' p193r002_7t19990711_z33_nn40.tif test.tif
...
> PROJCS["unnamed",
>     GEOGCS["GRS 1980(IUGG, 1980)",
>         DATUM["unknown",
>             SPHEROID["unnamed",6378137,298.2572221010042]],
>         PRIMEM["Greenwich",0],
>         UNIT["degree",0.0174532925199433]],
>     PROJECTION["Lambert_Azimuthal_Equal_Area"],
>     PARAMETER["latitude_of_center",52],
>     PARAMETER["longitude_of_center",10],
>     PARAMETER["false_easting",4321000],
>     PARAMETER["false_northing",3210000],
>     UNIT["metre",1,
...
> Can you tell me how i have to set the parameters for the DATUM and 
 > SPHEROID? Wenn i am using the EPSG:3035 ther will be no adjustments
 > for lat_0 and so on. Is there a possibility to use these parameters
 > as well with the EPSG Code?

Thomas,

GDAL translates:

  +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80

Into WKT as:

   PROJCS["unnamed",
     GEOGCS["GRS 1980(IUGG, 1980)",
         DATUM["unknown",
             SPHEROID["GRS80",6378137,298.257222101]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Lambert_Azimuthal_Equal_Area"],
     PARAMETER["latitude_of_center",52],
     PARAMETER["longitude_of_center",10],
     PARAMETER["false_easting",4321000],
     PARAMETER["false_northing",3210000]]

The name for +ellps=GRS80 comes from a hardcoded internal table
in the proj.4 string translation (adapted from PROJ.4 itself).
The spheroid name gets lost when the geotiff writer since there
is no code for capturing user provided spheroid names (assuming
there is a particular geokey for it).

EPSG:3035 translates as:

PROJCS["ETRS89 / ETRS-LAEA",
     GEOGCS["ETRS89",
         DATUM["European_Terrestrial_Reference_System_1989",
             SPHEROID["GRS 1980",6378137,298.257222101,
                 AUTHORITY["EPSG","7019"]],
             AUTHORITY["EPSG","6258"]],
         PRIMEM["Greenwich",0,
             AUTHORITY["EPSG","8901"]],
         UNIT["degree",0.01745329251994328,
             AUTHORITY["EPSG","9122"]],
         AUTHORITY["EPSG","4258"]],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]],
     PROJECTION["Lambert_Azimuthal_Equal_Area"],
     PARAMETER["latitude_of_center",52],
     PARAMETER["longitude_of_center",10],
     PARAMETER["false_easting",4321000],
     PARAMETER["false_northing",3210000],
     AUTHORITY["EPSG","3035"]]

I don't really understand your point about there being no
adjustements to +lat_0.  The latitude of center (+lat_0)
seems adequately addressed in this expansion.  But, you can actually
use WKT as input to gdalwarp.  So for instance, if you saved the
following as a file, you could use that filename as the argument
to -t_srs in gdalwarp.

PROJCS["ETRS89 / ETRS-LAEA",
     GEOGCS["ETRS89",
         DATUM["European_Terrestrial_Reference_System_1989",
             SPHEROID["GRS 1980",6378137,298.257222101,
                 AUTHORITY["EPSG","7019"]],
             AUTHORITY["EPSG","6258"]],
         PRIMEM["Greenwich",0,
             AUTHORITY["EPSG","8901"]],
         UNIT["degree",0.01745329251994328,
             AUTHORITY["EPSG","9122"]],
         AUTHORITY["EPSG","4258"]],
     PROJECTION["Lambert_Azimuthal_Equal_Area"],
     PARAMETER["latitude_of_center",52],
     PARAMETER["longitude_of_center",10],
     PARAMETER["false_easting",4321000],
     PARAMETER["false_northing",3210000],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]

This gives you the standard definition of ETRS89 while
still allowing you to modify the projection parameters if you
want.  It comes back from GeoTIFF looking like:

PROJCS["ETRS89 / ETRS-LAEA",
     GEOGCS["ETRS89",
         DATUM["European_Terrestrial_Reference_System_1989",
             SPHEROID["GRS 1980",6378137,298.2572221010002,
                 AUTHORITY["EPSG","7019"]],
             AUTHORITY["EPSG","6258"]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4258"]],
     PROJECTION["Lambert_Azimuthal_Equal_Area"],
     PARAMETER["latitude_of_center",52],
     PARAMETER["longitude_of_center",10],
     PARAMETER["false_easting",4321000],
     PARAMETER["false_northing",3210000],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]

I hope this helps though I fear I've missed part of your point.

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org