Skip to content

Instantly share code, notes, and snippets.

@BishopGIS
Last active September 29, 2021 11:16
Show Gist options
  • Save BishopGIS/073f3ad79babc1dc7099e2a02a556101 to your computer and use it in GitHub Desktop.
Save BishopGIS/073f3ad79babc1dc7099e2a02a556101 to your computer and use it in GitHub Desktop.
from osgeo import ogr, osr
pt = ogr.CreateGeometryFromWkt('POINT(37.0 55.0)')
epsg7683 = osr.SpatialReference()
epsg7683.ImportFromEPSG(7683)
epsg7683.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
# epsg7683.SetTOWGS84(23.57,-140.95,-79.8,0,0.35,0.79,-0.22)
epsg4326 = osr.SpatialReference()
epsg4326.ImportFromEPSG(4326)
epsg4326.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
epsg28408 = osr.SpatialReference()
epsg28408.ImportFromEPSG(28408)
# epsg28408.SetTOWGS84(23.57,-140.95,-79.8,0,0.35,0.79,-0.22)
epsg28408.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
pt.AssignSpatialReference(epsg7683)
pt.ExportToWkt()
pt_new = pt.Clone()
pt_new.TransformTo(epsg28408)
pt_new.ExportToWkt()
pt_old = pt_new.Clone()
pt_old.TransformTo(epsg4326)
pt_old.ExportToWkt()
print('Delta X {}/Y {}'.format(abs(pt.GetX() - pt_old.GetX()), abs(pt.GetY() - pt_old.GetY())))
pt_old1 = pt_new.Clone()
pt_old1.TransformTo(epsg7683)
pt_old1.ExportToWkt()
pt_old2 = pt_old.Clone()
pt_old2.TransformTo(epsg28408)
pt_old2.ExportToWkt()
print('Delta X {}/Y {}'.format(abs(pt_new.GetX() - pt_old2.GetX()), abs(pt_new.GetY() - pt_old2.GetY())))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment