pythongeopandasgeographic-lib

project local coordinates to global GPS with reference point


I have a bunch of shapes (e.g. shapely LineStrings or Polygons) in a geopandas GeoDataFrame. The shapes specify coordinates in a local 200x200 meters grid, i.e. all coordinates are between (0, 0) and (200, 200).

I now would like to "place" these lines globally. For this, I want to specify a GPS Point (with a given lat/lon) as a reference.

My first (naive) approach would be to use geographiclib, take all shapes' coords (in local X/Y) and apply the following transformation and "recreate" the shape:

# Convert coordinates to GPS location
from shapely.geometry import LineString
from geographiclib.geodesic import Geodesic
geod = Geodesic.WGS84  # the base geodesic (i.e. the world)

origin = (48.853772345870176, 2.350983211585546)  # this is somewhere in Paris, for example

def local_to_latlong(x, y, orientation=0, scale=1):
    """ Two step process.
    - First walk x meters to east from origin.
    - Then, from that point, walk y meters north from origin.
    
    Optional: 
    - orientation allows to "spin" the coordinates
    - scale allows to grow/shrink the distances
    """
    go_X = geod.Direct(*origin, orientation + 90, x * scale)  # x is East-coordinate
    go_Y = geod.Direct(go_X["lat2"], go_X["lon2"], orientation + 0, y * scale)  # y is North-coordinate
    return go_Y["lat2"], go_Y["lon2"]


original_line = LineString([(0,0), (100,100), (200,100)])
global_line = LineString([local_to_latlong(x, y) for y, x in original_line.coords])

However, I hope that this is not the smartest way to do it, and that there are smarter ways out there...

I would like to apply such a transformation onto any shape within a GeoDataFrame. Ideally, it would work using a "to_crs", but I am not sure how to transform the shapes so they are "in reference to a origin" and which crs to use.


Solution

  • import shapely.geometry
    import geopandas as gpd
    import pandas as pd
    import numpy as np
    
    # generate some polygons (squares), where grid is 200*200
    gdf = gpd.GeoDataFrame(
        geometry=pd.DataFrame(
            np.repeat(np.sort(np.random.randint(0, 200, [20, 2]), axis=1), 2, axis=1)
        ).apply(lambda d: shapely.geometry.box(*d), axis=1)
    )
    # chage to linestrings, clearer when we plot
    gdf["geometry"] = gdf["geometry"].exterior
    
    origin = (2.350983211585546, 48.853772345870176)  # this is somewhere in Paris, for example
    
    # work out utm crs of point.  utm is in metres
    gdf_o = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(origin)], crs="EPSG:4326")
    crs = gdf_o.estimate_utm_crs()
    # where is origin in utm zone
    xo,yo = gdf_o.to_crs(crs).loc[0,"geometry"].xy
    
    # translate custom zone to co-ordinates of utm zone
    # assume point is center of 200x200 grid (hence subtract 100)
    gdf_gps = gdf["geometry"].translate(xoff=xo[0]-100, yoff=yo[0]-100).set_crs(crs).to_crs("epsg:4326")
    
    # plot on map to show it has worked...
    m = gdf_gps.explore()
    m = gdf_o.explore(m=m, color="red", marker_kwds={"radius":20})
    m
    
    
    

    enter image description here