Opened 11 years ago

Closed 6 years ago

#5098 closed enhancement (duplicate)

GDAL/OGR in Python - Using new ntv2 datum shift file with EPSG codes

Reported by: tcgeophysics Owned by: warmerdam
Priority: normal Milestone:
Component: default Version: 1.9.2
Severity: normal Keywords:
Cc:

Description

Hi,

I have been testing both pyproj and osgeo libraries for a script that only takes EPSG codes (source and target CRS) as inputs to reproject sets of points.

I successfully modified (I think) the pyproj source to add support for additional NTv2 Datum shift grids (proj +nadgrids parameter) instead of the 3 or 7 parameters normally used to convert the coordinates to WGS84 (proj +towgs84 parameter).

I modified the pj_datums.c and epsg files from the pyproj source and compiled.

How could I get the osgeo library in python to support the same modification?

I tested my modification to support AGD66 to GDA94 transformation in Australia. CRS details are below:

AGD66 => EPSG:4202 => +proj=longlat +ellps=aust_SA +towgs84=-117.808,-51.536,137.784,0.303,0.446,0.234,-0.29 +no_defs

=> +proj=longlat +ellps=aust_SA +nadgrids=A66_National_13.09.01.gsb +no_defs

GDA94 => EPSG:4283 => +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

To support a direct use of the NTv2(.gsb)file when initializing with an EPSG code calling for EPSG:4202, cs2cs now returns;

+proj=longlat +datum=AGD66 +no_defs

and the datum is identified in pj_datums.c as;

"AGD66", "nadgrids=A66_National_13.09.01.gsb","aust_SA","Australian Geodetic Datum 1966",

What could be recommend to get the same EPSG support with osgeo?

This could be extended to the following datums as they are country specific:

AGD 1966, AGD 1984, NAD83(HARN), Cape, Lisbon, Datum 73, DHDN, Corrego Alegre 1961, Corrego Alegre 1970-72, NZGD 1949 and OSGB36.

Change History (5)

comment:1 by Even Rouault, 11 years ago

What you can do is edit the GDAL ressource files located in the GDAL_DATA directory.

In gcs.csv, comment the line starting with 4202, with a # character for example In epsg.wkt, add the following new line :

4202,GEOGCS["AGD66",DATUM["Australian_Geodetic_Datum_1966",SPHEROID["Australian National Spheroid",6378160,298.25,AUTHORITY["EPSG","7003"]],EXTENSION["PROJ4_GRIDS","A66_National_13.09.01.gsb"],AUTHORITY["EPSG","6202"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4202"]]

Once this is done, you can test:

$ gdalsrsinfo  EPSG:4202

PROJ.4 : '+proj=longlat +ellps=aust_SA +nadgrids=A66_National_13.09.01.gsb +no_defs '

OGC WKT :
GEOGCS["AGD66",
    DATUM["Australian_Geodetic_Datum_1966",
        SPHEROID["Australian National Spheroid",6378160,298.25,
            AUTHORITY["EPSG","7003"]],
        EXTENSION["PROJ4_GRIDS","A66_National_13.09.01.gsb"],
        AUTHORITY["EPSG","6202"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4202"]]

comment:2 by tcgeophysics, 11 years ago

Your suggestion worked and I got the same result with gdalsrsinfo EPSG:4202.

In Python the following returned the appropriate proj string using a nadgrids parameter;

from osgeo import osr
srs = osr.SpatialReference()
srs.ImportFromEPSG(4202)
print srs.ExportToProj4()

This code returned:

+proj=longlat +ellps=aust_SA +nadgrids=A66_National_13.09.01.gsb +no_defs

How could gdal default to the towgs84 method if a point cannot be transformed using the gsb file specified with nadgrids?

This is could be useful in Europe where ED50 (EPSG:4230) was commonly used and only Spain (I think) provides a grid shift file, PENR2009.gsb, from ED50 to ETRS89 (EPSG:4258). If the above modification is done for EPSG:4230, conversion outside of Spain will fail as no towgs84 parameter is available.

in reply to:  2 comment:3 by Even Rouault, 11 years ago

How could gdal default to the towgs84 method if a point cannot be transformed using the gsb file specified with nadgrids?

This is essentially a PROJ.4 issue. GDAL could pass both the +towgs84 and +nadgrids and let proj4 decide which one to use. But the main changes would be in PROJ.4

comment:4 by tcgeophysics, 11 years ago

Modifying GDAL's epsg.wkt and gcs.csv files works great in the python terminal and by calling the python script through windows batch files.

But it seems that none of these changes are taken into account when creating an executable of my python script using py2exe. The transformations are reverted to their original towgs84 method...

Are the information contained in epsg.wkt and gcs.csv included in some compiled files that would be used by py2exe when creating the executable?

Any suggestion about why I am experiencing that issue is welcome!

Thanks

comment:5 by Even Rouault, 6 years ago

Resolution: duplicate
Status: newclosed

Appears to bea duplicate of #5085

Note: See TracTickets for help on using tickets.