# Atelier SAGEO 2021 : Traitements spatiaux avec le plugin Python t4gpd dans le contexte d'un Jupyter Notebook

Thomas Leduc

## ► To cite this version:

Thomas Leduc. Atelier SAGEO 2021: Traitements spatiaux avec le plugin Python t4gpd dans le contexte d'un Jupyter Notebook. École thématique. La Rochelle, France. 2021. hal-03262455

HAL Id: hal-03262455

https://hal.science/hal-03262455

Submitted on 16 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Atelier SAGEO 2021

## Traitements spatiaux avec le plugin Python t4gpd dans le contexte d'un Jupyter Notebook

### Thomas Leduc

thomas.leduc@crenau.archi.fr

UMR 1563 AAU / CRENAU

<figure><img src="image_2.png" /></figure>

→ nantes

ensa

→ architecture

<figure><img src="image_3.png" /></figure>

CENTRALE
NANTES

5 mai 2021

<figure><img src="image_4.png" /></figure>
# Plan de la présentation

<figure><img src="image_5.png" /></figure>

## 1 Introduction à Pandas

<figure><img src="image_6.png" /></figure>

## 2 Introduction à Shapely

<figure><img src="image_7.png" /></figure>

## 3 Introduction à GeoPandas

<figure><img src="image_8.png" /></figure>

## 4 t4gpd

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>
# Pandas

* Bibliothèque Python de manipulation et analyse de données

* Les principales structures de données Pandas sont :

* les Series pour les données 1D

* les Frame pour les données 2D (ou plus)

* Ces structures de données sont indexées pour plus d'efficacité, couplées à des utilitaires de lecture/écriture de fichiers csv, tableurs xls(x), fichiers HDF5, bases de données SQL, etc.

* Pandas permet la concatenation et la fusion de gros volumes de données (*), l'analyse de séries temporelles (*), facilite la gestion des données manquantes (*), les tableaux croisés dynamiques (*), etc.

* Ce qui suit s'inspire largement de 10 minutes to pandas

<figure><img src="image_15.png" /></figure>
# Pandas - à propos de Series

* « Series is a one-dimensional labeled array capable of holding any data type (integers, strings, floating point numbers, Python objects, etc.). The axis labels are collectively referred to as the index. »

* Series acts similarly to a ` numpy.ndarray ` or a fixed-size ` dict `.

```python
1 from numpy import mean; from pandas import Series
2
3 s1 = Series(data=['a1', 'b2', 'c3', 'a4'])
4 print(s1[s1.str LIKE 'a'])
5
6 s2 = Series(data=range(0, 40, 10))
7 print(s2[s2 > mean(s2)])
8
9 print(s2.is_monotonic_increasing)
10
11 s3 = Series({c:'a{i}' for i,c in enumerate(['a','b','c'])})
12 print(s3)
13
14 print(s3 iloc[1], s3.loc['b'])
15
16 s4 = Series([1, 2, 3, 4, None, 2])
17 print(s4.hasnans, s4.is_unique)
```

<figure><img src="image_16.png" /></figure>
# Pandas - à propos de DFrame

* « **DataFrame** is a 2-dimensional labeled data structure with columns of potentially different types. You can think of it like a spreadsheet or SQL table, or a dict of Series objects. »

```python
1 from pandas import DFrame
2
df1 = DFrame(data=[['DE', 'Allemagne', 'Berlin'],
                 ['FR', 'France', 'Paris'], ['IT', 'Italie', 'Rome']],
             columns=['code', 'pays', 'capitale'])
3 print(df1)
4
df1.set_index(['code'], drop=True, inplace=True)
5 print(df1)
6
11 df1 reset_index inplace=True
12
14
df2 = df1[~df1.code.isin(['FR'])]
15 print(df2)
16
18
df2 reset_index drop=True, inplace=True)
17 print(df2)
19
```

```
1 %run scripts/_pandas.py

2 code pays capitale
3 DE Allemagne Berlin
4 1 FR France Paris
5 2 IT Italie Rome
6
7 code pays capitale
8 DE Allemagne Berlin
9 FR France Paris
10 IT Italie Rome
11
12 code pays capitale
13 DE Allemagne Berlin
14 1 FR France Paris
15 2 IT Italie Rome
16
17 code pays capitale
18 DE Allemagne Berlin
19 1 IT Italie Rome
```

<figure><img src="image_17.png" /></figure>
# Pandas - concatenation et fusion (jointure attributaire)

* Consulter notamment `pandas.concat(...)` et `pandas.merge(...)`

```python
1 from pandas import concat, DataFrame, merge
```

```python
2 df3 = DataFrame(data=[['Allemagne', 'Europe'],
3                     ['Italie', 'Europe']],
4                     columns=['pays', 'continent'])
```

```python
5 print(df3)
```

```python
6
```

```python
7
```

