Automating GIS Processes 2024 写経　Exercise 2 (Problem 3、Lesson 2振り返り)


https://github.com/Automating-GIS-processes-II-2024/exercise-2/blob/main/Exercise-2-problem-3.ipynb

> During this task, the aim is to calculate the (air-line) distance in meters that each social media user in the data set prepared in Problem 2 has travelled in-between the posts. We’re interested in the Euclidean distance between subsequent points generated by the same user.

> For this, we will need to use the userid column of the data set kruger_posts.shp that we created in Problem 2.

> Answer the following questions:
>    - What was the shortest distance a user travelled between all their posts (in meters)?
>    - What was the mean distance travelled per user (in meters)?
>    - What was the maximum distance a user travelled (in meters)?

このタスクでは、air-line distance、直線距離を計算します。単位はmです。

[Problem 2](https://qiita.com/t_t_japan/items/630096364b47082063ca)のデータセットから、投稿と投稿の間に移動した距離、ユーザーごとのEuclidean distance、ユークリッド距離（空間内の2点間の直線距離）を求めることができます。
今回の問いは以下です。

- ユーザーの最小移動距離(m)は？
- ユーザーあたりの平均移動距離(m)はどのくらいか?
- ユーザーの最大移動距離(m)は？


# a) Read the input file and re-project it

> - Read the input file kruger_points.shp into a geo-data frame kruger_points
> Transform the data from WGS84 to an EPSG:32735 projection (UTM Zone 35S, suitable for South Africa). This CRS has metres as units.

```
import geopandas as gpd
import pathlib
DATA_DIRECTORY = pathlib.Path().resolve() / "data"

kruger_points=gpd.read_file(DATA_DIRECTORY / "kruger_points.shp")

# Check the data
kruger_points.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/6a2a16b9-22ff-4489-9ef0-98e03e3e346b.png)

kruger_points.shpを`GeoDataFrame`に読み込み、CRSをWGS84から南アフリカに適したEPSG:32735、UTM Zone 35Sに変換しましょう。 EPSG:32735（というか投影座標系）の距離の単位はmです。

```
# NON-EDITABLE CODE CELL FOR TESTING YOUR SOLUTION
kruger_points.crs = "EPSG:32735"

# Check that the crs is correct after re-projecting (should be epsg:32735)
import pyproj
assert kruger_points.crs == pyproj.CRS("EPSG:32735")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/47e04f8d-e9e1-4ebe-b130-0070bc91742c.png)

https://epsg.io/32735

https://docs.qgis.org/3.40/ja/docs/gentle_gis_introduction/coordinate_reference_systems.html#universal-transverse-mercator-utm-crs-in-detail

> 南アフリカは歪みを最小にするために4つの UTM ゾーン に覆われています。これらの ゾーン は、 UTM 33S 、 UTM 34S 、 UTM 35S 、 UTM 36S と呼ばれています。ゾーンの後の S は、UTMゾーンが 赤道より南 に位置していることを意味します。

# b) Group the data by user id

> Group the data by userid and store the grouped data in a variable grouped_by_users

データをユーザーIDでグループ化し、変数`grouped_by_users`に保存します。
また、グループの数と`nunique()`メソッドで求めたユニークなユーザーIDの数が一致しているか、assertで確認します。
```
grouped_by_users = kruger_points.groupby("userid")

# Check the number of groups:
assert len(grouped_by_users.groups) == kruger_points["userid"].nunique(), "Number of groups should match number of unique users!"
```

# c) Create shapely.geometry.LineString objects for each user connecting the points from oldest to most recent

ユーザーごと、`Point`群から（時系列が古いデータから新しいデータに並びかえて）`LineString`オブジェクトを作成しましょう。

