[Gdal-dev] projection with ogr2ogr

[Mateusz Loskot]

Lijian Shi wrote:
> 
> Except WGS84 with epsg 4326, I also want to add the projection information:
> 
> PROJCS[" Projection Name = Lambert Azimuthal Equal-ar",
> ......
> PROJECTION["Lambert_Azimuthal_Equal_Area"],
> PARAMETER["latitude_of_center",45],
> PARAMETER["longitude_of_center",-100],
> PARAMETER["false_easting",0],
> PARAMETER["false_northing",0],
> UNIT["metre",1,
> AUTHORITY["EPSG","9001"]]]
> 
> So I tried the command: ogr2ogr -a_srs EPSG:4326 -t_srs "+proj=laea" out.shp
> point.shp
> 
> And then I get the "header" of shp file as followed:
> PROJCS["Lambert_Azimuthal_Equal_Area",
>     GEOGCS["GCS_WGS_1984",
>         DATUM["unknown",
>             SPHEROID["WGS84",6378137,298.257223563]],
>         PRIMEM["Greenwich",0],
>         UNIT["Degree",0.017453292519943295]],
>     PROJECTION["Lambert_Azimuthal_Equal_Area"],
>     PARAMETER["latitude_of_center",0],
>     PARAMETER["longitude_of_center",0],
>     PARAMETER["false_easting",0],
>     PARAMETER["false_northing",0],
>     UNIT["Meter",1]]
> 
> But there are still some difference between the "header" of shp file and
> geotiff file. I don't know whether my result is right. So look forward to
> your suggestion.

Lijian,

I see I've misunderstood you previously.
Let's summarize the problem:

1. You have a projected GeoTIFF
2. You have also a Shapefile that fits spatially the GeoTIFF,
   but without any projection assigned
3. You want to use OGR to assign the Shapefile with the same
   projection as in GeoTIFF.

Here is the procedure:

1. You can use gdalinfo to dump PROJECTION details in WKT format.
Actually, you have already done it.

2. Copy and paste the WKT string dumped from GeoTIFF into a text file,
for example mytiff.prj. Note, you don't need to format it in
a single line.

PROJCS[" Projection Name = Lambert Azimuthal Equal-ar",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",45],
    PARAMETER["longitude_of_center",-100],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

3. Use ogr2ogr to "translate" your unreferenced Shapefile with
assigning the projection from WKT you've already saved in mytiff.prj file:


ogr2ogr -a_srs mytiff.prj out.shp yourinputfile.shp

where yourinputfile.shp is the unreferenced Shapefile.

After these steps, you will get following files:
out.dbf
out.shp
out.shx
out.prj - here is the info about newly assigned projection


I hope it helps.

Cheers
-- 
Mateusz Loskot
http://mateusz.loskot.net