[gdal-dev] Reprojecting with Same Projection but Different Central_Meridians

[Tom Hayden]

You guys are the best - thank you. I thought I was going crazy. Brad - your
code worked great on version 3.6.0 but when I was running this locally on
3.0.4. I get the wrong answer at the end. No problem, though, because this
is fixed now. Thank you!!

>>> transform_func = osr.CoordinateTransformation(srs2, srs)
>>> u = ogr.Geometry(ogr.wkbPoint)
>>> u.AddPoint(0,0)
>>> u.Transform(transform_func)
0
>>> print(u.ExportToWkt())
POINT (0 0 0)
>>> import osgeo
>>> osgeo.__version__
'3.0.4'


Meanwhile, works great in 3.6.0

>>> import osgeo.osr as osr
>>> import osgeo.ogr as ogr
>>> p =
'PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-41.34804059746178],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]'
>>> p2 =
"""PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]"""
>>> srs2 = osr.SpatialReference(wkt=p2)
>>> srs = osr.SpatialReference(wkt=p)
>>> transform_func = osr.CoordinateTransformation(srs2, srs)
>>> u = ogr.Geometry(ogr.wkbPoint)
>>> u.AddPoint(0,0)
>>> u.Transform(transform_func)
0
>>> print(u.ExportToWkt())
POINT (4602842.82460905 0.0 0)
>>> osgeo.__version__
'3.6.0'




On Thu, Mar 16, 2023 at 3:14 AM Brad Hards <bradh at frogmouth.net> wrote:

> On Thursday, 16 March 2023 3:13:59 PM AEDT Tom Hayden wrote:
> > >>> u.Transfrom(transform_func)
> Did this line not give you an error?
>
> In any case, I used:
>
>
> import osgeo.osr as osr
> import osgeo.ogr as ogr
> p =
> 'PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-41.34804059746178],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]'
> p2 =
> """PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]"""
> srs2 = osr.SpatialReference(wkt=p2)
> srs = osr.SpatialReference(wkt=p)
> transform_func = osr.CoordinateTransformation(srs2, srs)
> u = ogr.Geometry(ogr.wkbPoint)
> u.AddPoint(0,0)
> u.Transform(transform_func)
> print(u.ExportToWkt())
>
> which resulted in:
> POINT (4602842.82460905 0.0 0)
>
>
> Brad
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230316/764a497d/attachment.htm>