[Gdal-dev] projection with ogr2ogr

[Lijian Shi]

Hi, Mateusz,

Thank you for your help. Now I have resolved the problem according your
suggestion. But I met a new problem: after reprojecting the shp file with
the command: ogr2ogr -t_srs WGS84 -s_srs NAD27 out.shp point.shp,
I overlay the out.shp and point.shp in ArcGIS. But there are some shift
between tow files. Is there any problem with my command? 

Best regards,

Lijian

-----Original Message-----
From: Mateusz Loskot [mailto:mateusz at loskot.net] 
Sent: Wednesday, March 07, 2007 1:30 AM
To: Lijian Shi
Cc: gdal-dev at lists.maptools.org
Subject: Re: [Gdal-dev] projection with ogr2ogr

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