Spatio-temporal analysis#
import datetime
import pylandstats as pls
We are often interested in the analysis of the temporal evolution of the configuration and composition of a particular landscape. To that end, we will use the three extracts of Veveyse district from the Swiss Land Statistics (SLS) datasets from the Swiss Federal Statistical Office for the years 1980, 1992, 2004 and 2013.
The land use/land cover (LULC) data used in this notebook ships with the docs in the data directory (see A03-swisslandstats-preprocessing.ipynb for how it is derived from the raw SLS data).
We can now use the class SpatioTemporalAnalysis, which we can instantiate with a temporally-ordered sequence of landscape snapshots.
URBAN_CLASS_VAL = 1
input_filepaths = [
"data/veveyse/LU85_4.tif",
"data/veveyse/LU97_4.tif",
"data/veveyse/LU09_4.tif",
"data/veveyse/LU18_4.tif",
]
years = ["1980", "1992", "2004", "2013"]
sta = pls.SpatioTemporalAnalysis(input_filepaths, dates=years)
Spatio-temporal data frames#
By now, SpatioTemporalAnalysis only supports class and landscape-level metrics, which can be computed by means of its methods compute_class_metrics_df and compute_landscape_metrics_df respectively. For instance, a data frame of the class-level metrics can be obtained as follows:
class_metrics_df = sta.compute_class_metrics_df()
class_metrics_df.head()
[ ] | 0% Completed | 213.61 us
[ ] | 0% Completed | 101.15 ms
[ ] | 0% Completed | 204.55 ms
[ ] | 0% Completed | 317.36 ms
[ ] | 0% Completed | 419.23 ms
[ ] | 0% Completed | 522.65 ms
[########## ] | 25% Completed | 625.00 ms
[#################### ] | 50% Completed | 727.60 ms
[#################### ] | 50% Completed | 829.79 ms
[#################### ] | 50% Completed | 937.33 ms
[#################### ] | 50% Completed | 1.05 s
[#################### ] | 50% Completed | 1.15 s
[########################################] | 100% Completed | 1.25 s
| total_area | proportion_of_landscape | number_of_patches | patch_density | largest_patch_index | total_edge | edge_density | total_core_area | core_area_proportion_of_landscape | number_of_disjunct_core_areas | ... | euclidean_nearest_neighbor_md | euclidean_nearest_neighbor_ra | euclidean_nearest_neighbor_sd | euclidean_nearest_neighbor_cv | disjunct_core_area_mn | disjunct_core_area_am | disjunct_core_area_md | disjunct_core_area_ra | disjunct_core_area_sd | disjunct_core_area_cv | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| class_val | dates | |||||||||||||||||||||
| 1 | 1980 | 702.0 | 5.225936 | 304 | 2.263083 | 0.454106 | 208400.0 | 15.514033 | 25.0 | 0.186109 | 9 | ... | 223.606798 | 1100.000000 | 169.300811 | 55.094159 | 2.777778 | 5.560000 | 1.0 | 9.0 | 2.779999 | 100.079968 |
| 1992 | 794.0 | 5.910817 | 309 | 2.300305 | 0.602993 | 227800.0 | 16.958237 | 28.0 | 0.208442 | 10 | ... | 223.606798 | 940.175425 | 156.817097 | 52.299458 | 2.800000 | 8.000000 | 1.0 | 13.0 | 3.815757 | 136.277029 | |
| 2004 | 906.0 | 6.744584 | 304 | 2.263083 | 0.699769 | 242000.0 | 18.015335 | 50.0 | 0.372218 | 12 | ... | 223.606798 | 1100.000000 | 162.962720 | 54.136494 | 4.166667 | 10.720000 | 1.5 | 16.0 | 5.225472 | 125.411323 | |
| 2013 | 1041.0 | 7.749572 | 287 | 2.136529 | 1.421872 | 256600.0 | 19.102211 | 90.0 | 0.669992 | 20 | ... | 223.606798 | 1100.000000 | 168.551251 | 54.504167 | 4.500000 | 20.400000 | 1.0 | 31.0 | 8.458723 | 187.971629 | |
| 2 | 1980 | 8351.0 | 62.167796 | 32 | 0.238219 | 52.170029 | 615000.0 | 45.782774 | 4076.0 | 30.343185 | 107 | ... | 200.000000 | 470.820393 | 85.842455 | 36.252869 | 38.093458 | 740.707066 | 3.0 | 1352.0 | 163.600067 | 429.470245 |
5 rows × 66 columns
Again, we can use the operations of any pandas data frame. For instance, we can get all the metrics for the urban class (class_val of 1) in 1992:
class_metrics_df.loc[(1, "1992")]
total_area 794.000000
proportion_of_landscape 5.910817
number_of_patches 309.000000
patch_density 2.300305
largest_patch_index 0.602993
...
disjunct_core_area_am 8.000000
disjunct_core_area_md 1.000000
disjunct_core_area_ra 13.000000
disjunct_core_area_sd 3.815757
disjunct_core_area_cv 136.277029
Name: (1, 1992), Length: 66, dtype: float64
Similarly, the data frame of landscape metrics can be obtained as follows:
sta.compute_landscape_metrics_df()
[ ] | 0% Completed | 146.65 us
[ ] | 0% Completed | 113.92 ms
[ ] | 0% Completed | 221.51 ms
[ ] | 0% Completed | 327.15 ms
[#################### ] | 50% Completed | 429.28 ms
[#################### ] | 50% Completed | 540.57 ms
[#################### ] | 50% Completed | 653.02 ms
[#################### ] | 50% Completed | 755.64 ms
[########################################] | 100% Completed | 860.04 ms
| total_area | number_of_patches | patch_density | largest_patch_index | total_edge | edge_density | total_core_area | number_of_disjunct_core_areas | landscape_shape_index | effective_mesh_size | ... | euclidean_nearest_neighbor_md | euclidean_nearest_neighbor_ra | euclidean_nearest_neighbor_sd | euclidean_nearest_neighbor_cv | disjunct_core_area_mn | disjunct_core_area_am | disjunct_core_area_md | disjunct_core_area_ra | disjunct_core_area_sd | disjunct_core_area_cv | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| dates | |||||||||||||||||||||
| 1980 | 13433.0 | 666 | 4.957939 | 52.170029 | 677000.0 | 50.398273 | 5667.0 | 200 | 16.711207 | 4074.760143 | ... | 223.606798 | 3124.154028 | 216.625926 | 66.540480 | 28.335000 | 533.272195 | 3.0 | 1352.0 | 124.855928 | 440.642061 |
| 1992 | 13433.0 | 662 | 4.928162 | 48.596739 | 686000.0 | 51.068265 | 5572.0 | 209 | 16.905172 | 3721.028140 | ... | 223.606798 | 3124.154028 | 222.286812 | 67.856637 | 26.660287 | 534.324217 | 3.0 | 1444.0 | 121.661679 | 456.340467 |
| 2004 | 13433.0 | 652 | 4.853718 | 47.829971 | 691400.0 | 51.470260 | 5476.0 | 219 | 17.021552 | 3711.910593 | ... | 223.606798 | 3124.154028 | 231.885409 | 69.562278 | 25.004566 | 777.260594 | 3.0 | 2049.0 | 143.674210 | 574.591890 |
| 2013 | 13433.0 | 636 | 4.734609 | 46.862205 | 703300.0 | 52.356138 | 5373.0 | 233 | 17.278017 | 3593.476141 | ... | 223.606798 | 3124.154028 | 236.550901 | 70.451170 | 23.060086 | 509.698816 | 3.0 | 1555.0 | 111.314332 | 482.714300 |
4 rows × 71 columns
Customizing your spatio-temporal analysis#
As within the Landscape analysis, we can also choose to compute a subset of metrics by passing them to the metrics keyword argument of the compute_class_metrics_df and compute_landscape_metrics_df methods, as in:
metrics = ["proportion_of_landscape", "edge_density", "fractal_dimension_am"]
sta.compute_class_metrics_df(metrics=metrics)
[ ] | 0% Completed | 160.03 us
[########################################] | 100% Completed | 101.03 ms
| proportion_of_landscape | edge_density | fractal_dimension_am | ||
|---|---|---|---|---|
| class_val | dates | |||
| 1 | 1980 | 5.225936 | 15.514033 | 1.073800 |
| 1992 | 5.910817 | 16.958237 | 1.083145 | |
| 2004 | 6.744584 | 18.015335 | 1.089937 | |
| 2013 | 7.749572 | 19.102211 | 1.113777 | |
| 2 | 1980 | 62.167796 | 45.782774 | 1.287254 |
| 1992 | 60.932033 | 46.281545 | 1.280678 | |
| 2004 | 59.941934 | 46.936649 | 1.281927 | |
| 2013 | 58.862503 | 47.606640 | 1.284123 | |
| 3 | 1980 | 29.531750 | 32.643490 | 1.206868 |
| 1992 | 30.313407 | 32.717933 | 1.218381 | |
| 2004 | 30.611181 | 32.330827 | 1.225019 | |
| 2013 | 30.715402 | 32.353160 | 1.224259 | |
| 4 | 1980 | 3.074518 | 6.856250 | 1.064272 |
| 1992 | 2.843743 | 6.178813 | 1.061158 | |
| 2004 | 2.702300 | 5.657709 | 1.055708 | |
| 2013 | 2.672523 | 5.650264 | 1.058115 |
At the class-level, we can choose to compute the metrics only for a subset of classes through the classes argument. We can simoultaneously choose a subset of metrics as well as a subset of classes by specifying both the metrics and classes arguments. For instance, we can choose to only compute the above metrics and only for the urban class (value of 1):
sta.compute_class_metrics_df(metrics=metrics, classes=[URBAN_CLASS_VAL])
[ ] | 0% Completed | 168.54 us
[########################################] | 100% Completed | 101.52 ms
| proportion_of_landscape | edge_density | fractal_dimension_am | ||
|---|---|---|---|---|
| class_val | dates | |||
| 1 | 1980 | 5.225936 | 15.514033 | 1.073800 |
| 1992 | 5.910817 | 16.958237 | 1.083145 | |
| 2004 | 6.744584 | 18.015335 | 1.089937 | |
| 2013 | 7.749572 | 19.102211 | 1.113777 |
In both the compute_class_metrics_df and compute_landscape_metrics_df methods, we can also customize how some metrics are computed through the metrics_kwargs argument:
metrics_kwargs = {
"proportion_of_landscape": {"percent": False},
"edge_density": {"count_boundary": True},
}
sta.compute_class_metrics_df(
metrics=metrics, classes=[URBAN_CLASS_VAL], metrics_kwargs=metrics_kwargs
)
[ ] | 0% Completed | 160.46 us
[########################################] | 100% Completed | 101.13 ms
| proportion_of_landscape | edge_density | fractal_dimension_am | ||
|---|---|---|---|---|
| class_val | dates | |||
| 1 | 1980 | 0.052259 | 15.692697 | 1.073800 |
| 1992 | 0.059108 | 17.151790 | 1.083145 | |
| 2004 | 0.067446 | 18.223777 | 1.089937 | |
| 2013 | 0.077496 | 19.325542 | 1.113777 |
On the other hand, the dates keyword argument might also be provided as string or datetime objects, e.g.:
dates = [datetime.date(int(year), 1, 1) for year in years]
sta = pls.SpatioTemporalAnalysis(input_filepaths, dates=dates)
sta.compute_landscape_metrics_df()
[ ] | 0% Completed | 172.04 us
[ ] | 0% Completed | 101.10 ms
[ ] | 0% Completed | 222.15 ms
[ ] | 0% Completed | 326.48 ms
[ ] | 0% Completed | 431.53 ms
[#################### ] | 50% Completed | 532.34 ms
[#################### ] | 50% Completed | 644.84 ms
[#################### ] | 50% Completed | 762.27 ms
[#################### ] | 50% Completed | 864.12 ms
[#################### ] | 50% Completed | 968.94 ms
[########################################] | 100% Completed | 1.07 s
| total_area | number_of_patches | patch_density | largest_patch_index | total_edge | edge_density | total_core_area | number_of_disjunct_core_areas | landscape_shape_index | effective_mesh_size | ... | euclidean_nearest_neighbor_md | euclidean_nearest_neighbor_ra | euclidean_nearest_neighbor_sd | euclidean_nearest_neighbor_cv | disjunct_core_area_mn | disjunct_core_area_am | disjunct_core_area_md | disjunct_core_area_ra | disjunct_core_area_sd | disjunct_core_area_cv | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| dates | |||||||||||||||||||||
| 1980-01-01 | 13433.0 | 666 | 4.957939 | 52.170029 | 677000.0 | 50.398273 | 5667.0 | 200 | 16.711207 | 4074.760143 | ... | 223.606798 | 3124.154028 | 216.625926 | 66.540480 | 28.335000 | 533.272195 | 3.0 | 1352.0 | 124.855928 | 440.642061 |
| 1992-01-01 | 13433.0 | 662 | 4.928162 | 48.596739 | 686000.0 | 51.068265 | 5572.0 | 209 | 16.905172 | 3721.028140 | ... | 223.606798 | 3124.154028 | 222.286812 | 67.856637 | 26.660287 | 534.324217 | 3.0 | 1444.0 | 121.661679 | 456.340467 |
| 2004-01-01 | 13433.0 | 652 | 4.853718 | 47.829971 | 691400.0 | 51.470260 | 5476.0 | 219 | 17.021552 | 3711.910593 | ... | 223.606798 | 3124.154028 | 231.885409 | 69.562278 | 25.004566 | 777.260594 | 3.0 | 2049.0 | 143.674210 | 574.591890 |
| 2013-01-01 | 13433.0 | 636 | 4.734609 | 46.862205 | 703300.0 | 52.356138 | 5373.0 | 233 | 17.278017 | 3593.476141 | ... | 223.606798 | 3124.154028 | 236.550901 | 70.451170 | 23.060086 | 509.698816 | 3.0 | 1555.0 | 111.314332 | 482.714300 |
4 rows × 71 columns
Plots#
One of the most important features of SpatioTemporalAnalysis is to plot the evolution of the metrics. We can plot the proportion o landscape occupied by the urban class (class_val of 1) as in:
sta.plot_metric("proportion_of_landscape", class_val=URBAN_CLASS_VAL)
[ ] | 0% Completed | 146.60 us
[########################################] | 100% Completed | 101.48 ms
<Axes: ylabel='PLAND'>
If we want to plot the evolution of a metric at the landscape level, we can do so by using the same plot_metric method, but without setting the class_val argument. Note however that we cannot compute the proportion_of_landscape at the landscape level (we could but it makes no sense, the landscape always occupies 100% of the landscape). Similarly, some metrics such as shannon_diversity_index cannot be computed at the class level but only at the landscape level. See the documentation of each metric for more details.
Let’s then plot the area-weighted fractal dimension, that is fractal_dimension_am, at both the class and landscape level:
ax = sta.plot_metric(
"fractal_dimension_am",
class_val=URBAN_CLASS_VAL,
plot_kwargs={"label": "class level (urban)"},
)
_ = sta.plot_metric(
"fractal_dimension_am", ax=ax, plot_kwargs={"label": "landscape level"}
)
ax.legend()
[ ] | 0% Completed | 131.84 us
[########################################] | 100% Completed | 100.98 ms
[ ] | 0% Completed | 162.89 us
[########################################] | 100% Completed | 102.30 ms
<matplotlib.legend.Legend at 0x736d65101010>
Note that we can pass keyword arguments to matplotlib’s plot method through the plt_kws argument of plot_metric. See the documentation of SpatioTemporalAnalysis.plot_metric for more details.