0

I would like to use pint to convert degrees (distance in a geographic CRS) into nautical miles.

https://geopandas.org/docs/reference/api/geopandas.GeoDataFrame.sjoin_nearest.html outputs distance in degree for epsg:4326.

Given distance (in nm) varies from equator to pole i'm not sure if this is possible.

I could use a rule of thumb of 1 deg ~= 111 km ~= 60 nm.

Perhaps it can be calculated using the starting point and distance using something like: https://github.com/anitagraser/movingpandas/blob/master/movingpandas/geometry_utils.py#L38

This code is also useful: https://geopy.readthedocs.io/en/stable/#module-geopy.distance

Here's some code to test:

import pandas as pd
import geopandas as gpd

df = pd.DataFrame({"lon": [0], "lat": [0]})
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="epsg:4326")

df2 = pd.DataFrame({"lon": [1, 2], "lat": [0, 0]})
gdf_pts = gpd.GeoDataFrame(df2, geometry=gpd.points_from_xy(df2["lon"], df2["lat"]), crs="epsg:4326")

value = gdf_pt.sjoin_nearest(gdf_pts, distance_col="distances")["distances"].values[0]


import pint

l = value * ureg.arcdegree
Ray Bell
  • 1,508
  • 4
  • 18
  • 45

2 Answers2

1

Probably best to throw it to Mercator and use that if you can

import pint_pandas

gdf = gdf_pt.to_crs("EPSG:3395").sjoin_nearest(gdf_pts.to_crs("EPSG:3395"), distance_col="distances")
gdf["distance"] = gdf["distance"].astype("pint[meter]").pint.to("nautical_mile")
Ray Bell
  • 1,508
  • 4
  • 18
  • 45
0

This function, which I've pulled from existing code, computes the distance in meters between two lat/long sets. "rlat" and "rlong" are expressed in radians; you'll have to do the conversion from degrees. To get nm instead of meters, just set R to 3440.

from math import *

# Radius of the earth, in meters.

R = 6371000

# Return distance between two lat/longs.

def distance( pt1, pt2 ):
    rlat1 = pt1.rlat
    rlat2 = pt2.rlat
    dlat = pt2.rlat - pt1.rlat
    dlong = pt2.rlong - pt1.rlong

    a = sin(dlat/2) * sin(dlat/2) + cos(rlat1) * cos(rlat2) * sin(dlong/2) * sin(dlong/2)
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    return R * c
Tim Roberts
  • 48,973
  • 4
  • 21
  • 30