> There are multiple ways to solve this problem (see the hints for this exercise. You can use, for instance, a dictionary or an empty GeoDataFrame to collect data that is generated using the steps below:

> - Use a for-loop to iterate over the grouped object. For each user’s data:
>   - sort the rows by timestamp
>   - create a shapely.geometry.LineString based on the user’s points

> **CAREFUL**: Remember that every LineString needs at least two points. Skip users who have less than two posts.

やりかたはいくつかあります。（[Hint](https://qiita.com/t_t_japan/items/c55c30931bc884265da5)を参照）

たとえば、`dictonary`や空の`GeoDataFrame`を使用するやり方です。グループごとにループで処理、ループ内ではタイムスタンプでソートし、`Point`群から`LineString`オブジェクトを生成します。

**気を付けるべき点**としては、`LineString`は少なくとも`Point`が２つ必要だということです。つまり`Point`が一つしかないユーザー（＝移動してないユーザー）の場合、`LineString`の生成をスキップする必要があります。

> Store the results in a geopandas.GeoDataFrame called movements, and remember to assign a CRS.

では、`GeoDataFrame`型の変数`movements`に格納した結果を見てみましょう。CRSの設定は忘れずに。

```
from shapely.geometry import LineString

movements_list = []

for key, group in grouped_by_users:

     points = group.sort_values(by="timestamp").geometry
     # Pointが一つの場合はスキップ
     if len(points) > 1:
        movements_list.append(LineString(points))

movements  = gpd.GeoDataFrame(movements_list, columns=["geometry"])
movements.crs = "EPSG:32735"

# Check the result
print(type(movements))
print(movements.crs)

movements
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/e0302be3-a4f6-4389-9ba3-8c6c1945637c.png)




# d) Calculate the distance between all posts of a user

> - Check once more that the CRS of the data frame is correct
> - Compute the lengths of the lines, and store it in a new column called distance

CRSが正しいか確認します。また、`LineString`オブジェクトの距離を計算し、その結果を`distance`列に格納します。

```
movements["distance"]= movements.apply(lambda row: row.geometry.length, axis=1)

# Check the result
movements.head()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/82cf347f-0fc5-45bc-9133-e8ece8ffdb21.png)

# e) Answer the original questions

> You should now be able to quickly find answers to the following questions:

> - What was the shortest distance a user travelled between all their posts (in meters)? (store the value in a variable shortest_distance)
> - What was the mean distance travelled per user (in meters)? (store the value in a variable mean_distance)
> - What was the maximum distance a user travelled (in meters)? (store the value in a variable longest_distance)

ここまでくれば、問いに対する答えを簡単に出すことができます。

- ユーザーの最小移動距離(m)は？ →変数`shortest_distance`に
- ユーザーあたりの平均移動距離(m)はどのくらいか?　→変数`mean_distance`に
- ユーザーの最大移動距離(m)は？　→変数`longest_distance`に

```
shortest_distance = movements["distance"].min()
mean_distance = movements["distance"].mean()
longest_distance = movements["distance"].max()

print(f"shortest_distance: {shortest_distance}")
print(f"mean_distance: {mean_distance}")
print(f"longest_distance: {longest_distance}")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/27c2ad07-e3e8-4342-8b55-d168075fe8e9.png)

最小は0m（移動していない）
平均は約1m
最大は約64mという結果になりました。

本当でしょうか。

まず、describe()メソッドで統計情報を見てみましょう。

```
movements.describe()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/063251e1-25c9-4424-8535-5a8dbdb1e82f.png)

50%は中央値（median）です。それが0.14m。つまりほとんど移動していません。
何となく想像はできますが、ヒストグラムにしてみましょう。最大値が約64mなので、ビン（棒の幅）はざっくりと65にしています。

```
movements.hist(column="distance", bins=65)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/a5e82ddd-cfe9-411e-b68b-6d343b86463b.png)

うん、ほとんどの人が移動してないですね。


# f) Save the movements in a file

> Save the movements into a new Shapefile called movements.shp inside the data directory.

変数`movements`をmovements.shpファイルとして保存します。

```
movements.to_file(DATA_DIRECTORY / "movements.shp")
```

ファイルが書き込めているかどうかを、shellコマンドで確認します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/23ceb057-d050-411a-a210-374a65238a59.png)

これで終わりです。

# Lesson 2振り返り

https://autogis-site.readthedocs.io/en/latest/lessons/lesson-2/overview.html

> Learning goals
> After this week, you should be able to:
> - read and write spatial data from and to common file formats,
> - filter and re-group data by spatial and non-spatial characteristics, and
> - manage and transform a data set’s coordinate reference system.

- 一般的なファイルフォーマットで空間データを読み書きします →できました
- データをフィルタしたり、再度グループ化します →できました
- データの座標参照系を変換します →できました 

次は Lession 3です。

5月ひと月で、Lesson1とLesson2を消化しました。
が、内容はだんだん難しくなっていくので、Lession 3は6月いっぱいかかるかもしれません。