2

I have 9888562 records in dataframe and I would like to convert my lat, long to UTM x,y. according to my code, I have used pyproj package but because my data are too much it takes a long time and finally, it doesn't work. I wonder whether you know another way or package that I can use for my data?

def rule(row):
    p = Proj(proj='utm',zone=10,ellps='WGS84', preserve_units=False)
    x,y = p(row["LON"], row["LAT"])
    return pd.Series({"X": x , "Y": y})
My_data = My_data.merge(My_data.apply(rule, axis=1), left_index= True, right_index= True)

My data is like this

  • what do you mean it doesn't work? can you give a pair of input points and output returned? – Ian Turton Sep 8 '19 at 9:52
  • Yes. I have tried my code for several rows and it works but when I try my code for large data(9888562 rows) it takes along time. I waited for 8 hours but it was running and I interrupted it. – Danial Sep 8 '19 at 10:05
3

You can directly use shapely or GeoPandas but with 9888562 records It will take a long time to do (if you want a Progress bar during the pandas operations, you can use tqdm:

enter image description here)

1) With your solution and the first 4 points

import pandas as pd
df = pd.DataFrame({'LAT':[47.9767,47.9803,47.9801,47.9798], 'LON':[-122.2450,-122.2458,-122.2472,-122.2465]})
My_data = df.merge(df.apply(rule, axis=1), left_index= True, right_index= True)
print(My_data)
     LAT       LON              X             Y
0  47.9767 -122.2450  556344.793287  5.313986e+06
1  47.9803 -122.2458  556281.177583  5.314386e+06
2  47.9801 -122.2472  556176.922041  5.314363e+06
3  47.9798 -122.2465  556229.484294  5.314330e+06

2) With shapely

from shapely.geometry import Point
df['geometry'] = df.apply(lambda row: Point(row.LON, row.LAT), axis=1)
print(df)
    LAT       LON                 geometry
0  47.9767 -122.2450  POINT (-122.245 47.9767)
1  47.9803 -122.2458  POINT (-122.2458 47.9803)
2  47.9801 -122.2472  POINT (-122.2472 47.9801)
3  47.9798 -122.2465  POINT (-122.2465 47.9798)
# change the projection of points
from functools import partial
from shapely.ops import transform
import pyproj   
project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'), # source coordinate system
    pyproj.Proj('+proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')) # destination coordinate system
sdf = df.copy()
sdf['project'] = sdf.geometry.apply(lambda row: transform(project,row))
print(sdf)
    LAT       LON                   geometry                  project
0  47.9767 -122.2450  POINT (-122.245 47.9767)   POINT (556344.7932874847 5313986.318933257)
1  47.9803 -122.2458  POINT (-122.2458 47.9803)  POINT (556281.1775833125 5314385.854768319)
2  47.9801 -122.2472  POINT (-122.2472 47.9801)  POINT (556176.9220411006 5314362.605105238)
3  47.9798 -122.2465  POINT (-122.2465 47.9798)  POINT (556229.4842943319 5314329.77191091)

3) Directly with GeoPandas

import geopandas as gpd
# change to GeoDataFrame
geodf = gpd.GeoDataFrame(df)
geodf.crs= {'init': 'epsg:4326'}
# change the projection of geodf
geodf = geodf.to_crs("+proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_def")
print(geodf)
       LAT       LON                          geometry
0  47.9767 -122.2450  POINT (556344.7932874847 5313986.318933257)
1  47.9803 -122.2458  POINT (556281.1775833125 5314385.854768319)
2  47.9801 -122.2472  POINT (556176.9220411006 5314362.605105238)
3  47.9798 -122.2465  POINT (556229.4842943319 5314329.77191091)

4) If you want the X and Y values, it is easy

geodf["x"] = geodf.geometry.apply(lambda row:row.x)
geodf["y"] = geodf.geometry.apply(lambda row:row.y)
print(geodf)
     LAT     LON                      geometry                       x                    y
0  47.9767 -122.2450  POINT (556344.7932874847 5313986.318933257)  556344.793287  5.313986e+06
1  47.9803 -122.2458  POINT (556281.1775833125 5314385.854768319)  556281.177583  5.314386e+06
2  47.9801 -122.2472  POINT (556176.9220411006 5314362.605105238)  556176.922041  5.314363e+06
3  47.9798 -122.2465  POINT (556229.4842943319 5314329.77191091)   556229.484294  5.314330e+06
1

UPDATE:

After thinking about it, the most efficient method for you to transform the coordinates is probably to not use apply but to use the column array.

from pyproj import Proj
pp = Proj(proj='utm',zone=10,ellps='WGS84', preserve_units=False)

xx, yy = pp(My_data["LON"].values, My_data["LAT"].values)
My_data["X"] = xx
My_data["Y"] = yy 

Using Transformer

from pyproj import Transformer

trans = Transformer.from_crs(
    "epsg:4326",
    "+proj=utm +zone=10 +ellps=WGS84",
    always_xy=True,
)
xx, yy = trans.transform(My_data["LON"].values, My_data["LAT"].values)
My_data["X"] = xx
My_data["Y"] = yy

ORIGINAL ANSWER:

This answer here is great: https://gis.stackexchange.com/a/334276/144357

The solution below is for the purposes of understanding the root of the problem a bit better.

Your code in its current form re-constructs the Proj object with each iteration. This is a costly operation and is why the pyproj.Transformer object was created. It assists with repeated transformations because you don't have to re-create it each time (see: https://pyproj4.github.io/pyproj/stable/advanced_examples.html#repeated-transformations).

So, to avoid re-creating the Proj object, you can modify your code like so:

from pyproj import Proj
from functools import partial

p = Proj(proj='utm',zone=10,ellps='WGS84', preserve_units=False)

def impartial_rule(row, proj):
    x,y = proj(row["LON"], row["LAT"])
    return pd.Series({"X": x , "Y": y})

rule = partial(impartial_rule, proj=p)
My_data = My_data.merge(My_data.apply(rule, axis=1), left_index= True, right_index= True)

This should improve your performance.

Here is the equivalent using the pyproj.Transformer:

from pyproj import Transformer
from functools import partial

trans = Transformer.from_crs(
    "epsg:4326",
    "+proj=utm +zone=10 +ellps=WGS84",
    always_xy=True,
)

def impartial_rule(row, proj):
    x,y = proj(row["LON"], row["LAT"])
    return pd.Series({"X": x , "Y": y})

rule = partial(impartial_rule, proj=trans.transform)
My_data = My_data.merge(My_data.apply(rule, axis=1), left_index= True, right_index= True)

Hopefully this is helpful. Good luck!

Also, I would recommend reading this about Proj: https://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-longitude-to-projection-converter

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.