[Gdal-dev] Re: [Proj] SRS ogr2ogr - cs2cs - postgis

[Jaime Mejía]

Paul Kelly wrote:

> On Sat, 5 Jun 2004, [ISO-8859-1] Jaime Mej?a wrote:
>
>> Hello.
>>
>> I've simplified my problem, in order to get some help:
>> I have a shape file (attached) containing a single point 
>> (994097.623824 1003980.4918) in DATUM BOGOTA, called punto_prueba.shp.
>>
>> When I try to transform it to WGS84 DATUM using ogr2ogr:
>> ----------------------
>> ogr2ogr -s_srs EPSG:21892 -t_srs EPSG:4326 
>> /home/jaime/resultado_wgs84 /home/jaime/punto_prueba.shp
>> ----------------------
>> The point now have the following coordinates: 
>> -74.1341093150,4.6350408295
>>
>> If I take the shape into postgis, and try to transform it, I get the 
>> following:
>>
>> ----------------------
>> SELECT num_radicacion,the_geom,transform(the_geom,1) AS wgs84 FROM 
>> expedientes WHERE num_radicacion=199931011;
>> num_radicacion |                 the_geom                 | wgs84
>> ----------------+------------------------------------------+-------------------------------------------------- 
>>
>>     199931011 | SRID=2;POINT(994097.623824 1003980.4918) | 
>> SRID=1;POINT(-74.1306989066892 4.63219422246886)
>> ----------------------
>>
>> Finally, using cs2cs:
>> ----------------------
>> jaime at fileserver[/maintask/clientes/creg]# cs2cs -v +proj=tmerc 
>> +lat_0=4.599047222 +lon_0=-74.080916667 +k=1.000000 +x_0=1000000.000 
>> +y_0=1tl +towgs84=307,304,-318,0,0,0,0 +to +proj=latlong +datum=WGS84
>> # ---- From Coordinate System ----
>> #Transverse Mercator
>> #       Cyl, Sph&Ell
>> # +proj=tmerc +lat_0=4.599047222 +lon_0=-74.080916667 +k=1.000000
>> # +x_0=1000000.000 +y_0=1000000.000 +ellps=intl
>> # +towgs84=307,304,-318,0,0,0,0
>> #--- following specified but NOT used
>> # +ellps=clrk66
>> # ---- To Coordinate System ----
>> #Lat/long (Geodetic)
>> #
>> # +proj=latlong +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>> 994097.623824 1003980.4918
>> 74d7'50.516"W   4d37'55.899"N 16.920
>> ----------------------
>> 74d7'50.516"W   4d37'55.899"N => -74.13069888888 4.63219416666
>>
>> As you can see, the results are different (+/- 500 m):
>> With ogr2ogr: -74.1341093150,4.6350408295
>> With postgis: -74.1306989066892 4.63219422246886
>> With cs2cs:-74.13069888888 4.63219416666
>>
>> So, I assume there is an error in the ogr2ogr command, or in ogr2ogr 
>> itself.
>
>
> No reason to assume an error---you gave ogr2ogr and cs2cs different 
> input.
>
> For ogr2ogr you said "EPSG:21892", which doesn't include any datum 
> transformation parameters, whereas for cs2cs you gave the parameters 
> specifically as +towgs84=307,304,-318
> It's not clear to me where the parameters are coming from in the 
> PostGIS example. The cs2cs example also looks like it's missing 
> ellipsoid information.
>
> If you want to reproduce the erroneous and correct result with cs2cs, 
> try:
>
> echo 994097.623824 1003980.4918 | cs2cs +init=epsg:21892 \
> +to +init=epsg:4326
>
> 74d8'2.794"W    4d38'5.669"N (wrong?)
>
> echo 994097.623824 1003980.4918 | cs2cs +init=epsg:21892 
> +towgs84=307,304,-318 \
> +to +init=epsg:4326
> 74d7'50.516"W   4d37'55.899"N (correct?)
>
> So something like -s_srs "+init=epsg:21892 +towgs84=307,304,-318" in 
> the ogr2ogr command line should probably work (if these parameters are 
> correct for the Bogota_1975 datum).
>
> Paul
> K

Paul,

We've tried with those parameters in the command line, but allways the 
command failed with error.

So we've tried with an external file:
-----------------------------------
ogr2ogr -s_srs merc_bogota -t_srs WGS84 resultado_wgs84 punto_prueba.shp
-----------------------------------
where the file merc_bogota contains the following line with the required
parameters to perform the transformation from BOGOTA DATUM to WGS84 DATUM:
+proj=tmerc +lat_0=4.599047222 +lon_0=-74.080916667 +k=1.000000
+x_0=1000000.000 +y_0=1000000.000 +ellps=intl +towgs84=307,304,-318,0,0,0,0
-----------------------------------
And this is the result:
-----------------------------------
jaime at centella[/pruebas]# ogrinfo -al resultado_wgs84/punto_prueba.shp
INFO: Open of `resultado_wgs84/punto_prueba.shp'
using driver `ESRI Shapefile' successful.

Layer name: punto_prueba
Geometry: Point
Feature Count: 1
Extent: (-74.134109, 4.635041) - (-74.134109, 4.635041)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
PROG_ID: Integer (5.0)
LABEL: String (50.0)
OGRFeature(punto_prueba):0
  PROG_ID (Integer) = 1
  LABEL (String) = Ramirez Elsa Graciela
  POINT (-74.13410932 4.63504083)

And the result (-74.13410932 4.63504083) is the same that with the command:
----------------------
ogr2ogr -s_srs EPSG:21892 -t_srs EPSG:4326 /home/jaime/resultado_wgs84 
/home/jaime/punto_prueba.shp
----------------------
The point now have the following coordinates: -74.1341093150,4.6350408295


In both cases, it seems that ogr ommits the 
+towgs84=307,304,-318,0,0,0,0 parameter.
What could I be doing wrong?

Thanks a lot.