Automating GIS Processes 2024 写経　Lesson 2 （Geopandas: an introduction）

https://autogis-site.readthedocs.io/en/latest/lessons/lesson-2/geopandas-an-introduction.html

> In this section, we will cover the basics of geopandas, a Python library to interact with geospatial vector data.

> Geopandas provides an easy-to-use interface to vector data sets. It combines the capabilities of pandas, the data analysis package we got to know in the Geo-Python course, with the geometry handling functionality of shapely, the geo-spatial file format support of fiona and the map projection libraries of pyproj.

`geopandas`によるベクターデータな基本的な処理をやっていきましょう。
`geopandas`は、[pandas](https://www.google.com/url?sa=t&source=web&rct=j&opi=89978449&url=https://pandas.pydata.org/&ved=2ahUKEwjS98Xhi66NAxU5rlYBHYw4BW8QFnoECBoQAQ&usg=AOvVaw3cD5ulu4AnZcNusojIyttY)に`shapely`・`fiona`・`pyproj`（これは別途）などを組み込んだものです。


> There is one key difference between pandas’s data frames and geopandas’ GeoDataFrames: a GeoDataFrame contains an additional column for geometries. By default, the name of this column is geometry, and it is a GeoSeries that contains the geometries (points, lines, polygons, …) as shapely.geometry objects.

`pandas`との違いは、`GeoDataFrames`であること。`GeoDataFrames`は、`GeoSeries`という型で、`shapely.geometry`オブジェクトを格納します。


まずは開始する前のお約束として、`geopandas`をインストールします。
```
! pip install geopandas
```
次に、読み込むファイル（m_L4132R_p.shp）をアップロードします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/356b48e1-f7b1-4269-9a95-d2fea7a8d4e7.png)

ファイルは　https://github.com/mopinfish/auto-gis-notebooks/tree/main/lessons/lesson-2/data/finland_topographic_database にあります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/095f8f34-45bd-4e7a-b9b4-4defbf8a4749.png)


```
import pathlib
import geopandas
import numpy
import pandas

DATA_DIRECTORY = pathlib.Path().resolve() / "data"
print (DATA_DIRECTORY)
HIGHLIGHT_STYLE = "background: #f66161;"

df = geopandas.read_file(DATA_DIRECTORY / "finland_topographic_database" / "m_L4132R_p.shp")

df["geom"] = df.geometry.to_wkt().apply(lambda wkt: wkt[:40] + " ...")

df = df[["RYHMA", "LUOKKA", "geom"]]
df = df.rename(columns={"RYHMA": "GROUP", "LUOKKA": "CLASS", "geom": "geometry"})

(
    df.head().style
        .applymap(lambda x: HIGHLIGHT_STYLE, subset=["geometry"])
        .apply_index(lambda x: numpy.where(x.isin(["geometry"]), HIGHLIGHT_STYLE, ""), axis=1)
)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/6affc120-03c9-4788-a8f9-a10da780b464.png)


コードが長めですが、あまり気にしなくてよいかと。
肝心なのは、`geometry`列の中身が、`POLYGON(....)`となっていることです。
確かに、`shapely.geometry`オブジェクトになっています。


# Input data: Finnish topographic database

> In this lesson, we will work with the National Land Survey of Finland (NLS)/Maanmittauslaitos (MML) topographic database. - The data set is licensed under the NLS’ open data licence (CC BY 4.0). - The structure of the data is described in a separate Excel file. - Further information about file naming is available at fairdata.fi (this link relates to the 2018 issue of the topographic database, but is still valid).

> For this lesson, we have acquired a subset of the topographic database as shapefiles from the Helsinki Region in Finland via the CSC’s Paituli download portal. You can find the files in data/finland_topographic_database/.

> The Paituli spatial download service offers data from a long list of national institutes and agencies.

以降で使うデータの源泉についてです。興味があれば読み解いてもらえれば。


# Read and explore geo-spatial data sets

> Before we attempt to load any files, let’s not forget to defining a constant that points to our data directory:

ファイルを読み込むディレクトリのフルパスを定義しておきます。

```
import pathlib
NOTEBOOK_PATH = pathlib.Path().resolve()
DATA_DIRECTORY = NOTEBOOK_PATH / "data"
```

> In this lesson, we will focus on terrain objects (Feature group: “Terrain/1” in the topographic database). The Terrain/1 feature group contains several feature classes.

> Our aim in this lesson is to save all the Terrain/1 feature classes into separate files.

この後は、“Terrain/1”という地形オブジェクト（でしょうか、自信なし）のデータを使用します。

ひとまずは、いろいろな地物（例：池、公園）を`feature class`という区分で格納している、と理解しておけばよいでしょう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/0209a73c-5752-4907-9de7-5260c85146ab.png)

> If you take a quick look at the data directory using a file browser, you will notice that the topographic database consists of many smaller files. Their names follow a strictly defined convention, according to this file naming convention, all files that we interested in (Terrain/1 and polygons) start with a letter m and end with a p.

過程はすっ飛ばします。以下のファイルをアップロードします。
m_L4132R_p.cpg
m_L4132R_p.dbf
m_L4132R_p.prj
m_L4132R_p.shp
m_L4132R_p.shx

ファイルは、https://github.com/mopinfish/auto-gis-notebooks/tree/main/lessons/lesson-2/data/finland_topographic_database にあります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/994a1a17-75a7-4cbd-9114-f1c3019f347b.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/21972b01-6ca8-411b-a7a3-2826da6bcf13.png)

そして、上記のファイルを読み込みます。

```
import geopandas
TOPOGRAPHIC_DATABASE_DIRECTORY = DATA_DIRECTORY / "finland_topographic_database"
input_filename = list(TOPOGRAPHIC_DATABASE_DIRECTORY.glob("m*p.shp"))[0]
data = geopandas.read_file(input_filename)
data.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/781ff163-2fde-4981-8d0c-37f328af64bc.png)

