ali kanbar
ali kanbar

Reputation: 1

Datum transformation using python

I'm trying to do a levant to wgs84 transformation using python and Helmert's rule but i'm always getting a wrong outcome,can someone please tell me the reason? I have a levant shapefile and the expected wgs84 shapefile,i'm trying yo transform the levant one and compare it to the wgs84 one, but it's always shifted and i can't figure out why.

import geopandas as gpd
from pyproj import CRS, Transformer
import pyproj
from shapely.ops import transform
from functools import partial
from shapely.geometry import Polygon, Point

bound_crs_wkt = """
   PROJCS["Stereographic_Leb_Clarke1880",
   GEOGCS["GCS_Clarke_1880",
   DATUM["D_Clarke_1880",
   SPHEROID["Clarke_1880",6378249.14480801,293.466307655625]],
   PRIMEM["Greenwich",0.0],
   UNIT["Degree",0.0174532925199433]],
   PROJECTION["Double_Stereographic"],
   PARAMETER["False_Easting",0.0],
   PARAMETER["False_Northing",0.0],
   PARAMETER["Central_Meridian",39.15],
   PARAMETER["Scale_Factor",0.9995341],
   PARAMETER["Latitude_Of_Origin",34.2],
   UNIT["Meter",1.0]]

"""

shapefile_path = "Tyre_withTileMerging/Tyre_withTileMering_Stereo.shp"
stereo_levant_data = gpd.read_file(shapefile_path)


stereo_levant_data.set_crs(CRS.from_wkt(bound_crs_wkt), inplace=True, allow_override=True)

source_crs = CRS.from_wkt(bound_crs_wkt)
target_crs = CRS.from_epsg(4326)  # WGS84

tx = .....
ty = .....
tz = .....
s = .....
rx = .....
ry = .....
rz = .....

helmert_crs = CRS.from_string(f"""
    +proj=latlong +datum=WGS84 
    +towgs84={tx},{ty},{tz},{rx},{ry},{rz},{s}
""")

transformer = Transformer.from_crs(source_crs, helmert_crs, always_xy=True)


def apply_transformation(geom):
    if geom.is_empty:
        return geom
    return transform(partial(transformer.transform), geom)

stereo_levant_data['geometry'] = stereo_levant_data['geometry'].apply(apply_transformation)

stereo_levant_data_wgs84 = stereo_levant_data.to_crs(epsg = 32636)

output_path_wgs84 = "transformed_file.shp"
stereo_levant_data_wgs84.to_file(output_path_wgs84,encoding='utf-8')

print(f"Data has been transformed to WGS84 and saved to {output_path_wgs84}")

I'm expecting a correct transformation but i'm not able to do it

Upvotes: 0

Views: 65

Answers (0)

Related Questions