Skip to content

Commit 3cdadae

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

20 files changed

+5265
-562
lines changed

‎README.md‎

Lines changed: 109 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -16,63 +16,129 @@ Example application to compute the proportion of tree canopy around (with multip
1616

1717
## Overview
1818

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:
19+
Compute multi-scale spatial predictors:
2020

2121
```python
22-
import swisstopopy
22+
import contextily as cx
23+
import geopandas as gpd
24+
import matplotlib.pyplot as plt
25+
from sklearn import ensemble
2326

2427
import focalpy
2528

26-
# parameters
27-
region = "EPFL"
28-
crs = "epsg:2056"
29-
buffer_dists = [10, 25, 50, 100]
30-
grid_res = 200
29+
species_richness_filepath = "data/bird-richness.gpkg"
30+
buildings_gdf_filepath = "data/buildings.gpkg"
31+
tree_canopy_filepath = "data/tree-canopy.tif"
3132

32-
# instantiate the multi-scale feature computer
33-
msfc = focalpy.MultiScaleFeatureComputer(region=region, crs=crs)
33+
buffer_dists = [100, 250, 500]
3434

35-
# generate a regular grid of points/sites within the region
36-
site_gser = msfc.generate_regular_grid_gser(grid_res, geometry_type="point")
35+
species_gdf = gpd.read_file(species_richness_filepath)
36+
y_col = "n.species" # species richness
3737

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
42-
43-
# generate a DEM raster from swisstopo data
44-
dem_filepath = "dem.tif"
45-
swisstopopy.get_dem_raster(region, dem_filepath)
46-
47-
# compute multi-scale features
48-
49-
# building areas from OpenStreetMap buildings (via osmnx)
50-
features_df = pd.concat(
38+
fa = focalpy.FocalAnalysis(
39+
[buildings_gdf_filepath, tree_canopy_filepath],
40+
species_gdf,
41+
buffer_dists,
5142
[
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-
),
43+
"compute_vector_features",
44+
"compute_raster_features",
6145
],
62-
axis="columns",
46+
feature_col_prefixes=["building", "tree"],
47+
feature_methods_args={
48+
"compute_vector_features": [{"area": "sum"}],
49+
},
50+
feature_methods_kwargs={
51+
"compute_raster_features": {"stats": "sum"},
52+
},
6353
)
64-
features_df.head()
54+
fa.features_df.head()
6555
```
6656

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 |
57+
<table border="1" class="dataframe">
58+
<thead>
59+
<tr style="text-align: right;">
60+
<th></th>
61+
<th>building_area_sum_100</th>
62+
<th>building_area_sum_250</th>
63+
<th>building_area_sum_500</th>
64+
<th>tree_sum_100</th>
65+
<th>tree_sum_250</th>
66+
<th>tree_sum_500</th>
67+
</tr>
68+
</thead>
69+
<tbody>
70+
<tr>
71+
<th>0</th>
72+
<td>13069.227511</td>
73+
<td>60218.251616</td>
74+
<td>207368.012055</td>
75+
<td>2016.0</td>
76+
<td>14875.0</td>
77+
<td>61452.0</td>
78+
</tr>
79+
<tr>
80+
<th>1</th>
81+
<td>7439.635337</td>
82+
<td>41645.546860</td>
83+
<td>131432.855040</td>
84+
<td>1331.0</td>
85+
<td>15760.0</td>
86+
<td>84520.0</td>
87+
</tr>
88+
<tr>
89+
<th>2</th>
90+
<td>8962.495280</td>
91+
<td>54251.129360</td>
92+
<td>146157.281494</td>
93+
<td>2385.0</td>
94+
<td>16725.0</td>
95+
<td>79704.0</td>
96+
</tr>
97+
<tr>
98+
<th>3</th>
99+
<td>8001.653873</td>
100+
<td>29735.393494</td>
101+
<td>102803.559740</td>
102+
<td>2512.0</td>
103+
<td>22892.0</td>
104+
<td>95945.0</td>
105+
</tr>
106+
<tr>
107+
<th>4</th>
108+
<td>10447.531020</td>
109+
<td>39405.263870</td>
110+
<td>110922.947475</td>
111+
<td>3886.0</td>
112+
<td>19860.0</td>
113+
<td>99111.0</td>
114+
</tr>
115+
</tbody>
116+
</table>
117+
118+
```
119+
# target area (for region-wide prediction/extrapolation)
120+
study_area_filepath = "data/study-area.gpkg"
121+
grid_res = 500
122+
123+
# train a model and spatially extrapolate it
124+
model = ensemble.GradientBoostingRegressor().fit(fa.features_df, species_gdf[y_col])
125+
pred_da = fa.predict_raster(model, study_area_filepath, grid_res, pred_label=y_col)
126+
127+
# plot the field data and predicted raster
128+
fig, ax = plt.subplots()
129+
cmap = "BuGn"
130+
vmin = min(pred_da.min().item(), species_gdf[y_col].min())
131+
vmax = max(pred_da.max().item(), species_gdf[y_col].max())
132+
pred_da.plot(ax=ax, alpha=0.7, vmin=vmin, vmax=vmax, cmap=cmap)
133+
species_gdf.plot(y_col, ax=ax, edgecolor="k", vmin=vmin, vmax=vmax, cmap=cmap)
134+
cx.add_basemap(ax, crs=species_gdf.crs, attribution=False)
135+
```
136+
137+
![pred-raster](https://github.com/martibosch/focalpy/raw/main/figures/pred-raster.png)
138+
139+
*(C) OpenStreetMap contributors, tiles style by Humanitarian OpenStreetMap Team hosted by OpenStreetMap France*
74140

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.
141+
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.
76142

77143
## Installation
78144

‎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)