列を見てみましょう。
```
data.columns
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/f00a5f63-6f42-4956-a2f9-d55e386c6837.png)

列でフィルターします。
サンプルコード違い、フィルター後のGeoDataFrameの名前は`data_modified`としています。
```
data_modified = data[["RYHMA", "LUOKKA", "geometry"]]
data_modified.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/46f77466-3011-453d-b19b-65a016b345e3.png)

列名を変更します。
```
data_modified = data_modified.rename(
    columns={
        "RYHMA": "GROUP",
        "LUOKKA": "CLASS"
    }
)
data_modified.head()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/206737c4-1a91-4ad3-bdc5-36cb5d2a7fa4.png)

### Explore the data set in a map:

> As geographers, we love maps. But beyond that, it’s always a good idea to explore a new data set also in a map. To create a simple map of a geopandas.GeoDataFrame, simply use its plot() method. It works similar to pandas (see Lesson 7 of the Geo-Python course, but draws a map based on the geometries of the data set instead of a chart.

`plot()`メソッドで`geopandas.GeoDataFrame`のデータを地図として表示できます。

```
data_modified.plot()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/f56effd6-fe80-4a6a-b493-1cd04db8b362.png)

### Geometries in geopandas

> Geopandas takes advantage of shapely’s geometry objects. Geometries are stored in a column called geometry.

`shapely`の`geometry`オブジェクトは、`geometry`と呼ばれる列に格納されています。

```
data_modified.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/5b36fd1d-fe35-457d-a1dc-214c272b551a.png)

> Lo and behold, the geometry column contains familiar-looking values: Well-Known Text (WKT) strings. Don’t be fooled, they are, in fact, shapely.geometry objects (you might remember from last week’s lesson) that, when print()ed or type-cast into a str, are represented as a WKT string).

表示はWell-Known Text (WKT) ですが、実際は`shapely`の`geometry`オブジェクトです。

> Since the geometries in a GeoDataFrame are stored as shapely objects, we can use shapely methods to handle geometries in geopandas.

> Let’s take a closer look at (one of) the polygon geometries in the terrain data set, and try to use some of the shapely functionality we are already familiar with. For the sake of clarity, first, we’ll work with the geometry of the very first record, only:

`at()`メソッドで１レコード目の`geometry`を参照し、shapelyとしての機能を確認しましょう。

```
# The value of the column `geometry` in row 0:
data_modified.at[0, "geometry"]
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/a3b5b979-8893-40af-8dce-bb94a1e5e307.png)

`area`プロパティから面積を求めることができます。
```
print(f"Area: {round(data_modified.at[0, 'geometry'].area)} m².")
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/798d33b7-5fe0-4449-ac3f-27c5e6c1d4b7.png)

> Let’s do the same for multiple rows, and explore different options of how to. First, use the reliable and tried iterrows() pattern we learned in lesson 6 of the Geo-Python course.

`iterrows() `メソッドで複数レコードまとめて処理しましょう。
```
for index, row in data_modified[:5].iterrows():
    polygon_area = row["geometry"].area
    print(f"The polygon in row {index} has a surface area of {polygon_area:0.1f} m².")
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/69492de2-6be3-4108-9bb7-25b85b9b18b6.png)

> As you see, all pandas functions, such as the iterrows() method, are available in geopandas without the need to call pandas separately. Geopandas builds on top of pandas, and it inherits most of its functionality.

> Of course the iterrows() pattern is not the most convenient and efficient way to calculate the area of many rows. Both GeoSeries (geometry columns) and GeoDataFrames have an area property:

`pandas`の機能は`geopandas`でも同じように使えます。

ただ、複数レコードで、`area`プロパティ処理するときに、`iterrows() `メソッドは最も有効な方法ではありません。

`GeoSeries`や`GeoDataFrames`そのものが`area`を持っています。
```
data_modified.area
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/b1f04fbb-fd74-44e2-8d72-492825f11b31.png)

上記は以下と同じ結果になります。
```
data_modified["geometry"].area
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/cb28d8b6-50d6-46cc-a399-59f8fb637835.png)

`area`プロパティを別の列に外だしすることもできます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/259a2b63-2780-4cd6-878c-f47595ac12ad.png)

長くなったのでここまでとします。続きます。