Skip to content

Commit bd75976

Browse files
martiboschmartibosch
committed
feat: landscape-species and land use regression examples
1 parent 662c45e commit bd75976

22 files changed

+5266
-564
lines changed

‎README.md‎

Lines changed: 110 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -8,71 +8,136 @@
88

99
Toolkit for focal site multi-scale studies in Python.
1010

11-
Example application to compute the proportion of tree canopy around (with multiple buffer radii) weather stations in Zurich, Switzerland:
12-
1311
![stations-tree-canopy](https://github.com/martibosch/focalpy/raw/main/figures/stations-tree-canopy.png)
1412

1513
*(C) OpenStreetMap contributors, tiles style by Humanitarian OpenStreetMap Team hosted by OpenStreetMap France*
1614

1715
## Overview
1816

19-
Start by instantiating a `MultiScaleFeatureComputer` for [a given region of interest](https://github.com/martibosch/pyregeon). Then, given a list of site locations, you can compute urban features at multiple scales, i.e., based on the landscape surrounding each site for multiple buffer radii:
17+
Compute multi-scale spatial predictors:
2018

2119
```python
22-
import swisstopopy
20+
import contextily as cx
21+
import geopandas as gpd
22+
import matplotlib.pyplot as plt
23+
from sklearn import ensemble
2324

2425
import focalpy
2526

26-
# parameters
27-
region = "EPFL"
28-
crs = "epsg:2056"
29-
buffer_dists = [10, 25, 50, 100]
30-
grid_res = 200
31-
32-
# instantiate the multi-scale feature computer
33-
msfc = focalpy.MultiScaleFeatureComputer(region=region, crs=crs)
34-
35-
# generate a regular grid of points/sites within the region
36-
site_gser = msfc.generate_regular_grid_gser(grid_res, geometry_type="point")
27+
species_richness_filepath = "data/bird-richness.gpkg"
28+
buildings_gdf_filepath = "data/buildings.gpkg"
29+
tree_canopy_filepath = "data/tree-canopy.tif"
3730

38-
# get a tree canopy raster from swisstopo data
39-
dst_filepath = "tree-canopy.tif"
40-
swisstopopy.get_tree_canopy_raster(region, dst_filepath)
41-
tree_val = 1 # pixel value representing a tree in the canopy raster
31+
buffer_dists = [100, 250, 500]
4232

43-
# generate a DEM raster from swisstopo data
44-
dem_filepath = "dem.tif"
45-
swisstopopy.get_dem_raster(region, dem_filepath)
33+
species_gdf = gpd.read_file(species_richness_filepath)
34+
y_col = "n.species" # species richness
4635

47-
# compute multi-scale features
48-
49-
# building areas from OpenStreetMap buildings (via osmnx)
50-
features_df = pd.concat(
36+
fa = focalpy.FocalAnalysis(
37+
[buildings_gdf_filepath, tree_canopy_filepath],
38+
species_gdf,
39+
buffer_dists,
5140
[
52-
msfc.compute_building_features(site_gser, buffer_dists),
53-
msfc.compute_tree_features(
54-
tree_canopy_filepath, site_gser, buffer_dists, tree_val
55-
),
56-
msfc.compute_topo_features_df(
57-
dem_filepath,
58-
site_gser,
59-
buffer_dists,
60-
),
41+
"compute_vector_features",
42+
"compute_raster_features",
6143
],
62-
axis="columns",
44+
feature_col_prefixes=["building", "tree"],
45+
feature_methods_args={
46+
"compute_vector_features": [{"area": "sum"}],
47+
},
48+
feature_methods_kwargs={
49+
"compute_raster_features": {"stats": "sum"},
50+
},
6351
)
64-
features_df.head()
52+
fa.features_df.head()
53+
```
54+
55+
<table border="1" class="dataframe">
56+
<thead>
57+
<tr style="text-align: right;">
58+
<th></th>
59+
<th>building_area_sum_100</th>
60+
<th>building_area_sum_250</th>
61+
<th>building_area_sum_500</th>
62+
<th>tree_sum_100</th>
63+
<th>tree_sum_250</th>
64+
<th>tree_sum_500</th>
65+
</tr>
66+
</thead>
67+
<tbody>
68+
<tr>
69+
<th>0</th>
70+
<td>13069.227511</td>
71+
<td>60218.251616</td>
72+
<td>207368.012055</td>
73+
<td>2016.0</td>
74+
<td>14875.0</td>
75+
<td>61452.0</td>
76+
</tr>
77+
<tr>
78+
<th>1</th>
79+
<td>7439.635337</td>
80+
<td>41645.546860</td>
81+
<td>131432.855040</td>
82+
<td>1331.0</td>
83+
<td>15760.0</td>
84+
<td>84520.0</td>
85+
</tr>
86+
<tr>
87+
<th>2</th>
88+
<td>8962.495280</td>
89+
<td>54251.129360</td>
90+
<td>146157.281494</td>
91+
<td>2385.0</td>
92+
<td>16725.0</td>
93+
<td>79704.0</td>
94+
</tr>
95+
<tr>
96+
<th>3</th>
97+
<td>8001.653873</td>
98+
<td>29735.393494</td>
99+
<td>102803.559740</td>
100+
<td>2512.0</td>
101+
<td>22892.0</td>
102+
<td>95945.0</td>
103+
</tr>
104+
<tr>
105+
<th>4</th>
106+
<td>10447.531020</td>
107+
<td>39405.263870</td>
108+
<td>110922.947475</td>
109+
<td>3886.0</td>
110+
<td>19860.0</td>
111+
<td>99111.0</td>
112+
</tr>
113+
</tbody>
114+
</table>
115+
116+
```
117+
# target area (for region-wide prediction/extrapolation)
118+
study_area_filepath = "data/study-area.gpkg"
119+
grid_res = 500
120+
121+
# train a model and spatially extrapolate it
122+
model = ensemble.GradientBoostingRegressor().fit(fa.features_df, species_gdf[y_col])
123+
pred_da = fa.predict_raster(model, study_area_filepath, grid_res, pred_label=y_col)
124+
125+
# plot the field data and predicted raster
126+
fig, ax = plt.subplots()
127+
cmap = "BuGn"
128+
vmin = min(pred_da.min().item(), species_gdf[y_col].min())
129+
vmax = max(pred_da.max().item(), species_gdf[y_col].max())
130+
pred_da.plot(ax=ax, alpha=0.7, vmin=vmin, vmax=vmax, cmap=cmap)
131+
species_gdf.plot(y_col, ax=ax, edgecolor="k", vmin=vmin, vmax=vmax, cmap=cmap)
132+
ax.set_axis_off()
133+
cx.add_basemap(ax, crs=species_gdf.crs, attribution=False)
65134
```
66135

67-
| grid_cell_id | buffer_dist | building_area | tree_canopy | slope | northness | tpi |
68-
| -----------: | ----------: | ------------: | ----------: | -------: | --------: | --------: |
69-
| 0 | 10 | 313.654849 | 0.000000 | 0.020963 | 0.180932 | -0.006561 |
70-
| | 25 | 3920.685613 | 0.014260 | 0.052408 | 0.023872 | 0.036682 |
71-
| | 50 | 31365.484905 | 0.047746 | 0.070575 | -0.006432 | -0.075104 |
72-
| | 100 | 250923.879244 | 0.043386 | 0.073637 | 0.006363 | -0.716217 |
73-
| 1 | 10 | 627.309698 | 0.000000 | 0.095521 | 0.228504 | 0.080963 |
136+
<img src="https://github.com/martibosch/focalpy/raw/main/figures/pred-raster.png" alt="pred-raster" style="max-width=600px">
137+
138+
*(C) OpenStreetMap contributors, tiles style by Humanitarian OpenStreetMap Team hosted by OpenStreetMap France*
74139

75-
See the [overview notebook](https://focalpy.readthedocs.io/en/latest/overview.html) and the [API documentation](https://focalpy.readthedocs.io/en/latest/api.html) for more details on the features of focalpy.
140+
See the [user guide](https://focalpy.readthedocs.io/en/latest/user-guide/introduction.html) and the [API documentation](https://focalpy.readthedocs.io/en/latest/api.html) for more details on the features of focalpy.
76141

77142
## Installation
78143

‎docs/api.md‎

Lines changed: 1 addition & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,20 +1,6 @@
11
# API
22

3-
## Computing features
4-
53
```{eval-rst}
6-
.. autoclass:: focalpy.MultiScaleFeatureComputer
4+
.. autoclass:: focalpy.FocalAnalysis
75
:members:
86
```
9-
10-
## Standalone topographic functions
11-
12-
```{eval-rst}
13-
.. autofunction:: focalpy.topo.compute_terrain_attribute
14-
15-
.. autofunction:: focalpy.topo.northness
16-
17-
.. autofunction:: focalpy.topo.comparative_height_at_center
18-
19-
.. autofunction:: focalpy.topo.flow_accumulation_at_center
20-
```

‎docs/conf.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,9 @@ def run(self, **kwargs):
7676

7777
# do NOT execute notebooks
7878
nbsphinx_execute = "never"
79+
# do NOT add ".ipynb" to source suffixes
80+
# see https://github.com/spatialaudio/nbsphinx/issues/310
81+
source_suffix = [".md", ".rst"]
7982

8083
# theme
8184
html_theme = "pydata_sphinx_theme"

‎docs/index.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,9 @@ hidden:
44
maxdepth: 1
55
---
66
7-
overview
87
user-guide/introduction
98
api
9+
planned-features
1010
contributing
1111
changelog
1212
references

‎docs/planned-features.md‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
# Planned features
2+
3+
- Support spatial regression models from [spreg](https://github.com/pysal/spreg) and the PySAL stack {cite:p}`rey2009pysal`.
4+
- Add a method to plot *scalograms*, i.e., plotting how the computed spatial predictors respond to changes in scale, which can reveal scale thresholds that maximize landscape heterogeneity {cite:p}`pasher2013optimizing` (and therefore the variance of the spatial predictors that act as independent variables).
5+
- Implement algorithms to sample locations for field data collection based on landscape heterogeneity {cite:p}`bowler2022optimising`.
6+
- Add methods to assess the "area of applicability" {cite:p}`meyer2021predicting` of the models (based on the latent space defined by the spatial predictors) as well as the "risk of spatial extrapolation" {cite:p}`gutzwiller2023using`.

0 commit comments

Comments
 (0)