```python
8
```

```python
9 df4 = DataFrame(data=[['Inde', 'Asie']],
                     columns=['pays', 'continent'])
```

```python
10
```

```python
11
```

```python
12
```

```python
13
```

```python
14 df5 = concat([df3, df4])
```

```python
15 print(df5)
```

```python
16
```

```python
17 df5.reset_index(drop=True, inplace=True)
```

```python
18 print(df5)
```

```python
19 df6 = merge(df5, df1[[pays, 'capitale']],
20             how='left', on='pays')
```

```python
20
```

```python
21
```

```python
22
```

```python
0 pays continent
```

```python
0 Allemagne Europe
1 Italie Europe
```

```python
1 pays continent
```

```python
0 Inde Asie
```

```python
0 pays continent
```

```python
0 Allemagne Europe
1 Italie Europe
0 Inde Asie
```

```python
1 pays continent
```

```python
0 Allemagne Europe
1 Italie Europe
0 Inde Asie
```

```python
1 pays continent
```

```python
0 pays continent
```

```python
0 pays continent
```

```python
0 pays continent
```

```python
0 Allemagne Europe
1 Italie Europe
2 Inde Asie
```

```python
1 pays continent
```

```python
0 pays continent
```

```python
0 pays continent capitale
```

```python
0 Allemagne Europe Berlin
1 Italie Europe Rome
2 Inde Asie NaN
```

```python
1 pays continent capitale
```

```python
0 pays continent capitale
```

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>
# Pandas - mécanisme de vue

* Attention au mécanisme de vue¹

```python
1 from pandas import DataFrame
```

```python
2 df1 = DataFrame(data=['a', 13], ['d', 7], ['z', 4],
```

```python
3 ['a', 9], ['z', 2]], columns=['colA', 'colB'])
```

```python
4 print(df1)
```

```python
5
```

```python
6 df2 = df1[(df1.colB >= 7) & (df1.colB <= 9)]
```

```python
7 print(df2) # df2 est une vue de df1
```

```python
8
```

```python
9 # df2.colB = 10 * df2.colB
```

```python
10 # A value is trying to be set on a copy of a slice
```

```python
11 # from a DataFrame.
```

```python
12 # Try using .loc(row_indexer, col_indexer) = value
```

```python
13 # instead
```

```python
14
```

```python
15 df3 = df1.loc[df1[(df1.colB >=7) & (df1.colB <=9)].index]
```

```python
16 df3.colB = 10 * df3.colB
```

```python
17 print(df3) # df3 est une copie profonde de df1
```

```python
18
```

```python
19 df4 = df2.copy(deep=True)
```

```python
20 df4.colB = 10 * df4.colB
```

```python
21 print(df3) # df4 est une copie profonde de df2
```

```python
10 in [1: %run scripts/54_pandas.py
```

```python
colA colB
0   a   13
1   d   7
2   z   4
3   a   9
4   z   2
```

```python
colA colB
1   d   7
2   a   9
3   a   9
4   z   2
```

```python
colA colB
1   d   70
2   a   90
3   a   90
```

```python
colA colB
1   d   70
2   a   90
3   a   90
```

```python
colA colB
1   d   7
2   a   9
3   a   9
```

```python
colA colB
1   d   70
2   a   90
3   a   90
```

```python
colA colB
1   d   70
2   a   90
3   a   90
```

```python
colA colB
1   d   70
2   a   90
3   a   90
```

```python
colA colB
1   d   70
2   a   90
3   a   90
```

```python
colA colB
1   d   70
2   a   90
3   a   90
```

1. Une modification de df2 ne sera pas reportée sur df1. < ▷ > < ▷ > < ▷ > < ▷ > < ▷ >
# Plan de la présentation

1 Introduction à Pandas

2 Introduction à Shapely

3 Introduction à GeoPandas

4 t4gpd

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>
<figure><img src="image_27.png" /></figure>

Extrait de https://shapely.readthedocs.io

