我正在对Geopanda库做一些工作,我在一个Excel文件中有一个包含多边形和shape的shapefile,然后将其转换为点。我想将两个DataFrame相交并将其导出到文件中。我还在两个投影(WGS84)上使用,以便可以对其进行比较。
至少应有一些与多边形相交的点。
我的相交的GeoSeries没有给我任何适合多边形的点,但是我不明白为什么...
我检查了shapefile的单位是否真的是公里,而不是其他东西。我不精通GeoPlot,因此无法真正确定GeoDataFrame的外观。
f = pd.read_excel(io = 'C:\\Users\\peilj\\meteo_sites.xlsx')
#Converting panda dataframe into a GeoDataFrame with CRS projection
geometry = [Point(xy) for xy in zip(df.geoBreite, df.geoLaenge)]
df = df.drop(['geoBreite', 'geoLaenge'], axis=1)
crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
gdf = GeoDataFrame(df, crs=crs, geometry=geometry)
#Reading shapefile and creating buffer
gdfBuffer = geopandas.read_file(filename = 'C:\\Users\\peilj\\lkr_vallanUTM.shp')
gdfBuffer = gdfBuffer.buffer(100) #When the unit is kilometer
#Converting positions long/lat into shapely object
gdfBuffer = gdfBuffer.to_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
#Intersection coordonates from polygon Buffer and points of stations
gdf['intersection'] = gdf.geometry.intersects(gdfBuffer)
#Problem: DOES NOT FIND ANY POINTS INSIDE STATIONS !!!!!!!
#Giving CRS projection to the intersect GeoDataframe
gdf_final = gdf.to_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
gdf_final['intersection'] = gdf_final['intersection'].astype(int) #Shapefile does not accept bool
#Exporting to a file
gdf_final.to_file(driver='ESRI Shapefile', filename=r'C:\\GIS\\dwd_stationen.shp
所需文件:
https://drive.google.com/open?id=11x55aNxPOdJVKDzRWLqrI3S_ExwbqCE9
最佳答案
两件事情:
创建指向以下内容的点时,需要交换geoBreite
和geoLaenge
:geometry = [Point(xy) for xy in zip(df.geoLaenge, df.geoBreite)]
这是因为在形状上遵循x,y逻辑,而不是lat,lon。
至于检查交叉点,您可以执行以下操作:
gdf['inside'] = gdf['geometry'].apply(lambda shp: shp.intersects(gdfBuffer.dissolve('LAND').iloc[0]['geometry']))
它会检测形状文件中的六个测站:
gdf['inside'].sum()
输出:
6
因此,连同其他一些小的修正,我们得到了:
import geopandas as gpd
from shapely.geometry import Point
df = pd.read_excel(r'C:\Users\peilj\meteo_sites.xlsx')
geometry = [Point(xy) for xy in zip(df.geoLaenge, df.geoBreite)]
crs = {'init': 'epsg:4326'}
gdf = gpd.GeoDataFrame(df, crs=crs, geometry=geometry)
gdfBuffer = gpd.read_file(filename = r'C:\Users\peilj\lkr_vallanUTM.shp')
gdfBuffer['goemetry'] = gdfBuffer['geometry'].buffer(100)
gdfBuffer = gdfBuffer.to_crs(crs)
gdf['inside'] = gdf['geometry'].apply(lambda shp: shp.intersects(gdfBuffer.dissolve('LAND').iloc[0]['geometry']))
关于python - 2 Geodataframe之间的交集,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/57430656/