GeoPandas GeoDataFrame 操作

GeoPandas 0.6.0
基本操作です。
※見やすさのため不要なインデントをつけています。

空間演算は 以下を参照。
 spatial overlays
https://github.com/geopandas/geopandas/blob/master/examples/overlays.ipynb
Spatial Joins
https://github.com/geopandas/geopandas/blob/master/examples/spatial_joins.ipynb

import geopandas as gpd
#Ver. 0.6.3
 
#読み込み
  #SHP、GeoJSON、GeoPackage に対応 
  gdf = gpd.read_file('./A.shp')
  #フォルダ内のSHPファイルを全て読み込み結合
  import glob
  import pandas as pd
  shpfiles = [
   shpfile for shpfile in glob.glob('./shp/*.shp')]
  gdf = pd.concat(
   [gpd.read_file(shpfile)for shpfile in shpfiles])
  #df(csv)からgdfを作成
  #https://gis.stackovernet.com/ja/q/47396
  import pandas as pd
  from shapely.geometry import Point
  geometry = [Point(xy) for xy in zip(df.lon, df.lat)]
  crs = {'init': 'epsg:6668'}
  gdf = gpd.GeoDataFrame(df, crs=crs, geometry=geometry)
 
#書き込みSHP、GeoJSON、GeoPackage 
  #shpは仕様上、属性255種までの制限あり。型の制限もあり。
  #geojsonはQGISで読みこみ可。Arcはtoolboxから変換。
  #geojsonはdaytime型不可。→事前にstr型変換
  gdf.to_file('./A.shp')
  gdf.to_file('./A.geojson', driver='GeoJSON')
 
#型変換(pandasと同操作)
  gdf['時刻']=gdf['時刻'].astype(str)     
     
#A列（属性）削除 (pandasと同操作)
  gdf.drop('A',axis=1)
 
#列名（属性名）変更 (pandasと同操作)
  gdf.columns = ['A','B','geometry'] 
 
#'都道府県'列（属性）から3県を含む行抽出 (pandasと同操作)
  gdf3=gdf[(gdf['都道府県']=='茨城')|
           (gdf['都道府県']=='埼玉')|
           (gdf['都道府県']=='群馬')|]
 
#型変換(pandasと同操作)
  gdf['時刻']=gdf['時刻'].astype(str)     
     
#属性の追加、演算(pandasと同操作)
  gdf['C']=gdf['A']-gdf['B']
 
#県単位でディゾルブ（地物をまとめる）
  gdf=gdf3.dissolve(by='都道府県')
  
#属性結合
  #gdfにA列をkeyとしてDataFrameを属性結合
  #df.merge だと df になる
  gdf=gdf.merge(df, on='A')
 
#Coordinate Reference Systems
  #CRS確認
  gdf.crs
  #gdf1のCRSをgdf2のCRSに変更（一致）
  #EPSG:4612 : JGD2000
  #EPSG:6668 : JGD2011
  gdf1 = gdf1.to_crs(gdf2.crs)
  
#地図プロット
  #matplotlibベース
  #1レイヤー
  gdf.plot(figsize=(10, 10),
           alpha=0.5,
           edgecolor='k')
     
  #3レイヤー：gdfを重ねて表示  
  　#Choropleth Map(ポリゴン1枚目)
  　base=gdf1.plot(figsize=(10, 10),
      　  column='A',
          cmap='OrRd',
          alpha=0.8,
          legend=True,
          legend_kwds={'label':'aaa',
          orientation:'horizontal'}) 
    #ポリゴンを重ねる（2枚目）
    ax=gdf2.plot(ax=base,
                 facecolor='none',
                 edgecolor='gray',
         　　    alpha=0.8)                
  　#ポイントを重ねる（3枚目）
  　ax=gdf3.plot(ax=base,
          　　   marker='o',
          　　   color='red',
          　　   markersize=5)
　　#仕上げ
    ax.set_xlabel('longitude')
    ax.set_ylabel('latitude')
    ax.grid()
    ax.set_axisbelow(True) #grid moves back
    #保存する場合は import matplotlib.pyplot as plt からの
    plt.savefig("./aaa.png")
  
 

20200220 Foliumでの図化は移動しました。
https://phreeqc.blogspot.com/2020/02/folium.html


番外編： rasterstats を用いたラスターからの統計量取得。
大量処理は速度面で不可。

#空間演算
  from rasterstats import zonal_stats
  zonal_stats(gdf, './aaa.tif', stats=['max'])
  #フォルダ内のgeotiffを全て読み込み、統計量をgdfに結合
  import glob
  import os
  path = './input/*.tif'
  files= glob.glob(path)
  for i, name in enumerate(files) :
    print(i)
    file_name = os.path.splitext(
                            os.path.basename(files[i]))[0]
    print(file_name)
    stats=pd.DataFrame(zonal_stats(gdf,
                                   files[i],
                                   stats=['max', 'count']))
    gdf_max[file_name]=stats['max']
    gdf_count[file_name]=stats['count']