I am trying to plot a point on top of a shapefile. the Coordinates of the point are accurate. I verified it using the go to XY tool of ArcMap. but when I try to plot overlay using Python the location of both is distorted. following is my code.
loc =gpd.read_file(r'C:\Users\MJ\Desktop\populationFunc\rentalRtes\resapart\gurgaonGISdata\locGISd.shp')
loc = loc.to_crs(epsg = 32643)
loc.plot(figsize = (10,10))
lat = 28.44198794960
long = 77.07818660830
p1 = Point((long, lat))
df = pd.DataFrame({'a':[lat,long]})
po = gpd.GeoDataFrame(geometry = [p1], crs = loc.crs)
po = po.to_crs(epsg = 32643)
po.plot()
fig,ax = plt.subplots(figsize = (15,7))
loc.plot(ax=ax, color = 'none', edgecolor = 'black')
po.plot(ax= ax)`
3 Answers 3
This is because you change the loc.crs
and then reference it. The following should work:
loc =gpd.read_file(r'C:\Users\MJ\Desktop\populationFunc\rentalRtes\resapart\gurgaonGISdata\locGISd.shp')
loc = loc.to_crs(epsg = 32643)
loc.plot(figsize = (10,10))
lat = 28.44198794960
long = 77.07818660830
p1 = Point((long, lat))
df = pd.DataFrame({'a':[lat,long]})
po = gpd.GeoDataFrame(geometry = [p1], crs = loc.crs)
loc = loc.to_crs(epsg = 32643)
loc.plot(figsize = (10,10))
po = po.to_crs(epsg = 32643)
po.plot()
fig,ax = plt.subplots(figsize = (15,7))
loc.plot(ax=ax, color = 'none', edgecolor = 'black')
po.plot(ax= ax)`
sidenote: Try to not use hard-coded paths in a string. Have a look at Pathlib
I dont't know the original projection of loc, nor the content of the GeoDataFrame, but if you fix the new projection of loc, the projection is epsg:32643
(in meters)
loc = loc.to_crs(epsg = 32643)
print(loc.crs)
epsg:32643
The unit of EPSG:32643 is meter and not degree (long,lat), therefore
lat = 28.44198794960 #in degree
long = 77.07818660830 #in degree
p1 = Point(long, lat)
po = gpd.GeoDataFrame(geometry = [p1],crs = "EPSG:4326") #in degrees
print(po)
geometry
0 POINT (77.07819 28.44199)
print(po.crs)
EPSG:4326
# conversion to meters
po = po.to_crs(epsg = 32643) # in meters
print(po)
geometry
0 POINT (703517.314 3147923.581)
print(po.crs)
EPSG:32643
With your solution
po = gpd.GeoDataFrame(geometry = [p1],crs = "EPSG:32643") #in meters
print(po)
geometry
0 POINT (77.07819 28.44199)
print(po.crs)
EPSG:32643
plot:
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(*po.geometry[0].xy,'o')
ax.axis('equal')
plt.show()
-
I have tried using GCS(WGS 84). plotting is accurate in that... but when I try using PCS(UTM Zone 43N) the output is distorted... I have used above bith methods... but not working for mealauddin– alauddin2021年02月11日 12:59:46 +00:00Commented Feb 11, 2021 at 12:59
-
look at my answer correctiongene– gene2021年02月11日 13:02:28 +00:00Commented Feb 11, 2021 at 13:02
Thank you, guys. the issue was with the projection. I changed the projection from GCS to PCS for all files and the issue resolved.