[GIS] Converting coordinates of a bunch of polygons in a shapefile which are in UTM to latitude and longitude

geopandaspython 3

My exact problem is that I have a file that report vehicular traffic that has 6 columns and two of these columns are latitude and longitude.for example this is the first row and 45. and 9. are latitude and longitude columns.

77c1b49236f442f37531551e2fc25e51    1430352000000000000 45.450824   9.211725    20  0
                         .
                         .
                         .

and on the other hand I have this shapefile of polygons

              id                                           geometry
0           3939_1_1  POLYGON ((504174.2673271392 5003118.268122713,...
1         3939_1_2_1  POLYGON ((508268.9938896392 5003118.268122713,...
2     3939_1_2_2_3_3  POLYGON ((511851.8796318267 5005064.825739901,...
3     3939_1_2_2_3_0  POLYGON ((511340.0388115142 5005064.825739901,...
4       3939_1_2_2_0  POLYGON ((510316.3571708892 5004415.973200838,...

with {'init': 'epsg:32632'}

As you see each polygon has an id I want to check each pair of latitude and longitude and see each one belongs to which polygon and add the column of polygon ids in my first file. So first I read my data file, then I converted the pair of latitude and longitude into a point and I created a geo data frame and then I read the shapefile and tried to check,
This is the piece of code that I have tried

import pandas as pd   
import geopandas as gpd
import geopandas.tools
from shapely.geometry import Point, mapping, shape
data_path = "/home/foroogh/PhD/milano-fcd/test"
filename = "2015-04-30-1430352000000000000"
mapDir = "/home/foroogh/PhD/milano-grid"
mapName = "intersection_Milano_W_GRIDIT_NEW.shp"
map_path = os.path.join(mapDir, mapName)
data = os.path.join(data_path, filename)
places = pd.read_csv(data, sep="\t",header = None, names =   ['id','dateandtime', 'latitude', 'longitude', 'vehiclecode', 'velocity'])
places = places[["id", "latitude", "longitude"]]
places["geometry"] = places.apply(lambda row: Point(row["longitude"],  row["latitude"]), axis = 1)
del(places["latitude"], places["longitude"])
places = gpd.GeoDataFrame(places, geometry = "geometry")
places.crs = {'init': 'epsg:32632'}
milano = gpd.GeoDataFrame.from_file(map_path)
print(milano)
milano = milano[["id", "geometry"]]
print(milano)
result = gpd.tools.sjoin(places, milano, how = "left")
print(result.head())

but I do not get the desired result and the problem is that my shapefile is in utm and my points in geo dataframe are in lat and lon and I did not know which one should be converted and how?

Best Answer

I had the same problem and realized geopandas being built on several other libraries makes it a bit convoluted.

What work for me was: 

#import some additional libraries: 
import pandas 
from pyproj import Proj, transform
from shapely.geometry import Point

#create your in and out projections with pyroj:
inProj = Proj({'init': 'epsg:xxxx'}) # use df.crs to get
outProj = Proj({'init': 'epsg:4326'})

#Iterate through your column of points and convert:
empty_list = []
for pt in geometry:
     #need to pull out data form Shapely object
     coords_obj = list(pt.coords)
      # transform points us pyroj
      x,y = transform(inProj,outProj, coords_obj[0][0], coords_obj[0][1])
      # put in your empty list as the desired Shapely object
      new_geo.append(Point(x,y))
 # convert to a pandas series
 empty_list =pandas.Series(empyt_list)
 # replace your geopandas column
 geometry = new_geo
Related Question