東京23区で一番緑が多いのはどこ
特に東京に行ったこともないのですが、東京で一番緑が多いのはどこかがふと気になった次第です。将来引っ越しすることがあったら参考になるかもしれないので調べます。
やりかた
- 行政区域データを集める
- 土地利用データを集める
- 各行政区域について土地利用ごとの面積を算出
- 各行政区域について、自然っぽい土地利用の割合を算出
当然このようにすればできます。文字にすると簡単そうです。
データ収集
行政区域データは国土数値情報ダウンロードサイトから取得します。平成29年版の東京都のジオメトリデータで、ESRI Shapefile形式をダウンロードします。
土地利用データはALOSの高解像度土地利用土地被覆図の10mメッシュのgeotiffをダウンロードし、必要そうなところだけ展開します。
収集したデータは以下になります。geotiffだけで3GBくらいあります。
$ tree . ├── geotiff │ ├── LC_N20E136.tif │ ├── LC_N24E141.tif │ ├── LC_N24E153.tif │ ├── LC_N25E141.tif │ ├── LC_N26E142.tif │ ├── LC_N27E140.tif │ ├── LC_N27E142.tif │ ├── LC_N30E140.tif │ ├── LC_N32E139.tif │ ├── LC_N33E136.tif │ ├── LC_N33E139.tif │ ├── LC_N34E136.tif │ ├── LC_N34E137.tif │ ├── LC_N34E138.tif │ ├── LC_N34E139.tif │ ├── LC_N35E136.tif │ ├── LC_N35E137.tif │ ├── LC_N35E138.tif │ ├── LC_N35E139.tif │ ├── LC_N35E140.tif │ └── ver1609_LC_GeoTiff.tgz └── shapefiles ├── KS-META-N03-17_13_170101.xml ├── N03-17_13_170101.dbf ├── N03-17_13_170101.prj ├── N03-17_13_170101.shp ├── N03-17_13_170101.shx └── N03-17_13_170101.xml
DockerでPostGISサーバを立てる
行政区域の土地利用ごとの面積を計算するには、行政区域のジオメトリを土地利用ごとに切り出したジオメトリを用意する必要があります。 ジオメトリ演算といえばPostGISを使うようです。環境構築は面倒くさそうなので、dockerを使います。
以下のようなDockerfileを用意します。
FROM mdillon/postgis:latest RUN rm -rf /docker-entrypoint-initdb.d/postgis.sh EXPOSE 5432
/docker-entrypoint-initdb.d/にはコンテナ起動直後に実行するポスグレDB初期化スクリプトが配置されます。postgis.shはmdillon/postgisが用意したPostGIS初期設定スクリプトなのですが、不要なextensionとかもインストールするので削除します。 ちなみに/docker-entrypoint-initdb.d/内のスクリプトは名前順に実行されるっぽいので、適当な名前でオレオレスクリプトを配置すると実行順序が狂って死んだりするので注意が必要です。
$ docker build --no-cache -t postgis_calcgreen . $ docker run \ -e POSTGRES_DB=postgres \ -e POSTGRES_USER=postgres \ -e POSTGRES_PASSWORD=postgres \ -v "$PWD":/usr/src/myapp \ -w /usr/src/myapp \ --name postgis_calcgreen \ -p 5432:5432 \ postgis_calcgreen
イメージつくってコンテナ起動します。パスワード等はテキトーです。/docker-entrypoint-initdb.d/のスクリプトがコンテナ起動直後に走るので、それを確認するために"-d"は付けないのが無難です。
動作確認します。psqlコマンド使います。postgresqlでは"\l"でデータベース一覧が取得できます。
$ export PGHOST=localhost
$ export PGDATABASE=postgres
$ export PGUSER=postgres
$ export PGPASSWORD=postgres
$ psql -c "\l"
List of databases
Name | Owner | Encoding | Collate | Ctype | Access privileges
-----------+----------+----------+------------+------------+-----------------------
postgres | postgres | UTF8 | en_US.utf8 | en_US.utf8 |
template0 | postgres | UTF8 | en_US.utf8 | en_US.utf8 | =c/postgres +
| | | | | postgres=CTc/postgres
template1 | postgres | UTF8 | en_US.utf8 | en_US.utf8 | =c/postgres +
| | | | | postgres=CTc/postgres
(3 rows)
ポスグレサーバが立ち上がりました。
データ投入とか
続いてpostgis extensionを追加します。
$ psql -c "create extension if not exists postgis;"
CREATE EXTENSION
$ psql -c "\dx;"
List of installed extensions
Name | Version | Schema | Description
---------+---------+------------+---------------------------------------------------------------------
plpgsql | 1.0 | pg_catalog | PL/pgSQL procedural language
postgis | 2.3.2 | public | PostGIS geometry, geography, and raster spatial types and functions
(2 rows)
できました。
行政区域シェープファイルはdockerコンテナ上でshp2pgsqlコマンドを使ってsqlに変換し、tmpという名のテーブルに登録します。-sはSRID、-iはInt4指定、-Iはgistインデックス作成、-Wは文字コード指定です。sqlに変換したら直接ローカルのpsqlコマンドに流して実行します。
tmpテーブルのカラムを調整し、全レコードを市区町村でgroup byしてST_Unionしてjapan_ksテーブルに登録します。tmpテーブルはお役御免とdropします。
$ docker exec -i $PGCONTAINER bash -c \
'shp2pgsql -s 4326 -D -i -I -W cp932 shapefiles/*.shp tmp' | psql >/dev/null
$ cat << EOS | psql
delete from tmp where n03_007 is NULL;
alter table tmp rename n03_001 to prefecture_name;
alter table tmp rename n03_007 to city_code;
alter table tmp add prefecture_code varchar(2);
update tmp set prefecture_code = substring(city_code,1,2);
alter table tmp add city_name varchar(60);
update tmp set city_name = coalesce(n03_003, '') || coalesce(n03_004, '');
create table japan_ks as
select
cast(prefecture_code as integer) as prefecture_code
, prefecture_name
, cast(city_code as integer) as city_code
, city_name
, sum(ST_Area(geography(geom))) as area
, ST_Union(geom) as geom
from tmp
group by
prefecture_code
, prefecture_name
, city_code
, city_name
;
drop table tmp;
create index on japan_ks using gist (geom);
analyze japan_ks;
EOS
次にalos土地利用データを登録します。raster2pgsqlコマンドを使います。-Iはインデックス作成、-Cは「SRIDやピクセルサイズ等のラスタ制約を適用して、raster_columnsビューで適切な登録ができるようにします。 」らしいです。-Nはヌル値指定、-tはgeotiffをレコード単位で分割するときの、1レコード当たりのタイルサイズです。タイルサイズが大きいとジオメトリ演算で時間がかかりそうですが、小さいと登録に時間がかかります。
$ ALOS_GEOTIF=./geotiff/LC*.tif $ docker exec -i $PGCONTAINER bash <<EOF | psql > /dev/null raster2pgsql -I -C -s 4326 -N 0 -t 100x100 \ $ALOS_GEOTIF \ alos_rast EOF
create table alos_code_name ( code integer primary key, name varchar(40) ); insert into alos_code_name (code, name) values (0, '未分類'), (1, '水域'), (2, '都市'), (3, '水田'), (4, '畑地'), (5, '草地'), (6, '落葉広葉樹'), (7, '落葉針葉樹'), (8, '常緑広葉樹'), (9, '常緑針葉樹'), (10, '裸地'), (11, '雪氷'), (253, 'その他'), (255, 'データなし') ;
alos_rastに格納したラスタデータをジオメトリに変換してalos_geoに格納します。行政区域ジオメトリと交差する土地利用ラスタのみ変換します。
create table alos_geo as select alos_geo.rid as rid , (alos_geo.geomval).val as code , (alos_geo.geomval).geom as geom from (select alos_rast.rid as rid , ST_DumpAsPolygons(alos_rast.rast) as geomval from (select ST_Union(geom) as geom from japan_ks) as whole_japan , alos_rast where ST_Intersects(whole_japan.geom, alos_rast.rast) ) as alos_geo ; create index on alos_geo using gist (geom); analyze alos_geo;
japan_ksとalos_geoのそれぞれのジオメトリのintersection(交差ジオメトリ)をjapan_alosテーブルに格納します。今回はalosのラスタをジオメトリに変換して(alos_geoテーブル)ST_Intersectionの引数に与えていますが、ST_Intersectionはラスタを直接引数にとることもできます(その場合は裏でラスタデータのST_DumpAsPolygonsが走る)。
create table japan_alos as select japan_ks.city_code as city_code , japan_ks.city_name as city_name , alos_geo.code as alos_code , ST_Intersection(japan_ks.geom, alos_geo.geom) as geom from japan_ks , alos_geo where ST_Intersects(japan_ks.geom, alos_geo.geom) ;
結果
準備できました。可視化したり緑の多いところを調べたりします。Pythonを使います。GISの可視化といえばgeopandasなのだそうです。
%matplotlib inline # ref. https://geohackweek.github.io/vector/06-geopandas-advanced/ import psycopg2 import matplotlib as mpl import matplotlib.pyplot as plt plt.style.use('bmh') import pandas as pd import geopandas as gpd
適当にmatplotlibの設定をします。
mpl.rc('font', family='Osaka') # Macだとこれで日本語が文字化けしないかもしれない mpl.rc('legend', loc="lower left", fontsize='large') # geopandas v0.3.0のlegend.locはbestでハードコードされているのでrcで指定しても無意味
PostGISサーバのアクセス情報です。
db_conn_dict = {
"host":"localhost",
"port":5432,
"dbname":"postgres",
"user":"postgres",
"password":"postgres"
}
試しに行政区域ジオメトリを可視化します。
with psycopg2.connect(**db_conn_dict) as conn: japan_ks = gpd.read_postgis( "SELECT * FROM japan_ks", conn, geom_col='geom', crs={'init': u'epsg:4326'} ) japan_ks.plot(column='city_code', cmap='tab20c', categorical=True, figsize=(14, 8)) plt.savefig("japan_ks.png", bbox_inches="tight", pad_inches=0.0, transparent=True
東京は広いですね。よくわかりません。島嶼を除きます。
japan_ks.plot(column='city_name', cmap='tab20c', categorical=True, figsize=(14, 8)) plt.ylim([35.45,35.95]) plt.xlim([138.85,140]) for idx, row in japan_ks.iterrows(): plt.annotate( s=row.city_name, xy=row.geom.representative_point().coords[0], horizontalalignment='center' ) plt.savefig("japan_ks.zoom.png", bbox_inches="tight", pad_inches=0.0)
ラベルが少し見にくいですがそれらしいような気がします。なぜかbbox_inches="tight"が無効になっています。
土地利用別の可視化もしてみます。つまりjapan_alosテーブルを可視化します。下はカラーマップです。
alos_color_palette = {
'00_未分類': [0 /255 ,0 /255 ,0 /255 ],
'01_水域': [0 /255 ,0 /255 ,100/255 ],
'02_都市': [255/255 ,0 /255 ,0 /255 ],
'03_水田': [0 /255 ,128/255 ,255/255 ],
'04_畑地': [255/255 ,193/255 ,191/255 ],
'05_草地': [255/255 ,255/255 ,0 /255 ],
'06_落葉広葉樹': [128/255 ,255/255 ,0 /255 ],
'07_落葉針葉樹': [0 /255 ,255/255 ,128/255 ],
'08_常緑広葉樹': [86 /255 ,172/255 ,0 /255 ],
'09_常緑針葉樹': [130/255 ,172/255 ,117/255 ],
'10_裸地': [128/255 ,100/255 ,0 /255 ]
}
def alos_colors(alos_names):
return mpl.colors.ListedColormap(
[alos_color_palette[x] for x in alos_names],
name="alos_colors",
N=len(alos_names))
query = """ select japan_alos.geom as geom , japan_alos.city_name as city_name , alos_code_name.code as alos_code , to_char(alos_code_name.code, 'FM00') || '_' || alos_code_name.name as alos_name from japan_alos left join alos_code_name on (japan_alos.alos_code = alos_code_name.code) ; """ with psycopg2.connect(**db_conn_dict) as conn: japan_alos = gpd.read_postgis( query, conn, geom_col='geom', crs={'init': u'epsg:4326'} )
japan_alos.plot(
column='alos_name',
cmap=alos_colors(sorted(japan_alos.alos_name.unique())),
categorical=True,
legend=True,
figsize=(14, 8)
)
plt.ylim([35.45,35.95])
plt.xlim([138.85,140.25])
plt.savefig("japan_alos.zoom.png", bbox_inches="tight", pad_inches=0.0)
ひとまず東京都できちんと切り出されていることは確認できました。
さて、市区町村ごとのalosの落葉広/落葉針/常緑広/常緑針葉樹の割合を算出します。
greenratio_query = """ select japan_alos.city_code , min(japan_alos.city_name) as city_name , sum(ST_Area(geography(japan_alos.geom))) / min(japan_ks.area) as green_ratio from japan_alos left join alos_code_name on (japan_alos.alos_code = alos_code_name.code) left join japan_ks on (japan_alos.city_code = japan_ks.city_code) where alos_code between 6 and 9 group by japan_alos.city_code order by green_ratio desc, japan_alos.city_code ;""" with psycopg2.connect(**db_conn_dict) as conn: greenratio = pd.io.sql.read_sql(greenratio_query, conn)
top10です。西多摩郡か島嶼です。西多摩郡檜原村というところがほとんど森です。
greenratio.head(n=10)
| city_code | city_name | green_ratio | |
|---|---|---|---|
| 0 | 13307 | 西多摩郡檜原村 | 0.992630 |
| 1 | 13308 | 西多摩郡奥多摩町 | 0.965666 |
| 2 | 13382 | 御蔵島村 | 0.915994 |
| 3 | 13361 | 大島町 | 0.836055 |
| 4 | 13362 | 利島村 | 0.832088 |
| 5 | 13305 | 西多摩郡日の出町 | 0.825971 |
| 6 | 13205 | 青梅市 | 0.763871 |
| 7 | 13364 | 神津島村 | 0.746610 |
| 8 | 13363 | 新島村 | 0.743949 |
| 9 | 13228 | あきる野市 | 0.735391 |
可視化してみますが、やはりほとんど森です。
query = """ SELECT japan_alos.geom as geom , japan_alos.city_name as city_name , alos_code_name.code as alos_code , to_char(alos_code_name.code, 'FM00') || '_' || alos_code_name.name as alos_name FROM japan_alos left join alos_code_name on (japan_alos.alos_code = alos_code_name.code) where city_name = '西多摩郡檜原村' ; """ with psycopg2.connect(**db_conn_dict) as conn: hinoharamura = gpd.read_postgis( query, conn, geom_col='geom', crs={'init': u'epsg:4326'} ) hinoharamura.plot( column='alos_name', cmap=alos_colors(sorted(hinoharamura.alos_name.unique())), categorical=True, legend=True, figsize=(14, 8) ) plt.xlim([139,139.22]) # plt.gca().get_xaxis().get_major_formatter().set_useOffset(False) plt.savefig("hinoharamura.png", bbox_inches="tight", pad_inches=0.0)
東京23区のランキングです。千代田区が1位ですね。
greenratio[greenratio.city_code < 13200].reset_index(drop=True)
| city_code | city_name | green_ratio | |
|---|---|---|---|
| 0 | 13101 | 千代田区 | 0.110995 |
| 1 | 13113 | 渋谷区 | 0.092015 |
| 2 | 13103 | 港区 | 0.075061 |
| 3 | 13105 | 文京区 | 0.061869 |
| 4 | 13112 | 世田谷区 | 0.052596 |
| 5 | 13120 | 練馬区 | 0.046121 |
| 6 | 13115 | 杉並区 | 0.045913 |
| 7 | 13104 | 新宿区 | 0.037274 |
| 8 | 13109 | 品川区 | 0.035927 |
| 9 | 13110 | 目黒区 | 0.034893 |
| 10 | 13106 | 台東区 | 0.031348 |
| 11 | 13119 | 板橋区 | 0.028702 |
| 12 | 13108 | 江東区 | 0.024454 |
| 13 | 13117 | 北区 | 0.019062 |
| 14 | 13116 | 豊島区 | 0.016745 |
| 15 | 13111 | 大田区 | 0.016584 |
| 16 | 13102 | 中央区 | 0.016507 |
| 17 | 13114 | 中野区 | 0.016115 |
| 18 | 13122 | 葛飾区 | 0.015348 |
| 19 | 13123 | 江戸川区 | 0.013803 |
| 20 | 13121 | 足立区 | 0.008889 |
| 21 | 13107 | 墨田区 | 0.005948 |
| 22 | 13118 | 荒川区 | 0.003172 |
千代田区も同様に可視化してみました。緑を囲むお堀のような構造が見えます。これはつまり旧江戸城現皇居ですね。
皇居はずるいように思えたので、23区で2番目に緑が多い渋谷区です。こちらも緑色のオーバルが見られますが、これは代々木公園のようです。
というわけで23区で一番緑が多いのは千代田区、皇居は自由に入れないので除外すべきと思えば渋谷区となりました。
