CODEWITHSUNDEEP
pythongridcoordinatesgispyproj
JaneDOEE
JaneDOEE

Reputation: 11

How to transform lat/lon coordinates from degree to km to create a 1 km grid

I have a dataframe containing lat/lon coordinates in decimal degree.

My goal is to aggregate the data on a rectangular grid of 1 km². For that matter, I transformed my coordinates into km based on the method described in Convert latitude, longitude to distance from equator in kilometers and round to nearest kilometer

The method consists in calculating the distance from a reference point to the points (lat=0, lon) and (lat, lon=0).

But it doesn't work, because it seems to depend on the reference point.

By taking my reference point as (lon_ref=mean(lon), lat_ref=mean(lat)), I end up aggregating in the same tile points that are 120km away from each other.

This is the code that I am using :

# get the coordinates of my reference point

lat_ref, lon_ref = data["lat"].mean() , data["lon"].mean()

# the distance function

from pyproj import Geod 
wgs84_geod = Geod(ellps='WGS84')

format = lambda x: wgs84_geod.inv(lon_ref,lat_ref,0,x)[2]/1000 #km 

format = lambda x: wgs84_geod.inv(lon_ref,lat_ref,x,0)[2]/1000 #km 

# Apply the function on my dataframe

data["londist"]=data['lon'].map(format)

data["latdist"]=data['lat'].map(format)

# round to the nearest km

step=1 # 1km

to_bin = lambda x: np.round(x / step) * step

data["latbin"] = data['latdist'].map(to_bin)

data["lonbin"] = data['londist'].map(to_bin)

This works for some lat/lon but not for others,

Example:

point1 (46.9574,4.29949) # lat,lon in °

point2( 46.9972 ,3.18153)

Calculate the distance and round using the above code:

point1 (latbin = 259 , lonbin=5205)

point2(latbin = 259 , lonbin=5205)

The two points will be aggregated together However, the distance between the two points is 85 km! 

dist=wgs84_geod.inv(4.29949,46.9574,3.18153,46.9972)[2]/1000

How can I solve this problem? Is there any other efficient method to make the aggregation given that I have 10 millions lat/lon in my dataframe?

Upvotes: 0

Views: 2228

Answers (1)

Michael Entin
Michael Entin

Reputation: 7764

It looks like you assign to format variable twice, second assignment erasing first one. Did you mean to have something like format_lon and format_lat?

Note that once you fix this, the result will still depend on reference point - this is inevitable when you project spherical coordinates to flat map. But the result should be reasonable.

Upvotes: 0

Related Questions