Automating GIS Processes 2024 写経　Lesson 3（Point-in-polygon queries）前半

https://autogis-site.readthedocs.io/en/latest/lessons/lesson-3/point-in-polygon-queries.html

# Point-in-polygon queries

> Finding out if a certain point is located inside or outside of an area, or finding out if a line intersects with another line or polygon are fundamental geospatial operations that are often used e.g. to select data based on location. Such spatial queries are one of the typical first steps of the workflow when doing spatial analysis. Performing a spatial join (will be introduced later) between two spatial datasets is one of the most typical applications where Point in Polygon (PIP) query is used.

> For further reading about PIP and other geometric operations, see Chapter 4.2 in Smith, Goodchild & Longley: Geospatial Analysis - 6th edition.

点があるエリアの内側か外側か、線分がほかの線やポリゴンと交差しているか否か、これらは基本的な地理空間操作で、特定の場所のデータを選択する際に使われます。

Point in Polygon(PIP)クエリは、２つの地理空間データを空間結合の典型例です。


## How to check if point is inside a polygon?

> Computationally, detecting if a point is inside a polygon is most commonly done using a specific formula called Ray Casting algorithm. Luckily, we do not need to create such a function ourselves for conducting the Point in Polygon (PIP) query. Instead, we can take advantage of Shapely’s binary predicates that can evaluate the topolocical relationships between geographical objects, such as the PIP as we’re interested here.

対象の点がとあるポリゴンの内側にあることを検出する方法として、[レイキャスティングアルゴリズム](https://en.wikipedia.org/wiki/Point_in_polygon#Ray_casting_algorithm)があります。

https://ja.scratch-wiki.info/wiki/%E3%83%AC%E3%82%A4%E3%82%AD%E3%83%A3%E3%82%B9%E3%83%86%E3%82%A3%E3%83%B3%E3%82%B0

ただ、このアルゴリズムを自分で実装する必要はありません。[Shapelyのbinary predicates](https://shapely.readthedocs.io/en/stable/manual.html#binary-predicates)、2つのものの関係をあらわす表現（←生成AIに聞いてみました）を使って、topological relationships、空間データの関係の仕方（←これも生成AIに聞いてみました）を判断することができます。


::: note
### Good to know: Ray-Casting Algorithm

> The ray-casting algorithm is a common method used to determine whether a given point lies inside or outside of a polygon. This algorithm is widely used in spatial analysis and GIS.

### How It Works

> - Cast a Ray: A ray (straight line) is extended from the point in question to infinity in one direction
> - Count Intersections: Count how many times this ray intersects with the edges of the polygon.
>   - If the ray crosses the edges an odd number of times, the point is considered to be inside the polygon.
>   - If it crosses an edge an even number of times, the point is outside the polygon. 

ここはざっくりとまとめます。

対象の点からRay、光線＝無限大の直線、を引きます。
対象の点とポリゴンの交差の数が**奇数ならば対象の点はポリゴンの内側**、**偶数ならばポリゴンの外側**、と判断します。（図を見てもらえばわかるかと思います。）

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/554742fb-a7d1-4a91-af50-82090209f1e9.png)

> Use Cases and Limitations
> - This method works well for both convex and concave polygons.
> - Special cases like points on the polygon edge may need additional handling.

convex（凸型）のポリゴンであっても、concave（凹面）のポリゴンであっても機能しますが、対象の点がポリゴンのエッジだったりする場合は、追加の処理が必要になることもあります。
:::

## Point-in-polygon queries on shapely geometries

> There are basically two ways of conducting PIP in Shapely:
> - using a function called within() that checks if a point is within a polygon
> - using a function called contains() that checks if a polygon contains a point

> **Note**
> Even though we are discussing a Point in Polygon operation, it is also possible to check if a LineString or Polygon is inside another Polygon.

Shapelyで、PIP、Point in Polygonを実施する場合、基本的には以下の２つの方法をとります。
`within()`関数：対象の点がとあるポリゴンの中にあるかを判定する
`contains()`関数：対象のポリゴンの中にとある点があるかを判定する

もちろん、「対象の点」（`Point`）については、「対象の線」（`LineString`）でも「対象のポリゴン」（`Polygon`）でもOKです。

では、やってみましょう。

ただ、その前に、位置の関係を図にしてみました。大雑把すぎですね。
とはいえ、なんとなく理解できると思います。`point1`は`polygon`の中にありますが、`point2`はそうではありません。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/d9ce3290-0211-44aa-aeac-a0c58f091c40.png)

まずは、`within()`関数です。

```
import shapely.geometry
point1 = shapely.geometry.Point(24.952242, 60.1696017)
point2 = shapely.geometry.Point(24.976567, 60.1612500)

polygon = shapely.geometry.Polygon(
    [
        (24.950899, 60.169158),
        (24.953492, 60.169158),
        (24.953510, 60.170104),
        (24.950958, 60.169990)
    ]
)

print(point1.within(polygon))
print(point2.within(polygon))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/1199e58e-9a7a-4fcd-952c-59fbbd68f712.png)


> It seems that the first point is inside the polygon, but the second one is not.
> We can turn the logic of the look-up around: Rather than check of the point is within the polygon, we can also ask whether the polygon contains() the point:

`point1`は`polygon`の中、`point2`は`polygon`の中ではない、という結果です。期待通りですね。

では、逆のパターン、`contains()`関数です。

```
print(polygon.contains(point1))
print(polygon.contains(point2))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/62d16cc3-7d50-4a15-8898-0851346fe1b3.png)

`polygon`は`point1`を含むが`point2`は含まない、という結果です。期待通りですね。


> **Hint**: The two ways of checking the spatial relationship are complementary and yield equivalent results; contains() is inverse to within(), and vice versa.
> Then, which one should you use? Well, it depends:
> - If you have many points and just one polygon and want to find out which points are inside the polygon, you might need to iterate over the points and check one at a time if it is within() the polygon.
> - If you have many polygons and just one point and want to find out which polygon contains the point, you might need to iterate over the polygons until you find one that contains() the point.

**ヒント**：この2つの関数はお互いを補完するものです。なので、結果は同じになります。使い分けは以下です。

- 点は複数×ポリゴンは一つ、どの点がポリゴンの内側にあるかを見たい：すべての点で`within()`関数を呼ぶ（反復処理）
- ポリゴンは複数×点は一つ、どのポリゴンが点を含んでいるかを見たい：すべてのポリゴンで`contains()`関数を呼ぶ（反復処理）


後半に続きます。