Species distributions¶

Preparing occurrence records and environmental predictors for species distribution models, fitting the models and projecting them: the ls.sdm module. Every step reproduces the R package the field relies on for it, and each function's notes name that package.

import landschaft as ls

occ = ls.load_table(ls.sdm.example_data("bradypus"))
env = ls.load_raster(ls.sdm.example_data("bioclim"))
data = ls.sdm.prepare(occ, env)          # clean, one per cell, predictors, background, spatial folds
res = ls.sdm.fit(data, models)           # cross-validated scikit-learn models (or ls.sdm.MaxNet)
suit = ls.sdm.predict(res, env)          # ensemble and its uncertainty

The whole preparation¶

prepare ¶

prepare(
    occ,
    env,
    *,
    species: str | None = None,
    clean=True,
    one_per_cell: bool = True,
    thin=None,
    area=None,
    predictors=10.0,
    background=10000,
    folds=5,
    weights: str | None = "balanced",
    seed: int = 1,
    x: str | None = None,
    y: str | None = None,
    crs=None,
)

Prepare occurrence records and predictors for a species distribution model, in one call.

The steps, in order, each one a function of its own that can be called alone (and each argument here takes False to skip the step, a number for its main setting, a dict of that function's arguments, or a function of yours):

clean – clean_occurrences() (drop records at centroids, capitals, institutions, in the sea, outliers...);
one_per_cell – one_per_cell();
thin – thin() (a distance in metres; off by default);
area – the accessible area: accessible_area() (a buffer distance in metres, its arguments, or your polygons; default the whole extent of env);
predictors – select_predictors() in the area (a VIF threshold, default 10);
background – background() points in the area (a number, default 10 000);
dataset() – predictor values at presences and background;
folds – cross-validation folds: a number k of spatial blocks (spatial_blocks(), default 5) or a dict with method one of "blocks", "clusters", "quadrants", "checkerboard", "random", "jackknife" and the method's arguments; blocks take size="auto" for the predictors' range of autocorrelation (autocorrelation_range()).

Every random step uses seed, so the same call gives the same data. data.report counts what each step removed, data.settings holds every choice (for the methods of a paper) and data.redo(...) repeats the preparation with some settings changed.

Parameters:

Name	Type	Description	Default
`occ`		Occurrence records: GeoDataFrame of points or a table with coordinate columns (e.g. `download_gbif()`'s or `read_occurrences()`').	required
`env`		Predictors: a raster stack (band, y, x), e.g. `ls.load_raster(files, stack=True)`.	required
`species`	`str \| None`	Keep only this species (the records' `species` column).	`None`
`clean`		The steps (see above).	`True`
`one_per_cell`		The steps (see above).	`True`
`thin`		The steps (see above).	`True`
`area`		The steps (see above).	`True`
`predictors`		The steps (see above).	`True`
`background`		The steps (see above).	`True`
`folds`		The steps (see above).	`True`
`weights`	`str \| None`	Sample weights of `dataset()` (`"balanced"` by default).	`'balanced'`
`seed`	`int`	Seed of every random step.	`1`
`x`	`str \| None`	Coordinate columns and CRS of a table of records.	`None`
`y`	`str \| None`	Coordinate columns and CRS of a table of records.	`None`
`crs`	`str \| None`	Coordinate columns and CRS of a table of records.	`None`

Returns:

Type	Description
`SDMData`

Examples:

>>> import landschaft as ls
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> data = ls.sdm.prepare(occ, env, clean=False, background=2000, area=1_000_000)
>>> data.report[["step", "presences", "removed"]]
               step  presences  removed
0             input        116        0
1      one per cell         94       22
2              area         94        0
3        predictors         94        0
4        background         94        0
5  predictor values         94        0
6             folds         94        0
>>> data.predictors
['bio5', 'bio7', 'bio16', 'bio17']
>>> again = data.redo(predictors=False)       # the same, keeping every predictor
>>> len(again.predictors)
9

Occurrence records¶

download_gbif ¶

download_gbif(
    species,
    *,
    country=None,
    area=None,
    year=None,
    basis_of_record=None,
    max_uncertainty: float | None = None,
    limit: int | None = None,
    has_geospatial_issue: bool = False,
    n_workers: int = 4,
) -> DataFrame

Occurrence records of a species from GBIF, with coordinates.

Uses GBIF's occurrence search API: no account is needed, up to 100 000 records per query (narrow it with country, area or year). For larger queries, and for the DOI a publication must cite, request a download on gbif.org and read it with read_occurrences().

Parameters:

Name	Type	Description	Default
`species`	`str or int`	Scientific name (matched to GBIF's backbone taxonomy; a synonym gives the records of the accepted name) or a GBIF taxon key.	required
`country`	`str or list of str`	ISO 3166-1 alpha-2 country codes (`"BR"`).	`None`
`area`	`GeoDataFrame or geometry`	Only the records in its bounding box (cut to the exact shape with `accessible_area()` or a spatial join afterwards).	`None`
`year`	`(int, str or tuple)`	A year, or a range `(1990, 2020)` / `"1990,2020"`.	`None`
`basis_of_record`	`str or list of str`	E.g. `"PRESERVED_SPECIMEN"`, `"HUMAN_OBSERVATION"`.	`None`
`max_uncertainty`	`float`	Drop the records whose coordinate uncertainty is above this (metres); records without one are kept.	`None`
`limit`	`int`	At most this many records (default: all, up to 100 000).	`None`
`has_geospatial_issue`	`bool`	False (default): leave out the records GBIF flags with geospatial issues (zero coordinates, country mismatch...).	`False`
`n_workers`	`int`	Pages requested at the same time.	`4`

Returns:

Type	Description
`DataFrame`	One row per record: `gbifID`, `species`, `scientificName`, `decimalLongitude`, `decimalLatitude`, `coordinateUncertaintyInMeters`, `countryCode`, `year`, `basisOfRecord`, `institutionCode`, `datasetKey`, `issue` and other Darwin Core terms. `attrs["gbif"]` keeps the query and the taxon match.

Examples:

>>> import landschaft as ls
>>> occ = ls.sdm.download_gbif("Tapirus terrestris", country="BR", limit=500)
>>> occ = ls.sdm.clean_occurrences(occ)

clean_occurrences ¶

clean_occurrences(
    occ,
    *,
    species: str | None = None,
    equal=True,
    zeros=0.5,
    capitals=10000,
    centroids=1000,
    seas=True,
    urban=False,
    countries=False,
    outliers=True,
    gbif=1000,
    institutions=100,
    duplicates=False,
    aohi=False,
    action: str = "drop",
    x: str | None = None,
    y: str | None = None,
    crs=None,
)

Flag or drop occurrence records with suspicious coordinates.

The tests of CoordinateCleaner (Zizka et al. 2019, Methods in Ecology and Evolution 10: 744-751), the standard for GBIF data: records at country or province centroids and capitals (coordinates filled in from the name of a place), at biodiversity institutions (the address of a museum or herbarium) or GBIF's headquarters, at 0/0 or with equal longitude and latitude, in the sea, in towns, outside the country they claim, far from the other records of their species, or duplicated. Records without valid coordinates are always flagged (invalid).

Each test takes False (skip it), True (its defaults), a number (its main option: a radius) or a dict of its options:

============== ======================================================================== test options (defaults) ============== ======================================================================== equal (none) zeros buffer 0.5 (degrees around 0/0) capitals buffer 10 000 m, geod True, ref None, verify False centroids buffer 1 000 m, detail "both" / "country" / "provinces", geod True, ref None, verify False seas scale 50 (Natural Earth land, 110 / 50 / 10), ref None urban ref None (Natural Earth urban areas) countries column (the records' country codes, required), ref None, ref_col None (iso_a2_eh for two-letter codes, else iso_a3_eh) outliers method "quantile" / "mad" / "distance", mltpl 5, tdi 1000 km, size 7 (species with more records are tested), min_occs 7 (distinct locations), thinning False, thinning_res 0.5 gbif buffer 1 000 m, geod True, verify False institutions buffer 100 m, geod False, ref None, verify False, verify_mltpl 10 duplicates additions None (more columns that must also be equal) aohi buffer 10 000 m, geod True, taxa (Aves, Insecta, Mammalia, Plantae) ============== ========================================================================

ref replaces the reference data (a table or GeoDataFrame of points, or of polygons for seas / urban / countries); passing it in place of the options also works (seas=my_land). The reference data are downloaded on first use and cached (see reference_data()).

The results are CoordinateCleaner's (clean_coordinates with the same options), except where it errs: the zeros radius is in degrees as documented (terra >= 1.8 made it 0.5 m), centroids "provinces" selects the provinces, circles that cross the antimeridian or hold a pole are drawn correctly, a record on a border between two countries is not misaligned, the country codes are matched against iso_a3_eh / iso_a2_eh (with iso_a3 France and Norway are "-99" and all their records failed), species too small for the outlier test are not flagged, the outlier test also runs on species with 10 000 records or more, and institutions' verify works. Invalid coordinates are flagged instead of stopping.

Parameters:

Name	Type	Description	Default
`occ`		Occurrence records: GeoDataFrame of points, or a table with coordinate columns (`decimalLongitude` / `decimalLatitude`, `lon` / `lat`...).	required
`species`	`str \| None`	Column with the species, for the outlier and duplicate tests and `verify` (default: `species`, `scientificName` or `label` when present; else all records are one species).	`None`
`equal`		The tests and their options (see above).	`True`
`zeros`		The tests and their options (see above).	`True`
`capitals`		The tests and their options (see above).	`True`
`centroids`		The tests and their options (see above).	`True`
`seas`		The tests and their options (see above).	`True`
`urban`		The tests and their options (see above).	`True`
`countries`		The tests and their options (see above).	`True`
`outliers`		The tests and their options (see above).	`True`
`gbif`		The tests and their options (see above).	`True`
`institutions`		The tests and their options (see above).	`True`
`duplicates`		The tests and their options (see above).	`True`
`aohi`		The tests and their options (see above).	`True`
`action`	`('drop', 'flag')`	`"drop"`: return the records that pass every test; `"flag"`: return all of them with a `flag_<test>` column per test run (True = suspicious) and `flagged` (any).	`"drop"`
`x`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`
`y`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`
`crs`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`

Returns:

Type	Description
`GeoDataFrame or DataFrame`	Same type and index as `occ`. `attrs["clean"]` counts the records each test flagged.

Examples:

>>> import landschaft as ls
>>> occ = ls.load_table(ls.sdm.example_data("acaule"))
>>> flags = ls.sdm.clean_occurrences(occ, action="flag", seas=False)
>>> int(flags["flag_invalid"].sum()), len(ls.sdm.clean_occurrences(occ, seas=False))
(284, 1050)

reference_data `cached` ¶

reference_data(name: str, *, refresh: bool = False)

A geographic reference data set for cleaning occurrence records.

Downloaded once (about 5 MB) from the CoordinateCleaner archive on CRAN and cached (see cache_dir()).

Parameters:

Name	Type	Description	Default
`name`	`('countryref', 'institutions', 'aohi', 'buffland', 'buffsea', 'land110', 'land50', 'land10', 'urban50', 'countries50')`	`"countryref"`: centroids of the countries and provinces and the country capitals (5305 rows: `iso3`, `iso2`, `name`, `type`, `centroid.lon` / `.lat`, `capital.lon` / `.lat`...); `"institutions"`: 11 601 biodiversity institutions (herbaria, museums, zoos, botanical gardens, universities) with coordinates; `"aohi"`: artificial hotspot occurrence inventory (Park et al. 2023); `"buffland"` / `"buffsea"`: the land / the sea of Natural Earth (1:50 million) widened by one degree; `"land110"`, `"land50"`, `"land10"`: Natural Earth land at 1:110, 1:50 and 1:10 million (4.1.0 / 4.1.0 / 5.1.1); `"urban50"`: urban areas (4.1.0); `"countries50"`: countries (5.1.1, lowercase field names).	`"countryref"`
`refresh`	`bool`	Rebuild the cached copy from the archive.	`False`

Returns:

Type	Description
`DataFrame or GeoDataFrame`

Examples:

>>> import landschaft as ls
>>> capitals = ls.sdm.reference_data("countryref").dropna(subset=["capital"])
>>> row = capitals.iloc[0]
>>> row["name"], row["capital"], float(row["capital.lon"]), float(row["capital.lat"])
('Afghanistan', 'Kabul', 69.18, 34.52)

one_per_cell ¶

one_per_cell(
    occ,
    raster,
    *,
    x: str | None = None,
    y: str | None = None,
    drop_missing: bool = False,
)

Keep one occurrence record per cell of the predictor grid.

Several records in the same cell carry the same environmental values: kept, they weigh that cell more than the others (pseudo-replication) and put the same values on both sides of a validation split. The first record of each cell is kept, in the order of the table, as occ[!duplicated(terra::cellFromXY(r, xy)), ] in R. Records without coordinates or outside the grid are dropped.

Parameters:

Name	Type	Description	Default
`occ`		Occurrence records: a GeoDataFrame of points, or a table with coordinate columns (`lon` / `lat`, `decimalLongitude` / `decimalLatitude`, `x` / `y`...).	required
`raster`		The predictors (or any raster on their grid): a DataArray, one band or a stack. The records are reprojected to its CRS to find their cells.	required
`x`	`str \| None`	Names of the coordinate columns of a table (default: found by name).	`None`
`y`	`str \| None`	Names of the coordinate columns of a table (default: found by name).	`None`
`drop_missing`	`bool`	Also drop the records whose cell has no value (NoData) in some band.	`False`

Returns:

Type	Description
`GeoDataFrame or DataFrame`	The records kept, with their original index.

Examples:

>>> import landschaft as ls
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> len(occ), len(ls.sdm.one_per_cell(occ, env))
(116, 94)

thin ¶

thin(
    occ,
    distance: float,
    *,
    reps: int = 10,
    seed: int | None = None,
    x: str | None = None,
    y: str | None = None,
    crs=None,
    return_all: bool = False,
)

Thin occurrence records so that no two kept records are closer than distance.

Records cluster where people sample (near roads, towns, universities), and a model fitted to them learns that bias along with the species' preferences. Thinning removes records until the remaining ones are at least distance apart while keeping as many as possible: the randomised heuristic of spThin (Aiello-Lammens et al. 2015, Ecography 38: 541-545), which repeatedly drops a record with the most neighbours (a random one among ties) and keeps the best of reps repetitions. The same records as spThin::thin for the same seed; distances on a sphere of radius 6378.388 km, as spThin's fields::rdist.earth.

Parameters:

Name	Type	Description	Default
`occ`		Occurrence records: GeoDataFrame of points, or a table with coordinate columns.	required
`distance`	`float`	Minimum distance between kept records, in metres.	required
`reps`	`int`	Repetitions of the heuristic; the one keeping the most records is returned (the first among equals). More repetitions keep more records, rarely many more beyond 10-100.	`10`
`seed`	`int \| None`	Seed of R's generator (`set.seed`); default random (recorded in `attrs["thin"]`).	`None`
`x`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`
`y`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`
`crs`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`
`return_all`	`bool`	Return the records kept by every repetition (a list, largest first) instead of the best one.	`False`

Returns:

Type	Description
`GeoDataFrame or DataFrame`	The records kept, with their original index; `attrs["thin"]` holds the seed and the number of records each repetition kept. Records without coordinates are always kept (they have no neighbours): drop them first with `clean_occurrences()`.

Examples:

>>> import landschaft as ls
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> kept = ls.sdm.thin(occ, 50_000, reps=20, seed=1)
>>> len(occ), len(kept)
(116, 85)

Study area¶

accessible_area ¶

accessible_area(
    occ,
    *,
    method: str = "buffer",
    distance: float | None = None,
    regions=None,
    within=None,
    resolution: int = 64,
    x: str | None = None,
    y: str | None = None,
    crs=None,
)

The area the species could have reached: where its background is drawn.

A model learns to tell the conditions where the species was found from those available to it. Drawing the background from areas the species never had access to (beyond a mountain range, on another continent) makes the model learn geography instead of ecology and inflates its scores (Barve et al. 2011, Ecological Modelling 222: 1810-1819). The accessible area bounds the background sampling, the predictor selection and the evaluation.

Parameters:

Name	Type	Description	Default
`occ`		Occurrence records: GeoDataFrame of points, or a table with coordinate columns.	required
`method`	`('buffer', 'convex_hull', 'regions')`	`"buffer"` (default): the land within `distance` metres of a record (geodesic circles on the WGS 84 ellipsoid, merged); `"convex_hull"`: the smallest convex polygon around the records, widened by `distance` metres if given; `"regions"`: the polygons of `regions` (ecoregions, biomes, watersheds) that hold at least one record, widened by `distance` if given.	`"buffer"`
`distance`	`float \| None`	Radius in metres. Required for `"buffer"`.	`None`
`regions`		GeoDataFrame of polygons, for `method="regions"`.	`None`
`within`		Polygons the result is cut to (a country, a continent, the land).	`None`
`resolution`	`int`	Vertices of each circle (`"buffer"`) or of the widening.	`64`
`x`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`
`y`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`
`crs`	`str \| None`	Coordinate columns and CRS of a table (default: found by name; longitude / latitude in EPSG:4326).	`None`

Returns:

Type	Description
`GeoDataFrame`	One (multi)polygon in longitude / latitude (EPSG:4326), with the method and distance as columns.

Examples:

>>> import landschaft as ls
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> m = ls.sdm.accessible_area(occ, distance=300_000)
>>> m.crs.to_epsg(), round(float(m.to_crs("EPSG:6933").area.iloc[0]) / 1e12, 1)    # million km2
(4326, 7.5)

Predictors¶

select_predictors ¶

select_predictors(
    env,
    *,
    vif: float | None = 10.0,
    cor: float | None = None,
    keep=None,
    area=None,
    method: str = "pearson",
    size: int | None = _SAMPLE_SIZE,
    seed: int | None = None,
)

Drop predictors that the others already explain (collinearity).

Climate layers are strongly correlated (the 19 bioclimatic variables carry far fewer than 19 independent signals); collinear predictors make a model's coefficients unstable and its projections to new conditions unreliable (Dormann et al. 2013, Ecography 36: 27-46). Two filters of the R package usdm, with the same results:

cor (usdm's vifcor): while the most correlated pair has |r| >= cor, drop the one of the two with the larger VIF;
vif (usdm's vifstep): while the largest VIF is >= vif, drop that variable.

With both, the correlation filter runs first. Variables in keep are never dropped. Where usdm errs, the correct form: the loops stop when one variable (or only kept ones) remains instead of failing, names with spaces work, rasters of 7000-10 000 cells are sampled at random (terra took their first cells) and tables are sampled from 7000 rows on in both filters (usdm's vifstep already from 6000).

Parameters:

Name	Type	Description	Default
`env`		Raster stack (band, y, x) or a table with one column per predictor.	required
`vif`	`float \| None`	VIF threshold (usdm's default 10; 5 is stricter); None: no VIF filter.	`10.0`
`cor`	`float \| None`	Correlation threshold (e.g. 0.7); None: no correlation filter.	`None`
`keep`		Names of predictors to keep whatever their collinearity.	`None`
`area`		Polygons the raster is cut to first: collinearity in the accessible area, not in the whole raster.	`None`
`method`	`('pearson', 'spearman', 'kendall')`	Correlation coefficient of the `cor` filter.	`"pearson"`
`size`	`int \| None`	Cells sampled from rasters of 7000 cells or more (usdm's 5000; None: all cells).	`_SAMPLE_SIZE`
`seed`	`int \| None`	Seed of R's generator for the sample (default random).	`None`

Returns:

Type	Description
`DataArray or DataFrame`	The predictors kept (the same type as `env`), with `attrs["selection"]`: the variables excluded in order, the VIF of the kept ones and the seed.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> kept = ls.sdm.select_predictors(env, vif=10, seed=1)
>>> [str(b) for b in kept.band.values]
['bio5', 'bio7', 'bio8', 'bio16', 'bio17']
>>> kept.attrs["selection"]["excluded"]
['bio1', 'bio6', 'bio9', 'bio12']

vif ¶

vif(
    env, *, area=None, size: int | None = _SAMPLE_SIZE, seed: int | None = None
) -> DataFrame

Variance inflation factor of each predictor: how well the others predict it.

VIF = 1 / (1 - R²) of each variable regressed on all the others: 1 for an independent variable, above 5-10 for one the others nearly determine (its coefficient in a linear model becomes unstable). The same values as the R package usdm's vif() (R's lm and summary.lm, to the last bit), on the same cells: a raster's cells with a value in every band (5000 of them, drawn as terra's spatSample with set.seed(seed), when it has 7000 or more; weighted by the cosine of the latitude for longitude / latitude grids), or the complete rows of a table.

Parameters:

Name	Type	Description	Default
`env`		Raster stack (band, y, x) or a table with one column per predictor.	required
`area`		Polygons the raster is cut to first (the accessible area).	`None`
`size`	`int \| None`	Cells sampled from rasters of 7000 cells or more (None: all).	`_SAMPLE_SIZE`
`seed`	`int \| None`	Seed of R's generator for the sample (default random).	`None`

Returns:

Type	Description
`DataFrame`	`variable` and `vif`.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> ls.sdm.vif(env, seed=1).round(1).head(3)    # bio7 = bio5 - bio6: exactly collinear
  variable      vif
0     bio1      inf
1     bio5   9952.0
2     bio6  40802.8

Background and data¶

background ¶

background(
    env,
    n: int = 10000,
    *,
    occ=None,
    area=None,
    bias=None,
    distance=None,
    envelope=None,
    exclude_occurrences: bool = True,
    tryf: float = 3,
    lonlat_correction: bool = True,
    seed: int | None = None,
    x: str | None = None,
    y: str | None = None,
    crs=None,
)

Background points: a sample of the conditions available to the species.

Presence-only models (MaxEnt, presence-background regressions, machine learning) contrast the presences with a sample of the cells the species could have been found in. The sample is drawn as dismo::randomPoints draws it, with the same cells for the same seed: cells with a value in every predictor, without replacement, on longitude / latitude grids weighted by the cells' east-west width (so that each square kilometre counts the same), leaving out the cells of the occurrences. Where to draw from is the choice that matters most (Barve et al. 2011; VanDerWal et al. 2009): restrict it to the accessible area, to a distance band around the records, or weight it by sampling effort.

Parameters:

Name	Type	Description	Default
`env`		Predictors (raster stack or one band): the grid and the cells with data.	required
`n`	`int`	Number of points (fewer, with a warning, when the area has fewer cells).	`10000`
`occ`		Occurrence records: their cells are left out (`exclude_occurrences`) and they anchor `distance` and `envelope`.	`None`
`area`		Polygons to draw from (the accessible area, see `accessible_area()`).	`None`
`bias`		A raster of sampling weights on the grid of `env` (e.g. the density of records of the whole taxonomic group: "target-group background", Phillips et al. 2009, which cancels the sampling bias of the presences). `dismo::randomPoints(prob = TRUE)`.	`None`
`distance`		Metres: only cells within this distance of an occurrence, or a pair `(min, max)` for cells farther than `min` and at most `max` (biomod2's "disk" pseudo-absences; longitude / latitude grids in metres on the sphere, projected ones in their units).	`None`
`envelope`		Leave out the cells inside the occurrences' environmental envelope (every predictor between its `envelope` and `1 - envelope` quantiles at the occurrences): biomod2's "SRE" pseudo-absences, e.g. 0.025.	`None`
`exclude_occurrences`	`bool`	Leave out the cells holding occurrences (dismo's `excludep`).	`True`
`tryf`	`float`	Candidates drawn per point before excluding the occurrence cells (dismo's `tryf`).	`3`
`lonlat_correction`	`bool`	Weight longitude / latitude cells by their area (dismo's `lonlatCorrection`).	`True`
`seed`	`int \| None`	Seed of R's generator (default random; recorded in `attrs`).	`None`
`x`	`str \| None`	Coordinate columns and CRS of a table of occurrences.	`None`
`y`	`str \| None`	Coordinate columns and CRS of a table of occurrences.	`None`
`crs`	`str \| None`	Coordinate columns and CRS of a table of occurrences.	`None`

Returns:

Type	Description
`GeoDataFrame`	Points at the centres of the cells drawn (CRS of `env`), with their `cell` (row-major, 0-based).

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> m = ls.sdm.accessible_area(occ, distance=500_000)
>>> bg = ls.sdm.background(env, 1000, occ=occ, area=m, seed=1)
>>> len(bg), bg.crs.to_epsg()
(1000, 4326)

dataset ¶

dataset(
    occ,
    background,
    env,
    *,
    weights: str | None = "balanced",
    folds=None,
    categorical=None,
    x: str | None = None,
    y: str | None = None,
    crs=None,
) -> SDMData

Put presences, background points and predictor values together.

Parameters:

Name	Type	Description	Default
`occ`		Presence records (GeoDataFrame of points, or a table with coordinate columns). A table with a column `y` / `presence` of 1 and 0 is taken as presences and absences, and `background` may be None.	required
`background`		Background points or absences (e.g. `background()`'s), or None.	required
`env`		Predictors (raster stack): their values at the points become `X`.	required
`weights`	`('balanced', 'equal', None)`	`"balanced"` (default): presences weigh 1 and the background `n_presences / n_background` each, so that both sides weigh the same (as recommended for presence-background data, Barbet-Massin et al. 2012); `"equal"` / None: all 1.	`"balanced"`
`folds`		Cross-validation folds, one per row (or set them later with `with_folds()`).	`None`
`categorical`		Names of categorical predictors (codes, not quantities).	`None`
`x`	`str \| None`	Coordinate columns and CRS of tables.	`None`
`y`	`str \| None`	Coordinate columns and CRS of tables.	`None`
`crs`	`str \| None`	Coordinate columns and CRS of tables.	`None`

Returns:

Type	Description
`SDMData`	Rows without a value in some predictor are dropped (counted in `report`).

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> bg = ls.sdm.background(env, 2000, occ=occ, seed=1)
>>> data = ls.sdm.dataset(occ, bg, env)
>>> data
SDMData(116 presences, 2000 background, 9 predictors: bio1, bio5, bio6, bio7, bio8, bio9, bio12, bio16, bio17; folds: none)

SDMData `dataclass` ¶

SDMData(
    X: DataFrame,
    y: ndarray,
    coords: object,
    sample_weight: ndarray | None = None,
    folds: ndarray | None = None,
    predictors: list = list(),
    categorical: list = list(),
    area: object = None,
    report: DataFrame = DataFrame(),
    settings: dict = dict(),
    _prepare: dict | None = None,
)

Bases: Saveable

Presences and background (or absences) with their predictor values, ready for a model.

Every scikit-learn estimator takes it as it is::

model.fit(data.X, data.y, sample_weight=data.sample_weight)
GridSearchCV(model, grid, cv=data.cv).fit(data.X, data.y)

Attributes:

Name	Type	Description
`X`	`DataFrame`	One row per record, one column per predictor.
`y`	`ndarray`	1 for a presence, 0 for a background point or absence.
`coords`	`GeoDataFrame`	The points (CRS of the predictors), with `y`.
`sample_weight`	`ndarray or None`	Weight of each row (see `dataset()`).
`folds`	`ndarray or None`	Cross-validation fold of each row (1..k; 0 = in every training set).
`predictors`	`list of str`	The predictor names (columns of `X`).
`categorical`	`list of str`	The categorical predictors.
`area`	`GeoDataFrame or None`	The accessible area the background was drawn from.
`report`	`DataFrame`	One row per preparation step: records before and after, and why.
`settings`	`dict`	Every choice made (thresholds, distances, seeds): enough to repeat the preparation or to write the methods of a paper.

cv `property` ¶

cv: list[tuple[ndarray, ndarray]]

The folds as scikit-learn takes them: a list of (train rows, test rows).

Rows with fold 0 (e.g. background points in a leave-one-out design) are always in training.

with_folds ¶

with_folds(folds) -> SDMData

A copy with these folds (an array, one fold per row).

redo ¶

redo(**changes) -> SDMData

Prepare again from the same inputs, with some steps changed.

Only for data made by prepare(). data.redo(thin=40_000) repeats every step with the same settings (and seeds) except the thinning distance; any argument of prepare() can be changed.

table ¶

table() -> DataFrame

Everything as one table: coordinates, y, weight, fold and the predictors.

Cross-validation folds¶

spatial_blocks ¶

spatial_blocks(
    data,
    *,
    size: float | None = None,
    rows_cols=(10, 10),
    k: int = 5,
    hexagon: bool = True,
    flat_top: bool = False,
    selection: str = "random",
    iteration: int = 100,
    column="y",
    offset=(0.0, 0.0),
    extend: float = 0.0,
    deg_to_metre: float = 111325.0,
    seed: int | None = None,
)

Folds made of spatial blocks: no fold is a neighbour of the records it tests.

Random cross-validation puts test records next to training records, where spatial autocorrelation makes them nearly copies: the scores overestimate how well the model predicts elsewhere (Roberts et al. 2017, Ecography 40: 913-929). Blocks larger than the range of autocorrelation are assigned to folds whole. The blocks and folds of blockCV::cv_spatial (Valavi et al. 2019), the same for the same seed.

Parameters:

Name	Type	Description	Default
`data`		`SDMData` (presences and background) or a GeoDataFrame of points.	required
`size`	`float \| None`	Block width in metres (on longitude / latitude data, divided by `deg_to_metre` to degrees). Without it, `rows_cols` blocks over the points' extent.	`None`
`rows_cols`		(rows, columns) of blocks when `size` is not given (hexagons: the rows set the size).	`(10, 10)`
`k`	`int`	Number of folds.	`5`
`hexagon`	`bool`	Hexagons (blockCV's default) or squares.	`True`
`flat_top`	`bool`	Hexagons with a flat top instead of a pointy one.	`False`
`selection`	`('random', 'systematic', 'checkerboard')`	How blocks go to folds: at random (the best of `iteration` tries: the most even folds), in turn, or as a checkerboard (squares, k = 2).	`"random"`
`iteration`	`int`	Random tries.	`100`
`column`		Classes kept balanced across folds in the random selection: `"y"` (presences and background, the default for SDMData), a column of a GeoDataFrame, or None.	`'y'`
`offset`		Shift of the blocks, as fractions of `size` (x, y).	`(0.0, 0.0)`
`extend`	`float`	Widen the extent of the blocks by this percentage (up to 5).	`0.0`
`deg_to_metre`	`float`	Metres per degree, to turn `size` into degrees.	`111325.0`
`seed`	`int \| None`	Seed of R's generator (default random).	`None`

Returns:

Type	Description
`SDMData or ndarray`	The data with `folds` (1..k) or the fold of each point. Points on no block get fold 0 (training only).

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 1000, occ=occ, seed=1), env)
>>> data = ls.sdm.spatial_blocks(data, size=500_000, k=5, seed=1)
>>> [len(test) for _, test in data.cv]
[235, 229, 214, 213, 225]

autocorrelation_range ¶

autocorrelation_range(
    env, *, n: int = 5000, seed: int | None = None, deg_to_metre: float = 111325.0
)

The range of spatial autocorrelation of each predictor (and their median), in metres.

A variogram is fitted to a random sample of each band's cells (the models tried by automap: spherical, exponential, Gaussian and Matern with Stein's parameterisation over a range of smoothness, the one with the smallest weighted squared error kept); its range parameter is how far values stay alike. Spatial blocks for cross-validation should be about as wide as the median range (Roberts et al. 2017; Valavi et al. 2019). As blockCV::cv_spatial_autocor: the same sample (R's generator, set.seed(seed)) and variogram classes, and gstat's Gauss-Newton fit step by step: the ranges agree to about 1e-6 (the Matern model's derivative is numerical). Climate layers on a continent often vary so smoothly that the fitted range exceeds the study area: use the result as a guide, not blindly.

Parameters:

Name	Type	Description	Default
`env`		Predictors (raster stack).	required
`n`	`int`	Cells sampled per band.	`5000`
`seed`	`int \| None`	Seed of R's generator (default random).	`None`
`deg_to_metre`	`float`	Metres per degree, for the block size on longitude / latitude grids.	`111325.0`

Returns:

Type	Description
`dict`	`range` (median, metres), `block_size` (the range in the units of the grid: degrees on longitude / latitude) and `layers` (a table: layer, range, sill, nugget, model, kappa).

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> r = ls.sdm.autocorrelation_range(env.sel(band=["bio1", "bio12"]), n=1000, seed=1)
>>> round(r["range"] / 1000)      # km: temperature varies smoothly across the continent
226192

spatial_clusters ¶

spatial_clusters(
    data,
    k: int = 5,
    *,
    env=None,
    scale: bool = True,
    raster_cluster: bool = False,
    num_sample: int = 10000,
    algorithm: str = "Hartigan-Wong",
    seed: int | None = None,
)

Folds from k-means clusters of the records' coordinates (or of their environment).

blockCV::cv_cluster: k-means (R's kmeans, Hartigan-Wong, 25 random starts, up to 500 iterations, the same clusters for the same seed) groups the records into k spatial clusters, each a fold. With env, the clusters are made in environmental space instead (predictors standardised as terra::scale; with raster_cluster, also on num_sample random cells, so that folds cover environmental gradients): the test folds then hold conditions the training data lack.

Parameters:

Name	Type	Description	Default
`data`		`SDMData` or a GeoDataFrame of points.	required
`k`	`int`	Number of folds.	`5`
`env`		Predictors to cluster in environmental space (default: coordinates).	`None`
`scale`	`bool`	Standardise the predictors.	`True`
`raster_cluster`	`bool`	Cluster the records together with a sample of the raster's cells.	`False`
`num_sample`	`int`	Cells sampled for `raster_cluster`.	`10000`
`algorithm`	`('Hartigan-Wong', 'Lloyd', 'Forgy', 'MacQueen')`	R's k-means algorithm.	`"Hartigan-Wong"`
`seed`	`int \| None`	Seed of R's generator (default random).	`None`

Returns:

Type	Description
`SDMData or ndarray`	The data with `folds` (1..k), or the fold of each point.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 500, occ=occ, seed=1), env)
>>> np.bincount(ls.sdm.spatial_clusters(data, 4, seed=2).folds).tolist()
[0, 201, 127, 197, 91]

buffer_folds ¶

buffer_folds(
    data, size: float, *, presence_background: bool = True, add_background: bool = False
)

Leave-one-out folds with a buffer: each presence is a test, nothing within size trains.

blockCV::cv_buffer: for each target record, the training set is every record farther than size (metres; Euclidean in projected coordinates, great-circle distance on longitude / latitude) and the test set is the record itself (and, with add_background, the background within the buffer). With presence_background the targets are the presences, else every record.

Returns:

Type	Description
`list of (train, test) index arrays`	One fold per target; pass it to scikit-learn as `cv`.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 300, occ=occ, seed=1), env)
>>> cv = ls.sdm.buffer_folds(data, 200_000)
>>> len(cv), len(cv[0][1])
(116, 1)

quadrant_folds ¶

quadrant_folds(data, orientation: str = 'lat_lon')

Four folds from the occurrences split at their median latitude and longitude (ENMeval "block").

ENMeval::get.block: the occurrences are split into two halves along the first axis of orientation, each half into two along the second; the background goes to the group whose quadrant it lies in. The same groups as ENMeval.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 500, occ=occ, seed=1), env)
>>> import numpy as np
>>> np.bincount(ls.sdm.quadrant_folds(data).folds[data.y == 1])[1:].tolist()
[29, 29, 29, 29]

checkerboard_folds ¶

checkerboard_folds(data, env, aggregation=2)

Folds from a checkerboard of aggregated cells (ENMeval "checkerboard").

ENMeval::get.checkerboard: the predictor grid aggregated by aggregation cells becomes a checkerboard; records on white squares form fold 1, on black ones fold 2. With two factors (hierarchical), the squares are aggregated again and the four combinations are folds 1-4. Presences and background follow the same grid.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 500, occ=occ, seed=1), env)
>>> sorted(set(ls.sdm.checkerboard_folds(data, env, aggregation=(4, 2)).folds.tolist()))
[1, 2, 3, 4]

random_folds ¶

random_folds(data, k: int = 5, *, seed: int | None = None)

Random k-fold of the presences, background always training (ENMeval "randomkfold").

predicts::folds: groups of (nearly) equal size, assigned in the order of runif() draws, the same as R for the same seed. Not spatial: use it as a non-spatial baseline.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 300, occ=occ, seed=1), env)
>>> import numpy as np
>>> np.bincount(ls.sdm.random_folds(data, 4, seed=3).folds[data.y == 1]).tolist()
[0, 29, 29, 29, 29]

jackknife_folds ¶

jackknife_folds(data)

Leave-one-out: each presence its own fold, background always training (ENMeval "jackknife").

For very few records (fewer than about 25, Shcheglovitova & Anderson 2013).

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus")).head(10)
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 100, occ=occ, seed=1), env)
>>> len(ls.sdm.jackknife_folds(data).cv)
10

Models¶

MaxNet ¶

MaxNet(
    classes: str = "default",
    regmult: float = 1.0,
    add_samples_to_background: bool = True,
    type: str = "cloglog",
    clamp: bool = True,
    categorical=None,
    threshold: float = 0.5,
    n_threads: int | None = None,
)

Bases: _Classifier, BaseEstimator

Maxent fitted as R's maxnet does it, natively: the same model, faster.

Maxent (Phillips et al. 2006) as an infinitely weighted logistic regression (Phillips et al. 2017): presences against background points with the features of maxnet (linear l, quadratic q, product p, hinge h, threshold t, categorical), maxnet's default regularisation and glmnet's lasso path, keeping its last lambda. Every step is a C++ port of R's maxnet 0.1.4 and glmnet 4.1-10 that gives R's results bit for bit on x86-64 (coefficients, alpha, entropy, the four output types with and without clamping), many times faster (hashing of the duplicate rows, the gradient checks screened by a bound and run in parallel, fused residual updates).

A scikit-learn classifier (fit, predict_proba, predict, decision_function, get_params / set_params, cloning and pickling), usable in GridSearchCV and pipelines; scikit-learn itself is not needed.

Parameters:

Name	Type	Description	Default
`classes`	`str`	Feature classes: `"default"` (maxnet's rule by the number of presences: `"l"` below 10, `"lq"` below 15, `"lqh"` below 80, else `"lqph"`) or letters among `l`, `q`, `p`, `h`, `t`. Categorical predictors always give one indicator per level.	`'default'`
`regmult`	`float`	Regularisation multiplier (larger: smoother, fewer features).	`1.0`
`add_samples_to_background`	`bool`	Add the presences to the background, unless the same values already are there (maxnet's `addsamplestobackground`).	`True`
`type`	`('cloglog', 'logistic', 'exponential', 'link')`	Output of `predict_proba()` (column 1) and `predict()`: cloglog (default) and logistic are suitabilities in [0, 1].	`"cloglog"`
`clamp`	`bool`	Clamp predictors and features to their training ranges when predicting.	`True`
`categorical`	`list of str or int`	Columns that are categorical (besides those of pandas `category` or text dtype). Their levels are those of the category dtype, else the sorted distinct values.	`None`
`threshold`	`float`	Suitability at or above which `predict()` gives 1 (e.g. one of `threshold()`'s).	`0.5`
`n_threads`	`int`	Threads (None or 0: all cores but one). The result does not depend on it.	`None`

Attributes:

Name	Type	Description
`betas_`	`dict`	Feature name (maxnet's column names, e.g. `"hinge(bio1):12.5:28"`) -> coefficient, for the features with a nonzero coefficient.
`alpha_`	`float`	Normalising constant (the exponential output sums to 1 over the background).
`entropy_`	`float`	Entropy of the exponential output over the background.
`lambda_`	`float`	The lambda of the model (the last of glmnet's path).
`lambdas_`	`ndarray`	The 200 lambdas of the path.
`penalty_factor_`	`ndarray`	Regularisation of every feature (`feature_names_`).
`feature_names_`	`list of str`	Every feature (column of maxnet's model matrix).
`feature_classes_`	`str`	The feature classes used (`"default"` resolved).
`feature_names_in_`	`ndarray`	The predictor names (when fitted on a DataFrame).
`n_features_in_`	`int`
`classes_`	`ndarray`	`[0, 1]`.
`n_passes_`	`int`	Coordinate-descent passes of the path (glmnet's `npasses`).
`var_min_, var_max_`	`dict`	Range of each continuous predictor (clamping).
`sample_means_`	`dict`	Mean of each continuous predictor and most frequent level of each categorical one over the presences (maxnet's `samplemeans`).
`levels_`	`dict`	Levels of each categorical predictor.

Notes

sample_weight is accepted for compatibility and ignored: maxnet weighs presences 1 and background points 100. Rows with a missing value are not allowed in fit; predict_* gives NaN for them (R drops them). Duplicate presences among the background are found by exact comparison of the values (R compares formatted text when there are categorical predictors, which can miss some).

References

Phillips, S. J., Anderson, R. P., Dudík, M., Schapire, R. E. & Blair, M. E. (2017). Opening the black box: an open-source release of Maxent. Ecography 40, 887-893.

Friedman, J., Hastie, T. & Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software 33, 1-22.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> bg = ls.sdm.background(env, 5000, occ=occ, seed=1)
>>> data = ls.sdm.dataset(occ, bg, env)
>>> model = ls.sdm.MaxNet().fit(data.X, data.y)
>>> model.feature_classes_, len(model.betas_)
('lqph', 29)
>>> suitability = model.predict_proba(data.X)[:, 1]   # cloglog, in [0, 1]
>>> bool(suitability[data.y == 1].mean() > suitability[data.y == 0].mean())
True

fit ¶

fit(X, y=None, sample_weight=None)

Fit to presences (y 1) and background points (y 0).

Parameters:

Name	Type	Description	Default
`X`	`DataFrame or array - like(n, p)`	Predictors; an `SDMData` (then `y` is its `y`).	required
`y`	`array - like(n)`	1 for a presence, 0 for a background point.	`None`
`sample_weight`	`ignored`	maxnet's weights are fixed (presences 1, background 100).	`None`

Returns:

Type	Description
`MaxNet`	self

predict_output ¶

predict_output(X, type: str | None = None, clamp: bool | None = None) -> ndarray

maxnet's predict(): the output of one type for every row.

Parameters:

Name	Type	Description	Default
`X`	`DataFrame or array - like(n, p)`	Predictors (the columns of `fit`, by name for a DataFrame).	required
`type`	`('cloglog', 'logistic', 'exponential', 'link')`	Default: the estimator's `type`.	`"cloglog"`
`clamp`	`bool`	Default: the estimator's `clamp`.	`None`

Returns:

Type	Description
`ndarray(n)`	NaN where a predictor the model uses is missing.

predict_proba ¶

predict_proba(X) -> ndarray

Suitability as a probability: column 1 the type output (cloglog or logistic), column 0 its complement.

Returns:

Type	Description
`ndarray(n, 2)`

predict ¶

predict(X) -> ndarray

1 where the type output reaches threshold, else 0 (NaN rows: 0).

Returns:

Type	Description
`numpy.ndarray (n,) of int`

decision_function ¶

decision_function(X) -> ndarray

The link (linear predictor plus alpha), clamped as clamp says.

Returns:

Type	Description
`ndarray(n)`

Evaluation¶

evaluate ¶

evaluate(presences, background, *, thresholds=None) -> Evaluation

Evaluate predictions at presences against background points (or absences).

The statistics of predicts::pa_evaluate / dismo::evaluate, the same values: AUC (Mann-Whitney, ties counted one half), the point-biserial correlation and, for every candidate threshold (the prediction values minus 0.0001), the confusion matrix with kappa, sensitivity, specificity, TSS (sensitivity + specificity - 1, Allouche et al. 2006) and the rest; plus the continuous Boyce index (see boyce()) and the usual thresholds.

Parameters:

Name	Description	Default
`presences`	Predictions at the presences and at the background points (or absences); missing values are dropped.	required
`background`	Predictions at the presences and at the background points (or absences); missing values are dropped.	required
`thresholds`	Candidate thresholds (default: from the predictions, as predicts).	`None`

Returns:

Type	Description
`Evaluation`

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> p = np.array([0.9, 0.8, 0.8, 0.7, 0.6, 0.55, 0.3])
>>> a = np.array([0.8, 0.5, 0.4, 0.4, 0.3, 0.2, 0.1, 0.05])
>>> e = ls.sdm.evaluate(p, a)
>>> round(e.stats["auc"], 4), round(e.stats["max_tss"], 4)
(0.8482, 0.7321)
>>> e.thresholds["max_spec_sens"]
0.5499

Evaluation `dataclass` ¶

Evaluation(stats: dict, table: DataFrame, thresholds: dict, boyce: dict = dict())

Bases: Saveable

How well predictions separate presences from background (or absences).

Attributes:

Name	Type	Description
`stats`	`dict`	`auc`, `cor` (point-biserial), `pcor` (its p-value), `np` and `na` (numbers of presences and background), `prevalence`, `max_tss` and `boyce` (continuous Boyce index).
`table`	`DataFrame`	One row per candidate threshold: the confusion matrix (`tp`, `fp`, `fn`, `tn`), `kappa`, `ccr`, `tpr` (sensitivity), `tnr` (specificity), `fpr`, `fnr`, `ppp`, `npp`, `mcr`, `or` and `tss`.
`thresholds`	`dict`	`max_kappa`, `max_spec_sens` (= max TSS), `no_omission`, `equal_prevalence`, `equal_sens_spec`, `p10` (10th percentile of the presences) and `mtp` (minimum training presence).

boyce ¶

boyce(
    presences,
    available,
    *,
    nclass: int = 0,
    window: float | None = None,
    res: int = 100,
    method: str = "spearman",
    rm_duplicate: bool = True,
) -> dict

Continuous Boyce index: does predicted suitability rise with the frequency of presences?

For moving windows of suitability, the ratio of the share of presences to the share of available cells (background) predicted in the window; the index is the rank correlation of that ratio with suitability: 1 for a model whose higher predictions always hold relatively more presences, 0 for a random one. The standard metric for presence-only data (Boyce et al. 2002; Hirzel et al. 2006). ecospat.boyce, the same windows, ratios and correlation, except for its "+1 trick": ecospat shifts its last window by one unit of suitability, which removes that window for predictions on 0-1 and keeps it on 0-1000, so the index depends on the scale of the predictions; here the last window is [max - width, max] on any scale.

Parameters:

Name	Type	Description	Default
`presences`		Predictions at the presences (test records).	required
`available`		Predictions at the background points (or every cell of the area).	required
`nclass`	`int`	0: moving windows (default); k > 0: k classes of equal width.	`0`
`window`	`float \| None`	Width of the moving window (default: a tenth of the range of `available`).	`None`
`res`	`int`	Number of window steps (100).	`100`
`method`	`('spearman', 'pearson', 'kendall')`	Correlation.	`"spearman"`
`rm_duplicate`	`bool`	Ignore successive equal ratios in the correlation.	`True`

Returns:

Type	Description
`dict`	`cor` (the index, rounded to 3 decimals as ecospat), `f_ratio` (predicted / expected ratios), `hs` (window centres) and `used` (the windows entering the correlation).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> rng = np.random.default_rng(1)
>>> bg = rng.uniform(0, 1, 2000)
>>> pres = rng.beta(4, 1.5, 200)            # presences where suitability is high
>>> ls.sdm.boyce(pres, bg)["cor"]
0.954

threshold ¶

threshold(
    presences, background, method: str = "max_spec_sens", *, sensitivity: float = 0.9
) -> float

A threshold that turns continuous predictions into presence / absence.

Parameters:

Name	Type	Description	Default
`presences`		Predictions at the presences and the background (or absences).	required
`background`		Predictions at the presences and the background (or absences).	required
`method`	`str`	`"max_spec_sens"` (maximum TSS, the recommended default; Liu et al. 2013), `"max_kappa"`, `"no_omission"` (the highest threshold that keeps every presence), `"equal_prevalence"`, `"equal_sens_spec"`, `"sensitivity"` (the highest threshold keeping at least `sensitivity` of the presences; dismo returns the first one below it), `"p10"` (10th percentile training presence) or `"mtp"` (minimum training presence).	`'max_spec_sens'`
`sensitivity`	`float`	For `method="sensitivity"`.	`0.9`

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> p = np.array([0.9, 0.8, 0.8, 0.7, 0.6, 0.55, 0.3])
>>> a = np.array([0.8, 0.5, 0.4, 0.4, 0.3, 0.2, 0.1, 0.05])
>>> ls.sdm.threshold(p, a, "sensitivity", sensitivity=0.8)
0.5499

omission_rate ¶

omission_rate(train, test, threshold: str = 'p10') -> float

Share of test presences predicted below a threshold set on the training presences.

ENMeval's or.mtp (threshold="mtp": below the lowest training prediction) and or.10p ("p10": below the 10th percentile of the training predictions). A well-fitted model omits about 0 and 0.1 of the test presences; much more means overfitting (Radosavljevic & Anderson 2014).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> train = np.linspace(0.2, 1.0, 50)
>>> test = np.array([0.1, 0.25, 0.5, 0.9])
>>> ls.sdm.omission_rate(train, test, "mtp"), ls.sdm.omission_rate(train, test, "p10")
(0.25, 0.5)

Extrapolation¶

mess ¶

mess(reference, env, *, full: bool = False)

Multivariate environmental similarity surface: how far conditions lie outside the reference.

For each variable, the similarity of a cell's value to the distribution of the reference values (positive inside the reference range, 100 at its median; negative outside it, in percent of the range); the MESS is the smallest over the variables, so negative values mark cells where at least one variable is outside anything the model saw (Elith et al. 2010, Methods in Ecology and Evolution 1: 330-342). predicts::mess, the same values.

Parameters:

Name	Type	Description	Default
`reference`		The conditions the model was fitted to: `SDMData` (its predictor values), a table, or a raster stack (its cells).	required
`env`		The conditions to project to: raster stack (e.g. future climate) or a table, with the reference's variables.	required
`full`	`bool`	Also return the similarity of every variable (bands `<name>_mess`).	`False`

Returns:

Type	Description
`DataArray or DataFrame`	Band `mess` (and the per-variable bands), NaN where every variable is missing.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 2000, occ=occ, seed=1), env)
>>> m = ls.sdm.mess(data, env + 30)      # 3 degrees warmer (and 30 mm wetter)
>>> bool((m.sel(band="mess") < 0).sum() > 0)
True

mop ¶

mop(
    reference,
    env,
    *,
    type: str = "basic",
    distance: str | None = "euclidean",
    where: str = "in_range",
    percentage: float = 1.0,
    scale: bool = False,
    center: bool = False,
    rescale: bool = False,
)

Mobility-oriented parity: where new conditions are outside or far from the reference ones.

The variables outside the reference range (simple: how many; basic: 1 where any) and, for each cell, the mean distance to its nearest reference conditions (the closest percentage percent of them): large distances mean dissimilar conditions even within the ranges (Owens et al. 2013, Ecological Modelling 263: 10-18). The R package mop's mop(), the same values; the Mahalanobis distance uses the covariance of the reference conditions (mop takes that of the cells tested, so a cell's distance depended on which other cells were passed).

Parameters:

Name	Type	Description	Default
`reference`		Conditions the model was fitted to (`SDMData`, table or raster stack).	required
`env`		Conditions to test (raster stack or table).	required
`type`	`('basic', 'simple', 'detailed')`	`"simple"` adds the number of variables outside the range; `"detailed"` also which variables are below / above it (`low_combined` / `high_combined`: sums of 2^j for variable j, 1-based).	`"basic"`
`distance`	`('euclidean', 'mahalanobis', None)`	Distance to the reference (None: none).	`"euclidean"`
`where`	`('in_range', 'out_range', 'all')`	Cells whose distance is computed.	`"in_range"`
`percentage`	`float`	Share of the nearest reference conditions averaged, in percent.	`1.0`
`scale`	`bool`	Standardise the variables (over reference and tested conditions together, as mop) before the Euclidean distance.	`False`
`center`	`bool`	Standardise the variables (over reference and tested conditions together, as mop) before the Euclidean distance.	`False`
`rescale`	`bool`	Rescale the distances to 0-1.	`False`

Returns:

Type	Description
`DataArray or DataFrame`	Bands `mop` (distance), `basic` (and `simple`, `low_combined`, `high_combined` with the other types).

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 2000, occ=occ, seed=1), env)
>>> out = ls.sdm.mop(data, env + 30)
>>> [str(b) for b in out.band.values]
['mop', 'basic']

exdet ¶

exdet(reference, env)

Extrapolation detection: univariate (NT1) and combinatorial (NT2) novelty.

Mesgaran et al. (2014, Diversity and Distributions 20: 1147-1159): NT1 < 0 where some variable is outside its reference range (more negative, farther out); inside every range, NT2 is the Mahalanobis distance to the reference centre relative to the largest one among the reference points, > 1 for a combination of values the reference never had. mic names the variable most responsible (1-based; 0 for analogue conditions). The ExDet of the R package dsmextra, except that its combinatorial mic pointed at the wrong variable (it numbered the variables backwards).

Returns:

Type	Description
`DataArray or DataFrame`	Bands `exdet` and `mic`.

Examples:

>>> import landschaft as ls
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> data = ls.sdm.dataset(occ, ls.sdm.background(env, 2000, occ=occ, seed=1), env)
>>> out = ls.sdm.exdet(data, env.sel(band=["bio1", "bio12", "bio16"]))
>>> [str(b) for b in out.band.values]
['exdet', 'mic']

Models and projections¶

fit ¶

fit(data: SDMData, models=None, *, cv: bool = True, on=None) -> SDMResult

Fit models to SDMData, cross-validate them on its folds, and fit them to all the data.

Parameters:

Name	Type	Description	Default
`data`	`SDMData`	From `prepare()` or `dataset()` (with folds for the cross-validated scores).	required
`models`		A dict of estimators (name -> estimator) with scikit-learn's API (`fit`; `predict_proba` or `decision_function`), a list, or one estimator; default `{"maxnet": ls.sdm.MaxNet()}`. Estimators whose `fit` takes `sample_weight` get `data.sample_weight` (in a scikit-learn Pipeline, its last step does).	`None`
`cv`	`bool`	Cross-validate on `data.folds` (when there are folds).	`True`
`on`		`ls.helpers()`: fit the models on the folds here and on the helpers of this computer's group (`parallel`), with the same scores; the models travel as their scikit-learn (or `ls.sdm`) class and parameters, which must be plain values (numbers, strings, lists, other estimators), and the helpers must have the same scikit-learn. The final fits on all the data run here.	`None`

Returns:

Type	Description
`SDMResult`

Examples:

>>> import landschaft as ls
>>> from sklearn.linear_model import LogisticRegression
>>> from sklearn.pipeline import make_pipeline
>>> from sklearn.preprocessing import StandardScaler
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> data = ls.sdm.prepare(occ, env, clean=False, background=2000)
>>> glm = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))
>>> res = ls.sdm.fit(data, {"glm": glm})
>>> round(float(res.summary.loc["glm", ("auc_test", "mean")]), 2)
0.87

SDMResult `dataclass` ¶

SDMResult(
    models: dict,
    metrics: DataFrame,
    summary: DataFrame,
    thresholds: dict,
    predictors: list,
    ranges: DataFrame,
    data: SDMData | None = None,
)

Bases: Saveable

Fitted models, their cross-validated scores and what is needed to project them.

Attributes:

Name	Type	Description
`models`	`dict`	Each model fitted to all the data (name -> estimator).
`metrics`	`DataFrame`	One row per model and fold: `auc_train`, `auc_test`, `auc_diff` (overfitting), `tss` (maximum, on the test fold), `boyce` (test presences against the background), `or_mtp` and `or_10p` (omission rates).
`summary`	`DataFrame`	Their mean and standard deviation per model.
`thresholds`	`dict`	Per model, the maximum-TSS threshold of its predictions on all data.
`predictors`	`list of str`	The predictors, in the order the models take them.
`ranges`	`DataFrame`	Minimum and maximum of each predictor in the data (for clamping).

weights ¶

weights(metric: str = 'auc_test') -> dict

Ensemble weights: each model's mean test score above chance (AUC - 0.5, TSS, Boyce), 0 if below.

predict ¶

predict(
    result: SDMResult,
    env,
    *,
    models=None,
    ensemble: str | None = "weighted",
    metric: str = "auc_test",
    clamp: bool = True,
    sd: bool = True,
    block_rows: int = 256,
)

Project fitted models onto rasters of predictors (present or future conditions).

Predicted in blocks of rows, so rasters larger than memory (lazy, dask-backed: ls.load_raster(..., chunks="auto")) work; the result is the same as predicting every cell at once.

Parameters:

Name	Type	Description	Default
`result`	`SDMResult`	From `fit()`.	required
`env`		Predictors with the bands the models were fitted with.	required
`models`		Names of the models to project (default all).	`None`
`ensemble`	`('weighted', 'mean', 'median', None)`	Combine the models: weighted by their cross-validated `metric` above chance (see `weights()`), the plain mean or the median; None: only each model.	`"weighted"`
`metric`	`str`	Score weighting the ensemble (`auc_test`, `tss`, `boyce`).	`'auc_test'`
`clamp`	`bool`	Clip each predictor to its range in the data first (projections do not extrapolate the response beyond what was observed).	`True`
`sd`	`bool`	Add the standard deviation across models: where they disagree.	`True`
`block_rows`	`int`	Rows of cells predicted at a time.	`256`

Returns:

Type	Description
`DataArray`	Bands: `ensemble` and `sd` (with several models), then one per model; NaN where some predictor is missing.

Examples:

>>> import landschaft as ls
>>> from sklearn.linear_model import LogisticRegression
>>> from sklearn.pipeline import make_pipeline
>>> from sklearn.preprocessing import StandardScaler
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> data = ls.sdm.prepare(occ, env, clean=False, background=2000)
>>> res = ls.sdm.fit(data, {"glm": make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))})
>>> suit = ls.sdm.predict(res, env)
>>> [str(b) for b in suit.band.values], suit.shape
(['glm'], (1, 192, 186))