► Shapely is a Python package for set-theoretic analysis and manipulation of planar features using (via Python's `ctypes` module) functions from the well known and widely deployed `GEOS` library. `GEOS`, a port of the `Java Topology Suite (JTS)`, is the geometry engine of the `PostGIS` spatial extension for the `PostgreSQL` `RDBMS`. The designs of `JTS` and `GEOS` are largely guided by the `Open Geospatial Consortium's` Simple Features Access Specification 1 (*) and Shapely adheres mainly to the same set of standard classes and operations. Shapely is thereby deeply rooted in the conventions of the geographic information systems (GIS) world, but aspires to be equally useful to programmers working on non-conventional problems.

<figure><img src="image_28.png" /></figure>
# Shapely

## Déclaration de géométries — représentation discrète du monde

```python
import matplotlib.pylot as plt
from shapelygeometry import LineString, Point, Polygon

red = Point(90, 90).buffer(10)
blue = LineString([ (10, 80), (40, 40), (60, 70), (80, 40) ])
green = Polygon([ (0, 0), (30, 0), (30, 30), (0, 30) ])

plt.subplots(figsize=(0.5*8.26, 0.5*8.26))
plt plot(*red.exterior.x, 'r-', linewidth=2)
plt plot(*blue.x, 'b-.', linewidth=2)
plt plot(*green.exterior.x, 'g--', linewidth=2)
plt savefig('../../../img/6_1_shapely.pdf')
```

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>
# Shapely

## Union (efficace) de géométries

```python
1 import matplotlib.pylab as plt
2 from shapely.geometry import MultiPoint
3 from shapely.ops import unary_union
4
5 geom = list()
6 for i in range(10, 100, 10):
7     geom.append(MultiPoint([ (0, i), (i, 10), (i, 90) ] ). buffer(10))
8
9 red = unary_union(geom)
10
11 plt.subplots(figsize=(0.5*8.26, 0.5*8.26))
12 plt plot(*red.exterior.xy, 'r-', linewidth=2)
13 plt savefig('../img/6_2_shapely.pdf')
```

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>
Transformations affines de géométries

1 import matplotlib.pylot as plt
2 from shapely.affinity import rotate, translate, scale
3 from shapely geometry import Point
4 from shapely.ops import unary_union
5
6 red = Point(10, 0).buffer(10, 1)
7 blue = translate(rotate(red, 45, origin='center'), xoff=40.0, yoff=0.0)
8 green = scale(translate(blue, xoff=40.0, yoff=0.0), xfact=1, yfact=5, origin='center')
9
10 plt.subplots(figsize=(0.5*8.26, 0.5*8.26))
11 plt plot(*red. exterior xy, 'r-', linewidth=2)
12 plt plot(*blue. exterior xy, 'b-.', linewidth=2)
13 plt plot(*green. exterior xy, 'g-', linewidth=2)
14 plt savefig('../img/6_3_shapely.pdf')

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>
Shapely

Prédicats topologiques

from shapelyOfficial geometry import LineString, Point, Polygon

red = Point((30, 30)).buffer(25)
green = Point((70, 30)).buffer(25)
blue = Point((50, 60)).buffer(30)
black = Polygon(((40, 80), (60, 80), (50, 100)))
magenta = Polygon(((40, 60), (60, 60), (50, 80)))
cyan = LineString(((20, 90), (40, 100), (50, 85), (60, 100)))

print(magenta.within(blue))  # True
print(blue.contAINS(magenta))  # True
print(red.difference(blue).intersects(green))  # True
print(magenta&touches(black))  # True
print(black.overlaps(blue))  # True
print(cyan.crosses(black))  # True
print(black.union(magenta).disjoint(green))  # True

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>
* A propos du Dimensionally Extended 9-Intersection Model (DE-9IM)

$$DE9IM(a, b) = \begin{bmatrix} \dim(I(a) \cap I(b)) & \dim(I(a) \cap B(b)) & \dim(I(a) \cap E(b)) \\ \dim(B(a) \cap I(b)) & \dim(B(a) \cap B(b)) & \dim(B(a) \cap E(b)) \\ \dim(E(a) \cap I(b)) & \dim(E(a) \cap B(b)) & \dim(E(a) \cap E(b)) \end{bmatrix}$$

* où dim est la dimension de l'intersection (∩) entre l'intérieur (I), la frontière (B) ou le complémentaire (E) des géométries a et b passées en paramètre.

* cette dimension prend la valeur F si l'intersection est l'ensemble vide (Ø), 0 si l'intersection est un ensemble de points, 1 si l'intersection est un ensemble de lignes et 2 si l'intersection est un ensemble de polygones.

* cette matrice peut être sérialisée, ligne après ligne, sous forme d'un *DE-9IM string code*. Ainsi, à partir de l'exemple précédent, l'instruction Shapely : black. relate (magenta) renvoie la chaîne de caractères : FF2F01212.

<figure><img src="image_37.png" /></figure>
## Opérations géométriques

1. `green = Point((5,5)).buffer(5)`

2. `blue = Point((10, 5)).buffer(5)`

3. `_intersection = green.intersection(blue)`

4. `_union = green.union(blue)`

5. `_symmetricDifference = green.symmetricDifference(blue)`

6. `_difference = green.difference(blue)`

Donnees en entree

<figure><img src="image_38.png" /></figure>

Intersection

<figure><img src="image_39.png" /></figure>

Union

<figure><img src="image_40.png" /></figure>

Difference symétrique

<figure><img src="image_41.png" /></figure>

Difference

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>
Shapely

<figure><img src="image_44.png" /></figure>

Simplification de géométrie

```python
import matplotlib.pyplot as plt
from shapelygeometry'import LineString

black = LineString([((i//2, i-i//2) for i in range(19))])
red = black.simplify(1.0, preserveTopology=True)

plt.subplots(figsize=(0.5*8.26, 0.5*8.26))

plt plot(*black.xy, 'k-', linewidth=1)

plt plot(*red.xy, '-.', linewidth=2)

plt savefig('../img/6_7_shapely.pdf')
```

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>
● Boîte englobante (aussi appelée bounding box)

import matplotlib.pylplot as plt
1
2
3 from random import randomint
4
5 from shapely.geography import box, MultiPolygon, Point
6
7 pts = [Point(1,2).buffer(0.5), Point(7,4).buffer(0.5)]
8
9 pts = MultiPolygon pts)

print pts.bounds

#~ (minx, miny, maxx, maxy) = (0.5, 1.5, 7.5, 4.5)
5
6
7
8
9

envelop = box(*pts.bounds)

11
12
13
14
15
16
17
18

ptl.subplots(figsize=(0.5*8.26, 0.5*8.26))
ptl.xlim(0, 8)
ptl.ylim(0, 8)
ptl plot(*envelop.exterior.xy, 'r-', linewidth=2)
for pt in pts:
    pt plot(*pt.exterior.xy, 'k-', linewidth=2)
ptl.savefig('../img/6_8_shapely.pdf')

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>
# Shapely

## Triangulation via shapely.ops.triangulate (...)

```python
1 from shapely geometry import MultiPoint
2 from shapely ops import triangulate
3
4 mpts = MultiPoint([[(4,0), (7,2), (2,4), (7,4),
5                  (5,7), (4,8), (2,9), (5,9), (10,9)]])
5 tin = triangulate(mpts) #~ retourne une liste de "Polygon"
6
```

<figure><img src="image_49.png" /></figure>

## Polygonisation via shapely.ops.polygonize (...)

```python
1 from shapely geometry import MultiLineString
2 from shapely ops import polygonize
3
4 mls = MultiLineString([[(2,9),(2,4)],[(2,4),(4,8)],
5                     [(4,8),(2,9)],[(4,8),(5,9)],[(5,9),(2,9)],[(4,8),(5,7)],
6                     [(5,7),(5,9)],[(5,7),(10,9)],[(10,9),(5,9)],[(5,7),(7,4)],
7                     [(7,4),(10,9)],[(7,4),(7,2)],[(7,2),(10,9)],[(4,0),(7,2)],
8                     [(7,2),(2,4)],[(2,4),(4,0)],[(7,4),(2,4)],[(5,7),(2,4)]]
9 polygons = polygonize(mls) #~ retourne une liste de "Polygon"
```

<figure><img src="image_50.png" /></figure>
Dilatation, érosion via buffer(...)

► Point((0, 0)). buffer (distance, resolution = 16, cap_style=1, join_style=1, mitre_limit=5.0)

<table><thead><tr><th></th><th>join_style=1<br>(round)</th><th>join_style=2<br>(mitre)</th><th>join_style=3<br>(bevel)</th></tr></thead><tbody><tr><th>cap_style=1<br>(round)</th><td></td><td></td><td></td></tr><tr><th>cap_style=2<br>(flat)</th><td></td><td></td><td></td></tr><tr><th>cap_style=3<br>(square)</th><td></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /></figure>

<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /></figure>
# Shapely

## Dilatation, érosion via buffer(...)

1

from shapely Official import LineString

2

from shapely Official import CAPCANVAS, JOINCANVAS

3

4

red = LineString([ (0,0), (10,5), (13,0), (20,4)])

5

\#~ DILATATION

6

dimgrey = red.buffer(4, cap_style=CAPCANVASspec flat, join_style=JOINCANVASspec mitre)

7

\#~ EROSION

8

lightgrey = dimgrey.buffer(-2, join_style=JOINCANVASspec round)

<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>
# Shapely

## linemerge(...), project(...) et interpolate(...)

```python
1 from shapely geometry import LineString, MultiPoint
2 from shapely.ops import linemerge
3
4 black = MultiPoint([[(2, 2), (7, 5), (6, 1), (10, 3), (11, 1), (15, 3), (19, 2)])
5
6 lines = [LineString([ (0, 0), (10, 5) ]), LineString([ (10, 5), (13, 0) ]),
7         LineString([ (13, 0), (20, 4) ])]
8
9
10
11 red = []
12 for _point in black.geoms:
13     # project(...) retourne une abscisse curviligne
14     # interpolate(...) retourne un Point shapely
15     red.append(greyinterpolate(grey.project(_point)))
16
17 red = MultiPoint(red)
```

<figure><img src="image_60.png" /></figure>
* **Import-export Well-Known Text (WKT)**

► **Export au format WKT :**
*from shapely Official import LineString*
*LineString([ (0,0), (10,0), (10,10) ] ). wkt*

► **Import du format WKT :**
*from shapely import wkt*
*wktbox('LINESTRING (0 0, 10 0, 10 0)')*

<figure><img src="image_61.png" /></figure>
# Shapely

* Attention, Shapely est une bibliothèque résolument 2D. La composante en Z des géométries est "décorative".

```python
1 import matplotlib.pyplot as plt
2 from shapely.affinity import translate
3 from shapelygeometry import Polygon
4
5 red = Polygon([((0,0,0)), (10,0,0)), (10,10,0)), (0,10,0)])
6 blue = translate(red, xoff=5.0, yoff=5.0, zoff=5.0)
7 green = redcap intersection(blue)
8 print(redcap intersects(blue)) # True
9 print(green.wkt)
10 #~ POLYGON Z ((5 10 2.5, 10 10 0, 10 5 2.5, 5 5 5, 5 10 2.5))
11
12 plt.subplots(figsize=(0.5*8.26, 0.5*8.26))
13 plt plot(*red .exterior .xy, 'r: ', linewidth=2)
14 plt plot(*blue .exterior .xy, 'b: ', linewidth=2)
15 plt plot(*green .exterior .xy, 'g- ', linewidth=2)
16 plt savefig('../img/6_d_shapely.pdf')
```

<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /></figure>
# Plan de la présentation

1 Introduction à Pandas

2 Introduction à Shapely

3 Introduction à GeoPandas

4 t4gpd

<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

<figure><img src="image_69.png" /></figure>
* Extraits de https://geopandas.org

► GeoPandas is an open source project to make working with geospatial data in python easier. GeoPandas extends the datatypes used by *pandas* to allow spatial operations on geometric types. Geometric operations are performed by *shapely*. Geopandas further depends on *fiona* for file access and *descartes* and *matplotlib* for plotting.

► The goal of GeoPandas is to make working with geospatial data in python easier. It combines the capabilities of *pandas* and *shapely*, providing geospatial operations in *pandas* and a high-level interface to multiple geometries to *shapely*. GeoPandas enables you to easily do operations in python that would otherwise require a spatial database such as *PostGIS*.

<figure><img src="image_70.png" /></figure>
GeoPandas

* `geopandas.GeoFrame` est une sous-classe de `pandas.DataFrame` (de même que `geopandas.GeoSeries` est une sous-classe de `pandas.Series`)

* Chaque `GeoSeries` embarque son propre `crs`

* Un `GeoFrame` est une combinaison d'une ou plusieurs `Series` (données attributaires) et d'une ou plusieurs `GeoSeries` (données spatiales), accompagnée d'un mécanisme d'indexation

<table><thead><tr><th></th><th>colA</th><th>colB</th><th>geometry</th></tr></thead><tbody><tr><td>0</td><td>12</td><td>34</td><td>POINT (12.00000 34.00000)</td></tr><tr><td>1</td><td>56</td><td>89</td><td>POINT (56.00000 89.00000)</td></tr><tr><td>2</td><td>23</td><td>45</td><td>POINT (23.00000 45.00000)</td></tr><tr><td>3</td><td>78</td><td>90</td><td>POINT (78.00000 90.00000)</td></tr></tbody></table>

<figure><img src="image_71.png" /></figure>

<figure><img src="image_72.png" /></figure>

<figure><img src="image_73.png" /></figure>

<figure><img src="image_74.png" /></figure>

<figure><img src="image_75.png" /></figure>

<figure><img src="image_76.png" /></figure>
* Table-Oriented Programming : manipulation de données (cf. Comparison with SQL)

* `type(countries.pop_est)` #~ `pandas.core series Series`

* `type(countries_geometry)` #~ `geopandas.geoseries GeoSeries`

* `SELECT continent, pop_est FROM countries WHERE countries LIKE 'Cameroon'`

* `countries [countries.name == 'Cameroon'][[name', 'pop_est']`

* `SELECT continent, pop_est FROM countries WHERE countries LIKE 'Cameroon' OR countries LIKE 'France'`

* `countries [countries.name.isin(['Cameroon', 'France'])][['name', 'pop_est']`

* `countries [~ countries.continent. isin (['Africa', 'Asia', 'Europe', 'North America', 'South America'])][['name', 'continent']`

* `sum(countries.pop_est)`

<figure><img src="image_77.png" /></figure>
## Table-Oriented Programming : manipulation de données

* `cities` [ `cities geometry.y` == `max(cities geometry.y)` ]

* `cities` [ `abs(cities.` `geometry.y`)` == `min(abs(cities.` `geometry.y))` ]

* `cities` [ `cities.name` == `Paris` ]. `geometry.y.squeeze()`

* `countries[(countries.continent == 'Europe') & (countries.bounds.maxx < 5)].name`

* `north = cities [ cities.` `geometry.y` >= 0]`

* `south = cities [ cities.` `geometry.y` < 0]`

* `lesCapitales = south.append(north)`

* `countries.continent.unique()`

<figure><img src="image_78.png" /></figure>

<figure><img src="image_79.png" /></figure>

<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /></figure>

<figure><img src="image_82.png" /></figure>

<figure><img src="image_83.png" /></figure>
## Table-Oriented Programming : manipulation de données

* **Jointure attributaire : `pandas.merge(...)`**

```python
1 import geopandas as gpd, pandas as pd
2 from numpy import round
3
4 countries = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
5
6 continents = countries.groupby(by='continent').sum().\_
7 rename colunms={ 'pop_est': 'NHab'}
8
9 result = pd.merge(countries, continents, how='left',\_
10 left_on='continent', right_index=True, validate='m:1')
11 result['ratio'] = round(100 * result.pop_est / result.NHab, 1)
12 result.sort_values('ratio', ascending=False, inplace=True)
13
14 print(result[[ 'continent', 'name', 'ratio']].head(7))
```

### Fichier Edition Affichage Rechercher Terminal Aide

```
In [1]: %run 7_3_geopandas.py
 continent      name      ratio
23   Seven seas (open ocean)  Fr. S. Antarctic Lands  100.0
159  Antarctica      Antarctica      100.0
137  Oceania         Australia      63.2
4    North America    United States of America  57.0
29  South America    Brazil         49.5
139  Asia            China         31.4
98  Asia            India         29.2
```

```
In [2]:
```
* Table-Oriented Programming : manipulation de données

► Jointure spatiale : geopandas.sjoin (...)

1 import geopandas as gpd
2
3 countries = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
4 cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))
5
6 #~ op: binary predicate, one of {'intersects', 'contains', 'within'}
7 result = gpd.sjoin(cities, countries, how='left', op='intersects',
8                   lsuffix='city', rsuffix='country')
9 result.sort_values('pop_est', ascending=True, inplace=True)
10 result.rename(columns={'pop_est': 'NHab_country'}, inplace=True)
11
12 print(result[[ 'name.city', 'continent', 'name.country', 'NHab.country']])
head(7))

<figure><img src="image_84.png" /></figure>

<figure><img src="image_85.png" /></figure>
* Table-Oriented Programming : manipulation de données

► Jointure spatiale : geopandas.sjoin (...)

```python
import geopandas as gpd, matplotlib.pyplot as plt

countries = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))

result = gpd.sjoin(cities, countries, how='left', op='intersects',
                   lsuffix='city', rsuffix='country')
print('Nombre de capitales sans pays de rattachement : %d' %
      len(result result.name_country.isnull()))
#- Nombre de capitales sans pays de rattachement : 30

countries = countries[countries.name == 'Libya']
cities = cities[cities.name == 'Tripoli']

_ , ax = plt.subplots(figsize=(0.5*8.26, 0.5*8.26))
ax = countries.plot(ax=ax, color='grey')
cities.plot(ax=ax, color='red', marker='o')
minx, miny, maxx, maxy = cities.buffer(0.5).totalBounds
plt.axis([minx, maxx, miny, maxy])
plt savefig('../img/7_5_geopandas.pdf')
```

<figure><img src="image_86.png" /></figure>
* Table-Oriented Programming : manipulation de données

► Overlay : geopandas overlays (...)

1 green1 = Polygon([ (0,0), (10,0), (10,10), (0,10) ])

2 green2 = Polygon([ (10,10), (20,10), (20,20), (10,20) ])

3 blue1 = Polygon([ (5,5), (15,5), (15,15), (5,15) ])

4 blue2 = Polygon([ (15,15), (25,15), (25,25), (15,25), (15,15) ])

5 green = gpd. GeoFrame( [{ 'geometry': green1}, { 'geometry': green2} ] )

6 blue = gpd. GeoFrame( [{ 'geometry': blue1}, { 'geometry': blue2} ] )

7 _intersection = gpd . overlay( green , blue , how='intersection')

8 _union = gpd . overlay( green , blue , how='union')

9 _symmetricDifference = gpd . overlay( green , blue , how='symmetricDifference')

10 _difference = gpd . overlay( green , blue , how='difference')

<figure><img src="image_87.png" /></figure>
* Table-Oriented Programming : manipulation de données

► Overlay : `geopandas overlays` (...)

1. `from geopandas import GeoDataFrame, overlay`

2. `from shapely geometry import Polygon`

3.

4. `green1 = Polygon([ (0,0), (10,0), (10,10), (0,10) ])`

5. `green2 = Polygon([ (10,10), (20,10), (20,20), (10,20) ])`

6. `blue1 = Polygon([ (5,5), (15,5), (15,15), (5,15) ])`

7. `blue2 = Polygon([ (15,15), (25,15), (25,25), (15,25), (15,15) ])`

8. `green = GeoDataFrame([{'geometry': green1, 'id': 'g1'}],`

9. `{'geometry': green2, 'id': 'g2'}`)

10. `blue = GeoDataFrame([{'geometry': blue1, 'id': 'b1'}],`

11. `{'geometry': blue2, 'id': 'b2'}`)

12.

13. `print(overlay(green, blue, how='intersection'))`

Fichier Edition Affichage Rechercher Terminal Aide

In [1]: %run 7.7_geopandas.py

td_1 td_2 geometry
0 g1 b1 POLYGON ((5.00000 10.00000, 10.00000 10.00000,...)
1 g2 b2 POLYGON ((10.00000 10.00000, 10.00000 15.00000,...)
2 g2 b2 POLYGON ((15.00000 20.00000, 20.00000 20.00000,...)

In [2]:

<figure><img src="image_88.png" /></figure>
<figure><img src="image_89.png" /></figure>

## Manipulation de données : geopandas clip (...)

1. import geopandas as gpd. matplotlib. pyplot as plt

2. from shapely. geometry. import box

3. from os. import path

4. if path. exists( 'c:/Users/tleduc'): HOMEDIR = 'c:/Users/tleduc'

5. elif path. exists( '/home/tleduc'): HOMEDIR = '/home/tleduc'

6. elif path. exists( '/home/tleduc'): HOMEDIR = '/home/tleduc'

7. btopo = HOMEDIR + '/data/btopo/'

8. BDTOP0_3_0_TOUSTHEMES_SHP_LAMB93_D044_2020-06-15/BDTOPO/1_DONNEES_LIVRAISON_2020-06-0047/

9. /BDT_3_0_SHP_LAMB93_D044-ED2020-06-15/BATI'

12. roi = 355019.3, 6689077.3, 355311.4, 6689528.0

13. red = gpd. read_file(btopo + '/BATIMENT.shp ', bbox=roi)

15. roi = box( *roi)

16. grey = gpd. clip(red, roi)

18. _ , basemap = plt. subplots( figsize=(0.5*8.26, 0.5*8.26))

19. plt. axis( 'off')

20. red. plot(ax=basemap, color='red')

21. red. plot(ax=basemap, color='red')

22. grey. plot(ax=basemap, color='grey')

23. plt. savefig( '../img/7_8_geopandas.pdf ', bbox=inches='tight')

<figure><img src="image_90.png" /></figure>
## Manipulation de données : aggédation de géométries via dissolve (...)

```python
import geopandas as gpd, matplotlib.pyplot as plt

from shapely.ops import unary_union

countries = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

continents = countries[[ 'continent', 'pop_est', 'geometry']].
    dissolve(by='continent', aggfunc='sum')

_ , basemap = plt.subplots(figsize=(1.2*8.26, 0.5*8.26))
# continents.reset_index(inplace=True)
# continents.plot(ax=basemap, column='continent', cmap='tab10')
continents.plot(ax=basemap, column='pop_est', cmap='YlOrBr', legend=True)
plt savefig('../img/7_9_geopandas.pdf', bbox_inches='tight')
```

<figure><img src="image_91.png" /></figure>

* l'argument cmap de la méthode plot permet de spécifier la **palette de couleurs** choisie.
## Table-Oriented Programming : modification de données

* `countries.drop colunns=[`iso_a3`, `gdp_md_est`]`, inplace=True)
  `cities.drop(cities [(abs(cities .geometry.x) > 5)].index, inplace=True)`

* `countries.rename colunns={`pop_est`: 'NHab.', 'name': 'NomPays'}`, inplace=True)

* `from shapelygeometry import Point`
  `cities.append(( 'name': 'Nantes', 'geometry': Point((47.2173, -1.5534))), ignore_index = True )`

* `cities ['hemisphere'] = cities .geometry .apply(lambda geom: 'Nord' if (geom. y > 0) else 'Sud')`

* `result .NHab_country.dtypes #~ de type float64`
* `result .NHab_country.astype(int) #~ Cannot convert non-finite values (NA or inf) to integer`
* `result .NHab_country.fillo na(-9999, inplace = True) #~ remplacement des valeurs NaN par -9999`
* `result .NHab_country = result .NHab_country .astype(int)`

<figure><img src="image_92.png" /></figure>
## Création de `GeoDataFrame`, ajout d'un champ et modification de la référence au champ géométrique

```python
import geopandas as gpd, matplotlib.pyplot as plt
from shapely geometry import Point

gdf = gpd.GeoDataFrame([{'geometry': Point((0, i)), 'gid': i}
                       for i in range(0, 50, 10)], crs='epsg:2154')
gdf['disc'] = gdf.buffer(2, resolution=1)
gdf.set_geometry('disc', inplace=True)
gdf plot(column='gid')
plt savefig('../img/7_a_geopandas.pdf', bbox_inches='tight')
```

<figure><img src="image_93.png" /></figure>

## Modification du système de coordonnées

► **import geopandas as gpd**

    cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))
    cities = cities [cities.name == 'Paris']
    cities .crs (~ EPSG:4326 - WGS 84)
    cities = cities .to_crs('epsg:2154')
    cities .crs (~ EPSG:2154 - RGF93 / Lambert-93)

<figure><img src="image_94.png" /></figure>
## Entrées-sorties (*)

* Chargement de SHP (l'usage d'un filtre - cf. arguments **bbox** ou **mask** - peut substantiellement accélérer le chargement)

* `import geopandas as gpd`
  * `myGrid = gpd.read_file('data/ville_numérique/insee/\Filosofi2015_carreaux_200m_metropole/Filosofi2015_carreaux_200m_metropole.shp')`

* Écriture de SHP (mais aussi GeoJSON, GPKG)

* `myGrid.to_file('data/grid200.shp', driver='ESRI Shapefile')`

* Écriture de CSV (avec encodage WKT des géométries)

* `myGrid.to_csv('data/grid200.csv')`

<figure><img src="image_95.png" /></figure>
* GeoPandas en quelques liens

► https://geopandas.org/

► https://geopandas.readthedocs.io/en/latest/

► https://automating-gis-processes.github.io/

<figure><img src="image_96.png" /></figure>

<figure><img src="image_97.png" /></figure>

<figure><img src="image_98.png" /></figure>

<figure><img src="image_99.png" /></figure>

<figure><img src="image_100.png" /></figure>

<figure><img src="image_101.png" /></figure>
# Plan de la présentation

1. Introduction à Pandas

2. Introduction à Shapely

3. Introduction à GeoPandas

4. t4gpd

<figure><img src="image_102.png" /></figure>

<figure><img src="image_103.png" /></figure>

<figure><img src="image_104.png" /></figure>

<figure><img src="image_105.png" /></figure>

<figure><img src="image_106.png" /></figure>

<figure><img src="image_107.png" /></figure>
Cinq registres de forme (Lévy, 2005) :

① l'approche de la forme urbaine comme **forme du paysage urbain** - l'espace saisi visuellement dans sa tridimensionnalité et dans sa matérialité plastique,

② l'approche de la forme urbaine comme **forme sociale** - l'espace étudié dans son occupation par les divers groupes sociaux ou la distribution des activités et fonctions dans la ville,

③ l'approche de la forme urbaine comme **forme bioclimatique** - l'espace étudié dans sa dimension environnementale, microclimatique (héliothermique, écologique, etc.),

④ l'approche de la forme urbaine comme **forme des tissus** - étude des interrelations parcellaire/viaire/libre/bâti, formes de mobilité/formes de ville,

⑤ l'approche de la forme urbaine comme **forme des tracés** - forme du plan (organique/planifié, orthogonal/radioconcentrique).

<figure><img src="image_108.png" /></figure>
* L'outil t4gpd permet d'analyser les formes d'espace construit dans différents registres, de **l'analyse bioclimatique** (orientation héliothermique, vue du ciel, etc.), à **l'analyse de tracés** (orientations, distances sur un graphe, etc.), en passant par des **analyses à connotations paysagères** (visibilités, études d'alignements d'arbres, etc.) ou des **analyses de composantes des tissus urbains** (identification de rues canyons, typologie d'intersections, etc.).

* Développé au sein de AAU-CRENAU, il bénéficie autant d'un ensemble de travaux conduits depuis plusieurs décennies à l'école nationale supérieure d'architecture de Nantes, que des développements récents autour de bibliothèques telles que GeoPandas ou Shapely.

<figure><img src="image_109.png" /></figure>
# t4gpd

* t4gpd s'inscrit dans une « filiation ». Il hérite d'un ensemble de connaissances développées dans :

* le plugin SketchUp t4su

* les scripts t4qg utilisables en console Python de QGIS 2 et 3

* mais aussi Solene, GearScape, OrbisGIS, etc.

* C'est un plugin Python 3 qui se télécharge sur SourceSup et se déploie, dans le contexte d'un environnement Miniconda3 (notamment), via un : pip install t4gpd-0.2.0.tar.gz

* Pour plus de détails concernant son installation, rendez-vous sur https://github.com/crenau/t4gpd.

<figure><img src="image_110.png" /></figure>
* https://github.com/crenau/t4gpd

* https://t4gpd-docs.readthedocs.io

<figure><img src="image_111.png" /></figure>

<figure><img src="image_112.png" /></figure>

<figure><img src="image_113.png" /></figure>

<figure><img src="image_114.png" /></figure>

<figure><img src="image_115.png" /></figure>

<figure><img src="image_116.png" /></figure>

<figure><img src="image_117.png" /></figure>