test_cs.py

from nomk import coord

from osgeo import ogr, osr

epsg7683 = osr.SpatialReference()
epsg7683.ImportFromEPSG(4284)
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)
# epsg4326.SetTOWGS84(23.57,-140.95,-79.8,0,0.35,0.79,-0.22)

epsg28408 = osr.SpatialReference()
epsg28408.ImportFromEPSG(28408)
epsg28408.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
# epsg28408.SetTOWGS84(23.57,-140.95,-79.8,0,0.35,0.79,-0.22)

epsg3857 = osr.SpatialReference()
epsg3857.ImportFromEPSG(3857)
epsg3857.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
# epsg3857.SetTOWGS84(23.57,-140.95,-79.8,0,0.35,0.79,-0.22)

x = 43.5 + (44.0 - 43.5) / 2
y = 54.333333 + (54.666667 - 54.333333) / 2

# letter, number, last_letter_num, is_south = nomk_parser.parse100k('N-38-052')
# _ , min_x, max_x, min_y, max_y = nomk_ptext.text_to_100k(letter, number, last_letter_num, is_south)

_ , min_x, max_x, min_y, max_y = coord.coords_to_100k(x, y)

steps = 30

ring = ogr.Geometry(type=ogr.wkbLinearRing)
step_size_x = abs((max_x - min_x) / steps)
step_size_y = abs((max_y - min_y) / steps)

for y in range(steps + 1):
    ring.AddPoint(min_x, min_y + step_size_y * y)

for x in range(steps + 1):
    ring.AddPoint(min_x + step_size_x * x, max_y)

for y in range(steps + 1):
    ring.AddPoint(max_x, max_y - step_size_y * y)

for x in range(steps + 1):
    ring.AddPoint(max_x - step_size_x * x, min_y)

polygon = ogr.Geometry(type=ogr.wkbPolygon)
polygon.AddGeometryDirectly(ring)
# polygon.AssignSpatialReference(epsg7683)
out_geom = ogr.Geometry(type=ogr.wkbMultiPolygon)
out_geom.AddGeometryDirectly(polygon)
out_geom.AssignSpatialReference(epsg7683)

out_geom_1_1 = out_geom.Clone()
out_geom_1_1.TransformTo(epsg28408)

out_geom_1_2 = out_geom_1_1.Clone()
out_geom_1_2.TransformTo(epsg3857)
# out_geom_1_2.ExportToWkt()

out_geom_2_1 = out_geom.Clone()
out_geom_2_1.TransformTo(epsg3857)
# out_geom_2_1.ExportToWkt()

ring_1_2 = out_geom_1_2.GetGeometryRef(0).GetGeometryRef(0)
ring_2_1 = out_geom_2_1.GetGeometryRef(0).GetGeometryRef(0)

points_1_2 = ring_1_2.GetPoints()
points_2_1 = ring_2_1.GetPoints()

for i in range(steps):
    delta_x = abs(points_1_2[i][0] - points_2_1[i][0])
    delta_y = abs(points_1_2[i][1] - points_2_1[i][1])

    print('Delta {}/{}'.format(delta_x, delta_y))

Delta 2.79396772385e-09/0.0
Delta 2.79396772385e-09/0.0
Delta 9.31322574615e-10/0.0
Delta 0.0/5.58793544769e-09
Delta 9.31322574615e-10/0.0
Delta 1.86264514923e-09/3.72529029846e-09
Delta 0.0/0.0
Delta 1.86264514923e-09/3.72529029846e-09
Delta 0.0/0.0
Delta 1.86264514923e-09/0.0
Delta 0.0/2.79396772385e-09
Delta 1.86264514923e-09/0.0
Delta 1.86264514923e-09/3.72529029846e-09
Delta 0.0/0.0
Delta 9.31322574615e-10/0.0
Delta 9.31322574615e-10/0.0
Delta 1.86264514923e-09/0.0
Delta 2.79396772385e-09/0.0
Delta 9.31322574615e-10/0.0
Delta 9.31322574615e-10/0.0
Delta 1.86264514923e-09/0.0
Delta 0.0/0.0
Delta 1.86264514923e-09/3.72529029846e-09
Delta 0.0/0.0
Delta 1.86264514923e-09/0.0
Delta 0.0/0.0
Delta 0.0/0.0
Delta 1.86264514923e-09/0.0
Delta 9.31322574615e-10/5.58793544769e-09
Delta 1.86264514923e-09/0.0