response_curves ¶

response_curves(
    result: SDMResult, *, predictors=None, n: int = 100, fixed: str = "mean"
) -> DataFrame

How each model's suitability changes along one predictor, the others held at their mean (or median).

Returns:

Type	Description
`DataFrame`	Long table: `model`, `predictor`, `value`, `suitability`.

Examples:

>>> import landschaft as ls
>>> from sklearn.linear_model import LogisticRegression
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> data = ls.sdm.prepare(occ, env, clean=False, background=2000, folds=False)
>>> res = ls.sdm.fit(data, {"glm": LogisticRegression(max_iter=5000)})
>>> curves = ls.sdm.response_curves(res, n=20)
>>> curves.groupby("predictor").size().tolist()
[20, 20, 20, 20, 20]

importance ¶

importance(result: SDMResult, *, n_repeats: int = 5, seed: int = 1) -> DataFrame

Permutation importance: how much AUC drops when a predictor's values are shuffled.

On the data the models were fitted to; n_repeats shuffles per predictor (the mean and standard deviation of the drop).

Returns:

Type	Description
`DataFrame`	`model`, `predictor`, `importance` (mean drop in AUC), `sd`.

Examples:

>>> import landschaft as ls
>>> from sklearn.linear_model import LogisticRegression
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> data = ls.sdm.prepare(occ, env, clean=False, background=2000, folds=False)
>>> res = ls.sdm.fit(data, {"glm": LogisticRegression(max_iter=5000)})
>>> imp = ls.sdm.importance(res)
>>> imp.columns.tolist()
['model', 'predictor', 'importance', 'sd']

