[gdal-dev] Confused about datum transformations

Simon Lyngby Kokkendorff silyko at gmail.com
Wed Oct 19 05:54:53 EDT 2011


Hi Marc,

  Just to understand your problem fully:  You have data in your database in
WGS84, in geographical coordinates, or what? That would correspond to a
srs=osr.SpatialReference(), srs.SetWellKnownGeogCS("WGS84").

   What kind of projection or geographical coordinate system do you want to
project to? It seems to be: EPSG:3557 which is: NAD83(NSRS2007) / Maine
East. Is this correct? If you have the GDAL_DATA env-variable pointing to
the right directory (containing EPSG definitions etc), gdal should be able
to also do:
srs2=osr.SpatialReference(), srs2.ImportFromEPSG(3557), which becomes:

>>> srs2.ExportToPrettyWkt()
'PROJCS["NAD83(NSRS2007) / Maine East",\n
GEOGCS["NAD83(NSRS2007)",\n
DATUM["NAD83_National_Spatial_Reference_System_2007",\n
SPHEROID["GRS 1980",6378137,298.257222101,\n
AUTHORITY["EPSG","7019"]],\n
       TOWGS84[0,0,0,0,0,0,0],\n
AUTHORITY["EPSG","6759"]],\n        PRIMEM["Greenwich",0,\n
AUTHORITY["EPSG","8901"]],\n
UNIT["degree",0.01745329251994328,\n
AUTHORITY["EPSG","9122"]],\n        AUTHORITY["EPSG","4759"]],\n
UNIT["metre",1,\n        AUTHORITY["EPSG","9001"]],\n
PROJECTION["Transverse_Mercator"],\n
PARAMETER["latitude_of_origin",43.66666666666666],\n
PARAMETER["central_meridian",-68.5],\n
PARAMETER["scale_factor",0.9999],\n
PARAMETER["false_easting",300000],\n    PARAMETER["false_northing",0],\n
AUTHORITY["EPSG","3557"],\n    AXIS["X",EAST],\n    AXIS["Y",NORTH]]'
'

If you want to perform the actual coordinate transformation, besides just
having descriptions of the systems, you'll definitely need to have proj4.dll
somewhere in your system, where gdal can find it (I guess besides gdal.dll
or somewhere in your PATH).

Cheers,
Simon


On Wed, Oct 19, 2011 at 11:21 AM, <dev4cx4m9z at snkmail.com> wrote:

> I have a big clue - I have inspected the output of ExportToWKT, and the
> results are interesting
>
>        With WGS84SRS
>            GDALError = .SetProjCS("WGS84 Datum")
>            Debug.Assert(GDALError = OGRError.OGRERR_NONE,
> GetLastErrorMsg())
>            GDALError = .SetWellKnownGeogCS("WGS84")          ' Define WGS84
> coordinate system
>            Debug.Assert(GDALError = OGRError.OGRERR_NONE,
> GetLastErrorMsg())
>        End With
>
> Produces PROJCS["WGS84 Datum",GEOGCS["WGS
> 84",DATUM["WGS_1984",SPHEROID["WGS
> 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]]
>        With UTMSRS
>            GDALError = .SetProjCS("UTM coordinate system")
>            Debug.Assert(GDALError = OGRError.OGRERR_NONE,
> GetLastErrorMsg())
>            GDALError = .SetWellKnownGeogCS("WGS84")              ' Define a
> UTM coordinate system
>             Debug.Assert(GDALError = OGRError.OGRERR_NONE,
> GetLastErrorMsg())
>        End With
>
> Produces PROJCS["UTM coordinate system",GEOGCS["WGS
> 84",DATUM["WGS_1984",SPHEROID["WGS
> 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]]
>
>       With MapSRS
>            GDALError = .SetProjCS("Map coordinate system")
>            Debug.Assert(GDALError = OGRError.OGRERR_NONE,
> GetLastErrorMsg())
>             GDALError = .ImportFromEPSG(i)
>             Debug.Assert(GDALError = OGRError.OGRERR_NONE,
> GetLastErrorMsg())
>        End With
>
> Produces
> GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS
> 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]]
>
> The error when I attempt to create the transformation is always with the
> SRS that seems to be a projection definition, not a datum. ImportFromEPSG
> seems to be correctly producing a GEOGCS, but SetWellKnownGeogCS seems to be
> producing a "PROJCS"
>
> I feel I'm very close to the answer, but I still can't see it.
>
> -----Original Message-----
> From: Peter J Halls P.Halls-at-york.ac.uk |GDAL| [mailto:
> 6wyv643jct at sneakemail.com]
> Sent: Wednesday, October 19, 2011 7:54 PM
> To: marc.hillman at tpg.com.au; dev4cx4m9z at snkmail.com
> Subject: Re: [gdal-dev] Confused about datum transformations
>
> A comment from another GDAL user, rather than the development team.
>
> GDAL and OGR are actually C++ tools that have C, Java, Python and other
> language bindings.  What you have picked up is an artefact of the fact that
> GDAL is not just C.
>
> The proj4 library for your platform is automatically loaded with GDAL and
> OGR:
> it must be present but no explicit action is required by the user.
>
> Chaitanya's advice was to double check that the transformation parameters
> are being passed correctly.  Datum transformations *are* fiddly (full stop)
> to specify and this advice is sound.  You appear to be specifying two parts
> of the transformation, presumably being content with the default values for
> the other parameters.  Ideally, you need to get the system to regurgitate
> its understanding of your parameters in the WKT version of the final
> transformation specification - ie the complete string.  This is what you are
> being advised to do.
>
> Problems like this are not confined to GDAL / OGR: I've lost count of the
> number of times I've handled projection problems here.  What sounds to be a
> simple operation is not trivial and unpicking problems most often requires
> an understanding of the transformation method and all the parameters
> applied.
>
> Hoping that helps,
>
> Peter
>
> dev4cx4m9z at snkmail.com wrote:
> > I have not used org2ogr. I have just realised that the C API for GDAL
> > uses a mixture of STD and CDECL. How annoying. Anyway – I have
> > carefully checked and all are correct. I get a stack balance problem
> > if they are not anyway. It’s hard to see how I could get the arguments
> wrong – one is the string “WGS84”, and the other is an integer.
> >
> >
> >
> > I still don’t know whether proj.dll is required.
> >
> >
> >
> > From: Chaitanya kumar CH chaitanya.ch-at-gmail.com |GDAL|
> > [mailto:5rm5yuc82t at sneakemail.com]
> > Sent: Wednesday, October 19, 2011 6:36 PM
> > To: marc.hillman at tpg.com.au; dev4cx4m9z at snkmail.com
> > Cc: gdal-dev at lists.osgeo.org
> > Subject: Re: [gdal-dev] Confused about datum transformations
> >
> >
> >
> > Marc,
> >
> > You don't need to rebuild proj.4; you are already using it.. Since you
> were able to use ogr2ogr to do the transformation, it works.
> >
> > What you need to do is find out if the arguments you are passing to the
> transformation builder are good. Especially 'MapSRS'. If it doesn't show any
> problem, plug some debug messages in the OSR code in ogrct.cpp and rebuild
> gdal.
> >
> > On Wed, Oct 19, 2011 at 12:46 PM, <dev4cx4m9z at snkmail.com> wrote:
> >
> > Thanks for all the suggestions – they have given me many possibilities to
> eliminate, but I still haven’t really got an answer to the original
> question.
> >
> >
> >
> > Am I required to have proj.dll built with support for the
> datums/projections I want to use, or is it sufficient to set GDAL_DATA such
> that gcs.csv and pcs.csv are visible? I’m working in VB.NET, and I’m not
> confident I can rebuild PROJ with the required datums, so I wish to avoid if
> possible. I’d even be prepared to use SetGeogCS explicitly as I only use a
> few datums. If I use SetGeogCS, I can’t see what purpose proj.dll serves.
> >
> >
> >
> > I have set GDAL_DATA to the correct folder.
> >
> > I have checked that the datum numbers I use are in the gcs.csv file.
> >
> >
> >
> > From: gdal-dev-bounces at lists.osgeo.org
> > [mailto:gdal-dev-bounces at lists.osgeo.org
> > <mailto:gdal-dev-bounces at lists.osgeo..org> ] On Behalf Of
> > dev4cx4m9z at snkmail.com
> > Sent: Tuesday, October 18, 2011 4:12 AM
> > To: gdal-dev at lists.osgeo.org
> > Subject: RE: [gdal-dev] Confused about datum transformations
> >
> >
> >
> > I’m not sure the code will help too much, but here are the relevant
> > code fragments. I’m using a hand crafted VB.net wrapper which explains
> > some of the syntax
> >
> >
> >
> > Private WGS84SRS As New OGRSpatialReference     ' Coordinate system in
> WGS84
> >
> > Private MapSRS As OGRSpatialReference       ' Coordinate system of map
> >
> > Private ctWGS84toMap As OGRCoordinateTransformation       ' Coordinate
> transform from WGS84 to map
> >
> >
> >
> > With WGS84SRS
> >
> >     GDALError = .SetProjCS("WGS84 Datum")
> >
> >     Debug.Assert(GDALError = OGRError.OGRERR_NONE, GetLastErrorMsg())
> >
> >     GDALError = .SetWellKnownGeogCS("WGS84")          ' Define WGS84
> coordinate system
> >
> >     Debug.Assert(GDALError = OGRError.OGRERR_NONE, GetLastErrorMsg())
> >
> > End With
> >
> >
> >
> > With MapSRS
> >
> >      GDALError = .SetProjCS("Map coordinate system")
> >
> >      Debug.Assert(GDALError = OGRError.OGRERR_NONE, GetLastErrorMsg())
> >
> >      GDALError = .ImportFromEPSG(3577)
> >
> >      Debug.Assert(GDALError = OGRError.OGRERR_NONE, GetLastErrorMsg())
> >
> > End With
> >
> >
> >
> > ctWGS84toMap = GDAL.CreateCoordinateTransformation(WGS84SRS, MapSRS)
>    ' Create a coordinate transformation - WGS84 to map
> >
> > Debug.Assert(ctWGS84toMap IsNot Nothing, GetLastErrorMsg())
> >
> >
> >
> > Marc,
> >
> > Can you provide the code leading to the error?
> >
> > On Tue, Oct 18, 2011 at 1:21 PM, Marc Hillman <dev4cx4m9z at snkmail.com>
> wrote:
> >
> > I have a GIS database, and all data is stored in WGS84 datum. I have
> developed an application using GDAL, and I now need to translate my WGS84
> coordinates to several different datums, depending what map is loaded. I
> seem to be going around in circles, and need a nudge in the right direction.
> I am defining an OGRSpatialReference for each datum (WGS84 and current map)
> and then attempting to create an OGRCoordinateTransformation for the
> transformation.. This last step fails with a PROJ.4 error (“No PROJ.4
> translation for destination SRS, coordinate transformation initialization
> failed”
> >
> >
> >
> > My question is - Do I need PROJ.4 to do the translation? Can I not just
> use ImportFromWKT to define my OGRSpatialReference? If so, is there a nice
> canned list of all the WKT somewhere?  I’m keen to avoid PROJ.4 if possible
> as I’m not sure I can build it.
> >
> >
> >
> >
> > _______________________________________________
> > gdal-dev mailing list
> > gdal-dev at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
> >
> >
> >
> >
> >
> >
> > ----------------------------------------------------------------------
> > --
> >
> > _______________________________________________
> > gdal-dev mailing list
> > gdal-dev at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> --
>
> --------------------------------------------------------------------------------
> Peter J Halls, GIS Advisor & Team Leader Applications Support and Training,
>                Information Directorate, University of York
> Telephone: 01904 323806     Fax: 01904 323740
> Snail mail: Harry Fairhurst Building, University of York,
>             Heslington, York YO10 5DD
> This message has the status of a private and personal communication
>
> --------------------------------------------------------------------------------
>
>
>
> -----
> No virus found in this message.
> Checked by AVG - www.avg.com
> Version: 2012.0.1831 / Virus Database: 2092/4560 - Release Date: 10/18/11
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20111019/d6a898e3/attachment-0001.html


More information about the gdal-dev mailing list