Area of habitat¶

area_of_habitat ¶

area_of_habitat(
    ranges,
    land_cover,
    natural: Iterable[str],
    *,
    suitability=None,
    original=None,
    elevation=None,
    elevation_limits=None,
    include_unlisted: bool = True,
    cell_area=None,
    loss: bool = False,
    n_threads: int | None = None,
) -> HabitatAreas

Area of habitat of species: their ranges refined by habitats and elevations.

The area of habitat (AOH; Brooks et al. 2019) of a species is the part of its range covered by the land-cover classes it lives in, within its elevation limits. With land cover given as shares of classes per cell (as plangea's restoration model uses it), every cell of the range counts for its area times the share of the species' natural classes. The potential habitat is the same with the original vegetation, and the gain map says how much habitat a species wins per unit of a cell restored: what a restoration plan needs (plangea's harmonize_bd).

Parameters:

Name	Type	Description	Default
`ranges`		The species' ranges: a `(band, y, x)` raster whose band names are the species (values above 0 inside the range), or a mapping or list of named rasters.	required
`land_cover`		The current share of every class in each cell (as `original_vegetation()` takes it); NaN counts as 0.	required
`natural`	`Iterable[str]`	The natural classes: a species' habitat is among them.	required
`suitability`		The classes each species lives in: a table of two columns (species, class), a mapping species -> classes. Classes that are not natural are left out.	`None`
`original`		The original vegetation: a `OriginalVegetation` (its original cover gives the potential habitat, its restored proportions and anthropic share the gain) or the original share of every natural class (then restored land returns to those proportions). Without it the potential habitat is the current one.	`None`
`elevation`		An elevation raster on the same grid.	`None`
`elevation_limits`		A table with the species and their lower and upper elevations (columns `species`, `lower`, `upper`, or the first three columns in that order); a missing limit is no limit. Cells without an elevation stay in the range.	`None`
`include_unlisted`	`bool`	Species without a row in `suitability` live in every natural class (True) or are left out (False, `valid` False).	`True`
`cell_area`		The area of a cell: a number, a raster, or None: from the grid in km² (see `cell_area_km2()`).	`None`
`loss`	`bool`	Also the habitat lost per unit of natural land converted (`loss`): what conservation and conversion plans need.	`False`
`n_threads`	`int \| None`	Threads (default: all cores but one).	`None`

Returns:

Type	Description
`HabitatAreas`	The areas by species and the gain matrix.

Notes

The sums are plangea's, bit for bit: the shares of a species' classes added in their order, times the cell's area, summed over its range in R's 80-bit sum().

Examples:

>>> import numpy as np
>>> import xarray as xr
>>> import landschaft as ls
>>> cover = xr.DataArray(np.array([[[0.6, 0.0]], [[0.2, 0.1]], [[0.2, 0.9]]]),
...                      dims=("band", "y", "x"), coords={"band": ["forest", "savanna", "crops"]})
>>> ranges = xr.DataArray(np.array([[[1.0, 1.0]], [[0.0, 1.0]]]),
...                       dims=("band", "y", "x"), coords={"band": ["toucan", "rhea"]})
>>> veg = ls.original_vegetation(cover, ["forest", "savanna"], ["crops"])
>>> aoh = ls.sdm.area_of_habitat(ranges, cover, ["forest", "savanna"], original=veg, cell_area=1.0,
...                              suitability={"toucan": ["forest"], "rhea": ["savanna"]})
>>> aoh.table[["species", "current", "potential"]].round(3)
  species  current  potential
0  toucan      0.6       0.75
1    rhea      0.1       1.00

HabitatAreas `dataclass` ¶

HabitatAreas(
    table: DataFrame,
    gain: csr_matrix | None = None,
    _template: DataArray | None = None,
    loss: csr_matrix | None = None,
)

Bases: Saveable

Result of area_of_habitat().

Attributes:

Name	Type	Description
`table`	`DataFrame`	One row per species: `species`, its `current` and `potential` area of habitat (km², or the unit of `cell_area`), the number of cells of its range (`range_cells`, after the elevation limits), its natural `classes` and whether it is `valid` (some suitable natural class, current and potential habitat above 0).
`gain`	`csr_matrix \| None`	The habitat each species gains per unit of a cell restored (km² per unit of the cell's area restored): a sparse matrix, species (in the order of `table`) x cells of the grid (row-major), over the cells of its range with anthropic land. None without the original vegetation.
`loss`	`csr_matrix \| None`	The habitat each species loses per unit of a cell's area of natural land lost (converted, or not conserved): the cell's area times the share of the species' classes in its current natural cover, over the cells of its range with natural land (the same layout as `gain`). None unless asked for (`loss=True`).

gain_map ¶

gain_map(species: str) -> DataArray

The habitat a species gains per unit of each cell restored, as a raster (0 elsewhere).

cell_area_km2 ¶

cell_area_km2(template: DataArray) -> ndarray

The area of every cell of a raster's grid in km²: the product of the resolutions in a projected coordinate system (in its linear unit), the cell on the ellipsoid in a geographic one.

Example data¶

example_data ¶

example_data(name: str = 'bradypus') -> Path

Path of an example data set for species distribution modelling.

Parameters:

Name	Type	Description	Default
`name`	`('bradypus', 'acaule', 'bioclim', 'biome')`	`"bradypus"`: 116 occurrence records of the brown-throated sloth (Bradypus variegatus), a CSV table with `species`, `lon`, `lat`; `"acaule"`: 1366 GBIF records of the wild potato Solanum acaule (with country, locality, coordinate uncertainty, institution and basis of record; 284 without coordinates), the classic data for cleaning; `"bioclim"`: nine WorldClim bioclimatic variables (bio1, bio5, bio6, bio7, bio8, bio9, bio12, bio16, bio17; temperatures in tenths of a degree), 192 x 186 cells of 0.5 degree over Central and South America; `"biome"`: the WWF biomes (1-14) on the same grid.	`"bradypus"`

Returns:

Type	Description
`Path`

Examples:

>>> import landschaft as ls
>>> occ = ls.load_vector(ls.sdm.example_data("bradypus"), crs="EPSG:4326")
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> len(occ), env.shape
(116, (9, 192, 186))

Species distributions¶

The whole preparation¶

prepare ¶

Occurrence records¶

download_gbif ¶

clean_occurrences ¶

reference_data cached ¶

one_per_cell ¶

thin ¶

Study area¶

accessible_area ¶

Predictors¶

select_predictors ¶

vif ¶

Background and data¶

background ¶

dataset ¶

SDMData dataclass ¶

cv property ¶

with_folds ¶

redo ¶

table ¶

Cross-validation folds¶

spatial_blocks ¶

autocorrelation_range ¶

spatial_clusters ¶

buffer_folds ¶

quadrant_folds ¶

checkerboard_folds ¶

random_folds ¶

jackknife_folds ¶

Models¶

MaxNet ¶

fit ¶

predict_output ¶

predict_proba ¶

predict ¶

decision_function ¶

Evaluation¶

evaluate ¶

Evaluation dataclass ¶

boyce ¶

threshold ¶

omission_rate ¶

Extrapolation¶

mess ¶

mop ¶

exdet ¶

Models and projections¶

fit ¶

SDMResult dataclass ¶

weights ¶

predict ¶

response_curves ¶

importance ¶

Area of habitat¶

area_of_habitat ¶

HabitatAreas dataclass ¶

gain_map ¶

cell_area_km2 ¶

Example data¶

example_data ¶

reference_data `cached` ¶

SDMData `dataclass` ¶

cv `property` ¶

Evaluation `dataclass` ¶

SDMResult `dataclass` ¶

HabitatAreas `dataclass` ¶