Sorghastrum nutans, and the Comanche National Grassland

Project description

In this project will assess habitat suitability for a particular grass species, over a geographical area. Based on the necessary growing conditions for the species selected, Raster data will be downloaded representing the ideal and tolerance ranges for variables such as soil pH and elevation. The raster data will then be used to determine which areas are suitable or not suitable for the grass species. The study areas selected were the Comanche National Grassland and Pawnee National Grassland, and the species is Sorghastrum nutans. The study will use data from the POLARIS (soil) and SRTM (elevation) databases, and a global climate model will be used from the MACA GCM compilation. These data sources will be used together to create a derived slope raster and habitat suitability rasters for Sorghastrum nutans.

The Comanche National Grassland

The Comanche National Grassland in southeast Colorado consists of about 450,000 acres of land managed by the USDA Forest Service. It is split into two units. The Timpas region is further West, and possibly named after a Spanish word for a bar of iron, which is a reference to geological formations containing iron in that area. The name Carrizo refers to common reed grass, Phragmites australis. The Comanche National Grassland is named in honor of the Comanche tribe, which was centrally located here for 150 years, from approximately 1700-1850. The Comanche arrived in this area (indirectly) from Mexico, which is evident partly from words in their language such as "tecolotl", a word meaning "owl" which stems from the Aztec language Nahuatl. The Comanche competed over this area with the Ute and Apache tribes, which were eventually pushed out or migrated elsewhere. They eventually began to compete with Mexico as well, raiding Mexican settlements. These grasslands became the Comanche heartland, until they were asked to sign a treaty in 1853, and by 1880 they had ceased to have a presence there. (1)

The Pawnee National Grassland

The Pawnee National Grassland contains 193,060 acres of land located in Weld County, CO. Before European-American settlement, the land was inhabited by several indigenous groups including the Pawnee, Arapaho, Comanche, Kiowa, and Cheyenne, who lived semi-nomadic lives following the migration of the bison. These first nations came under pressure due to the decreasing numbers of bison and the encroachment of European-American settlers. Eventually, the Treaty of Medicine Lodge in 1867 resulted in the forced removal of many of the first nations people. Following this, the area was occupied by European-Americans who established a cattle trade route between Texas and Wyoming, known as the Goodnight-Loving trail. Many homesteaders who settled in this area tried to plow the shortgrass prarie ecosystem in order to plant crops, however the replacement of drought-resistant and deep-rooted native grasses by shallow-rooted crop species resulted in the cycle of topsoil dessication and erosion which caused the Dust Bowl of the 1930's. Today the area is managed by the US Forest Service, and contains grazing land for cattle as well as oil and gas wells. (2, 3)

Species Descriptions

Sorghastrum nutans, colloquially known as indiangrass is a "tall, bunching sod-former, 3-8 ft. in height, with broad blue-green blades and a large, plume-like, soft, golden-brown seed head. This showy perennial’s fall color is deep orange to purple" (4). The optimal environmental conditions and tolerance ranges for Sorghastrum nutans are as follows:

• Indiangrass is adapted to deep, moist soils ranging from heavy clays to sand.
• pH range of 4.8 to 8.0. For purposes of this model, we'll assume that 4.8 and 8.0 are the tolerance range, and that optimal growing conditions are represented by the median value between these two, which is 6.4. To make this a range, I added an arbitrary buffer of 0.4, or 6.0-6.8.
• medium tolerance to salinity and drought
• adapted to periodic burning and survives by sprouting from underground rhizomes
• Occurs in areas receiving between 11-45 inches of annual precipitation. This is the tolerance range. I couldn't find a source for the optimal range, so I took the median value between 11 and 45, which is 28. I made this into a range spanning 6 inches of precipitation, or 25-31 inches for the tolerance range.
• Monocultures of Sorghastrum nutans do not respond well to prescribed burning, however in mixed stands it shows increased productivity.
• Optimal temperatures for growth are between 85-95 degrees farenheit, and lower-range tolerance is 60.

(4, 5, 6)

Changing climate, changing conditions

Plant species like Sorghastrum nutans and Phragmites australis rely on particular parameters for the suitability of their habitat, for example temperature, precipitation, and soil or ecosystem type. As the climate shifts, these parameters also shift, which changes the habitat characteristics and range of the species. In this study, I plan to use a climate model to project the distribution of these habitat parameters across the Comanche National Grassland for a future scenario.

Methods:

The purpose of the project was to assess the suitabilty of habitat within the Comanche and Pawnee national grasslands for Sorghastrum nutans. To do this, several datasets were used. Soil pH was obtained from the POLARIS 30-m probabilistic maps. Elevation data was obtained from NASA via the Earthaccess python search tool, from the Shuttle Radar Topography Mission Global 1 arc second V003. Finally, a MACA (Multivariate Adaptive Constructed Analog) model was used to project precipitation values for the region into 2096 given an "extreme" rcp85 climate scenario. These datasets were harmonized to create a combined prediction of habitat viability for Sorghastrum nutans, given researched values for the plant's requirements for optimal growth and its survival tolerances.

Writing the code for this process proceeded in several steps. In notebook 1, National Grasslands boundaries vector data was downloaded. Two study areas were chosen from the National Grasslands: The Pawnee and Comanche National Grasslands. These were separated from the larger National Grasslands dataset and were plotted on top of ESRI imagery files. In notebook 2, code was written to access POLARIS soil pH rasters. First, a template URL was constructed that could be easily edited to code for different variables and study areas. A for loop was created to identify the bounds of the study areas, and compile a list of all rasters necessary within the bounds. In notebook 2a, the rasters within the list were then printed within a combined plot.

In notebook 3, code was written to login to earthaccess, from where the SRTMGL1 NASA Shuttle Radar Topography Mission Global 1 arc second V003 dataset was selected. This was then clipped to the Comanche and Pawnee bounds, and printed with the boundary overlaid on top. The SRTM for the Comanche grassland was then reprojected into UTM 13 N along with the grassland boundary, and rioxarray tools were used to derive a slope raster. The derived slope raster and Comanche grassland bounds were then plotted together. The same would have been done for the Pawnee grassland, however as this stage a bug was encountered which could not be fixed until much later. The Pawnee grassland raster could not be plotted, and only a blank plot with an axis in the middle would print. This was eventually fixed by editing the code used to select the soil data tiles, however due to time constraints, additional analysis of the Pawnee grassland was never completed.

In notebook 4, a code loop was written to access the MACA data. The code was written to be interoperable, and easily editable for testing additional variables or climate scenarios. For this model run, the variable selected was precipitation, the start year was 2096, the end year was 2099, and the climate scenario was rcp85, a 'more extreme' climate scenario. This selection was then clipped to the bounds of the Comanche National Grassland. Notebook 5 contains experimental starter code from my professor Lily Jones. It is intended to run a fuzzy logic model. Some parameters for the model were input, and some edits were made, however the model could not be made to run successfully.

Notebook 6 harmonizes the SRTM DEM with the MACA model. First, the precipitation values from the MACA model were plotted as a raster. These values were then assigned a suitability score of 0 (unsuitable) or 1 (suitable) based on the optimal rainfall and drought survival tolerances of Sorghastrum nutans. Similarly, a suitability raster was created from the DEM raster, giving a 1 or a 0 based on the elevation ranges at which nutans can grow. Finally, these rasters were multiplied to give a combined suitability raster for the 2096 rcp85 projected climate scenario. The resulting suitability raster showed that no location within the study bounds will be habitable for nutans in 2096, if the rcp85 climate scenario does indeed come to pass.

Given time for future study, there are a variety of ways in which this work could be expanded upon. First, it would be useful to complete the entire goal of the original study, which involved plotting suitability rasters for the Pawnee National Grassland and testing the model on at least one other variable. This would not be difficult, since the code to do this already exists: It would simply have to be edited to include the new study areas and variables, or copied into a new notebook from the original study. Second, the POLARIS soil pH rasters should be used to create suitability rasters. These could then be combined with the other variables or additonal scenarios to create composite suitability rasters. Third, if the fuzzy logic model could be completed, it would provide a more detailed and nuanced version of the suitability score. Instead of providing a binary raster with values of either suitable or unsuitable, the fuzzy logic model would create a decimal suitability score between one and zero. Finally, additional grass species, additional study areas, variables, and climate scenarios could all be tested and compared.

Data Description

• Polaris: Polaris is a "probabilistic soil classification and property database over the contiguous United States at a 30-meter spatial resolution" (7). The database is an aggregation of several other sources of data including historical soil survey data and "readily available high-resolution environmental covariates". It uses modern high-speed computing to calculate soil rasters at 30 meter resolution, based on other available data. For a more detailed explanation of the probabilistic methods used to create these rasters, see "POLARIS Soiil Properties: 30-m Probabilistic Maps of Soil Properties over the Contiguous United States". (8), (Data Source: 9)

• SRTM: The Shuttle Radar Topography Mission (SRTM). The SRTM data was downloaded via Earthaccess, a Python API library that provides access to NASA's vast quantities of data (10). SRTM stands for Shuttle Radar Topography Mission, and it consisted of two radar antennas boosted to space via the space shuttle Endeavor. Together, the antennae managed to collect topographic data for 80% of the Earth's surface for the year 2000 (11). The dataset was 8.6 Terabytes, which was unprecedented at the time, and covered the areas between 60 degrees North and 56 degrees South at 30m resolution. To acquire the data, a technique called radio interferometry was used, which separated the two antennae by 60 meters and used the slight differences in their perspectives to calculate elevation. (12), (Data source: 13)

• MACA: MACA (Multivariate Adaptive Constructed Analog) datasets are based on GCM (Global Climate Model) data compiled within the CIMP5 framework. Global Climate Models are computational models that project future climate scenarios based on inputs such as ocean and atmospheric circulation, radiative heat transfer, and other factors. CIMP5 stands for the stands for phase 5 of the Climate Model Intercomparison Project. The CIMP5 project combined Global Climate Models from over 40 countries worldwide, and is currently the most extensive CIMP. However, due to the global nature of CIMP data, the resolution on these models is very low, often 100-300 km. MACA uses a statistical process to increase the resolution on these models, refining them to 4 or 6km resolution. (14), (Data source: 15), (16)

• RCP85 model: The RCP85 (Representative Concentration Pathways 8.5) represents a 'very high baseline emissions' climate scenario. The scenarios are based on radiative forcing, rather than on complex socioeconomic predictions. (model description: 17).

• USFS National Grasslands Boundaries: Boundaries for the Comanche and Pawnee National Grasslands were downloaded from the US Forest Service Enterprise Data Warehouse (EDW) (18)

Sources:

1. Donald L. Hazlett, Vascular Plant Species of the Comanche National Grassland in Southeatern Colorado, https://www.fs.usda.gov/rm/pubs/rmrs_gtr130.pdf
2. Colorado Encyclopedia, Pawnee National Grassland https://coloradoencyclopedia.org/article/pawnee-national-grassland#:~:text=History,sustained%20themselves%20by%20hunting%20bison.
3. US Forest Service, Pawnee National Grassland History https://www.fs.usda.gov/detail/arp/learning/history-culture/?cid=fsm91_058308
4. Plant Database, Sorghastrum nutans, https://www.wildflower.org/plants/result.php?id_plant=sonu2
5. Gemma Johnstone, 3/26/24, "How to Grow and Care for Wood Grass", https://www.thespruce.com/indian-grass-plant-profile-5069914
6. Natural Resources Conservation Service Plant Guide, Indiangrass, https://www.nrcs.usda.gov/plantmaterials/etpmcpg13196.pdf
7. Duke University, "Polaris - a probabilistic soil classification and property database over the contiguous United States at a 30-meter spatial resolution", https://otc.duke.edu/technologies/polaris-a-probabilistic-soil-classification-and-property-database-over-the-contiguous-united-states-at-a-30-meter-spatial-resolution/#:~:text=POLARIS%20%E2%80%93%20a%20probabilistic%20soil%20classification,meter%20spatial%20resolution%20%E2%80%94%20Duke%20OTC
8. Chaney et al., AGU Water Resources Research, "POLARIS Soil Properties: 30-m Probabilistic Maps of Soil Properties Over the Contiguous United States", https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018WR022797
9. Polaris data source, Duke University: http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/
10. Josh Blumenfeld, NASA EarthData, "Earthaccess: Earth Science Data Simplified" https://www.earthdata.nasa.gov/news/blog/earthaccess-earth-science-data-simplified
11. NASA EarthData, SRTM, https://www.earthdata.nasa.gov/data/instruments/srtm#:~:text=The%20Shuttle%20Radar%20Topography%20Mission,global%20dataset%20of%20land%20elevations.
12. NASA's Earth Observing System, Shuttle Radar Topography Mission (SRTM), https://eospso.nasa.gov/missions/shuttle-radar-topography-mission
13. SRTMGL1 NASA Shuttle Radar Topography Mission Global 1 arc second V003, Earthaccess, 12/15/2024
14. Climatology Lab, MACA, https://www.climatologylab.org/maca.html
15. Northwest Knowledge Network, http://thredds.northwestknowledge.net:8080/thredds/reacch_climate_CMIP5_macav2_catalog2.html
16. Earth Lab, 11/19/21, https://www.earthdatascience.org/courses/use-data-open-source-python/hierarchical-data-formats-hdf/intro-to-MACAv2-cmip5-data/
17. Carbon Brief, https://www.carbonbrief.org/explainer-the-high-emissions-rcp8-5-global-warming-scenario/
18. USFS National Grasslands boundaries, Enterprise Data Warehouse (EDW) 12/15/24, https://data.fs.usda.gov/geodata/edw/

# Import tools
import os
import pathlib #finds the home folder
import pandas as pd
import geopandas as gpd
import matplotlib
import hvplot.pandas
import zipfile
import numpy
import xarray as xr
import rioxarray as rio
import cartopy as ccrs
import earthaccess

Define a study area: Download the USFS National Grassland Units and select study sites. I selected Pawnee and Comanche National Grasslands.

Pseudocode outline:

• make the project directories
• conditional download grassland data only once
• use .plot to create maps of the grassland ecosystems

# make the project directories

data_dir = os.path.join(
    pathlib.Path.home(),
    'earth-analytics',
    'data',
    'habitat_suitability'
)
os.makedirs(data_dir, exist_ok=True)
data_dir

'C:\\Users\\moenc\\earth-analytics\\data\\habitat_suitability'

# Define info for grasslands download
grasslands_url = (
    "https://data.fs.usda.gov/geodata/edw/edw_resources/"
    "shp/S_USA.NationalGrassland.zip"
)
grasslands_dir = os.path.join(data_dir, 'grasslands')
os.makedirs(grasslands_dir, exist_ok=True)
grasslands_path = os.path.join(grasslands_dir, 'grasslands.shp')

# Only download once
if not os.path.exists (grasslands_path):
    grasslands_gdf = gpd.read_file(grasslands_url)
    grasslands_gdf.to_file(grasslands_path)

# Load from file
grasslands_gdf = gpd.read_file(grasslands_path)

# Check the data
grasslands_gdf.plot()

<Axes: >

Figure 1: above is a vector dataset containing all of the national grasslands. The Pawnee and Comanche national grasslands will need to be selected, and separated out from within this dataset.

# Print the full grasslands_gdf. Note the NATIONALGR column. 
grasslands_gdf

	NATIONALGR	GRASSLANDN	GIS_ACRES	SHAPE_AREA	SHAPE_LEN	geometry
0	281771010328	Fort Pierre National Grassland	209044.225	0.095149	1.455518	POLYGON ((-100.08409 44.28162, -100.08409 44.2...
1	295507010328	Butte Valley National Grassland	19489.170	0.008557	0.853736	MULTIPOLYGON (((-121.996 41.84049, -121.996 41...
2	295508010328	Kiowa National Grassland	144281.321	0.058543	9.858642	MULTIPOLYGON (((-104.30414 36.08063, -104.3041...
3	295509010328	Sheyenne National Grassland	70428.175	0.033356	4.097398	MULTIPOLYGON (((-97.31081 46.51457, -97.30559 ...
4	295510010328	Cedar River National Grassland	6717.517	0.003157	0.999947	MULTIPOLYGON (((-101.82221 45.95896, -101.8170...
5	295511010328	Black Kettle National Grassland	33103.349	0.013340	4.561060	MULTIPOLYGON (((-99.91659 35.71892, -99.91661 ...
6	295512010328	Rita Blanca National Grassland	94127.091	0.038271	6.652113	MULTIPOLYGON (((-102.62993 36.44072, -102.6298...
7	295513010328	Thunder Basin National Grassland	626249.208	0.282888	44.088050	MULTIPOLYGON (((-105.46005 43.31908, -105.4601...
8	295514010328	McClellan Creek National Grassland	1401.715	0.000562	0.115902	POLYGON ((-100.86003 35.20951, -100.86008 35.2...
9	295515010328	Caddo National Grassland	68479.549	0.026940	1.159342	MULTIPOLYGON (((-95.85492 33.79814, -95.85494 ...
10	295516010328	Lyndon B. Johnson National Grassland	115408.461	0.045238	1.177130	POLYGON ((-97.67964 33.43396, -97.67951 33.433...
11	295517010328	Grand River National Grassland	154665.301	0.072336	9.796405	MULTIPOLYGON (((-102.08142 45.68619, -102.0814...
12	295518010328	Buffalo Gap National Grassland	654877.177	0.295196	33.006274	MULTIPOLYGON (((-102.33839 43.93678, -102.3383...
13	295519010328	Crooked River National Grassland	173593.026	0.079499	2.815873	POLYGON ((-120.98025 44.72935, -120.9803 44.72...
14	295520010328	Little Missouri National Grassland	1025317.517	0.492748	60.421675	MULTIPOLYGON (((-104.01688 47.51793, -104.0115...
15	295521010328	Oglala National Grassland	215804.927	0.096279	1.970612	POLYGON ((-103.72477 43.001, -103.72007 43.000...
16	295522010328	Comanche National Grassland	444413.904	0.183064	26.658022	MULTIPOLYGON (((-104.02263 37.69224, -104.0225...
17	295523010328	Pawnee National Grassland	208424.885	0.089972	15.341594	MULTIPOLYGON (((-104.58106 40.82664, -104.5810...
18	295524010328	Cimarron National Grassland	109101.348	0.044765	5.539623	MULTIPOLYGON (((-101.98977 37.11933, -101.9852...
19	295525010328	Curlew National Grassland	74783.572	0.032983	1.771642	MULTIPOLYGON (((-112.6995 42.2475, -112.6995 4...

# Define and print the boundary for the Pawnee National Grassland, with ESRI imagery as the background. 
# For future runs of this code on different study areas, simply replace the "NATIONALGR" value 
# on the first line below with the corresponding value for the new study area from the table above.

pawnee_gdf = grasslands_gdf[grasslands_gdf.NATIONALGR=='295523010328']
pawnee_gdf.dissolve().hvplot(
    geo=True, tiles='EsriImagery',
    title='Pawnee National Grassland',
    fill_color=None, line_color='black', line_width=1.5,
    frame_width=800
)

c:\Users\moenc\miniconda3\envs\earth-analytics-python\Lib\site-packages\dask\dataframe\__init__.py:42: FutureWarning: 
Dask dataframe query planning is disabled because dask-expr is not installed.

You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.

  warnings.warn(msg, FutureWarning)

Figure 2: I successfully printed the boundary for the Pawnee National Grassland.

pawnee_gdf.plot()

<Axes: >

# Define and print the Comanche National Grassland boundary, with ESRI imagery as the background.

comanche_gdf = grasslands_gdf[grasslands_gdf.NATIONALGR=='295522010328']
comanche_gdf.dissolve().hvplot(
    geo=True, tiles='EsriImagery',
    title='Comanche National Grassland',
    fill_color=None, line_color='black', line_width=1.5,
    frame_width=800
)

Figure 3: I successfuly printed out the boundary for the Comanche National Grassland.

comanche_gdf.plot()

<Axes: >

%store comanche_gdf pawnee_gdf

Stored 'comanche_gdf' (GeoDataFrame)
Stored 'pawnee_gdf' (GeoDataFrame)

Download soil variable from the Polaris dataset

In the next two sections, I will download a soil variable (pH) from the Polaris dataset. To do this, I will write a code defining the URLs for the raster tiles that cover the extent of the area in question, so that only the correct tiles are downloaded. These raster tiles will be compiled into a list, and will then be plotted in a composite plot. This first section will define the bounds of the two grasslands, and compile the URLs for the appropriate soil rasters into a list.

# Import tools
import os
import pathlib #finds the home folder
import pandas as pd
import geopandas as gpd
import matplotlib
import hvplot.pandas
import zipfile
import numpy
import xarray as xr
import rioxarray as rxr
import cartopy as ccrs
import earthaccess
import glob
from glob import glob
from math import floor, ceil

# Load stored gdf's from the previous notebook
%store -r pawnee_gdf comanche_gdf

pawnee_gdf.bounds

	minx	miny	maxx	maxy
17	-104.791442	40.609563	-103.573286	41.001847

comanche_gdf.bounds

	minx	miny	maxx	maxy
16	-104.059313	36.994609	-102.313615	37.913807

# Define and format a url template for downloading Pawnee soil tiles
p_soil_url_template = ("http://hydrology.cee.duke.edu/POLARIS/"
                     "PROPERTIES/v1.0/ph/mean/30_60/"
                     "lat{min_lat}{max_lat}_lon{min_lon}{max_lon}.tif")
p_soil_url=p_soil_url_template.format(
    min_lat=40, max_lat=41, min_lon=-105, max_lon=-104)
p_soil_url

# Define and format a url template for downloading comanche soil tile
c_soil_url_template = ("http://hydrology.cee.duke.edu/POLARIS/"
                     "PROPERTIES/v1.0/ph/mean/30_60/"
                     "lat{min_lat}{max_lat}_lon{min_lon}{max_lon}.tif")
c_soil_url=p_soil_url_template.format(
    min_lat=36, max_lat=38, min_lon=-104, max_lon=-102)
c_soil_url

'http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ph/mean/30_60/lat3638_lon-104-102.tif'

# Open the Pawnee soil raster as a data array
soil_da = rxr.open_rasterio(
    p_soil_url,
    mask_and_scale=True
).squeeze()
soil_da

<xarray.DataArray (y: 3600, x: 3600)> Size: 52MB
[12960000 values with dtype=float32]
Coordinates:
    band         int64 8B 1
  * x            (x) float64 29kB -105.0 -105.0 -105.0 ... -104.0 -104.0 -104.0
  * y            (y) float64 29kB 41.0 41.0 41.0 41.0 ... 40.0 40.0 40.0 40.0
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area

xarray.DataArray

• y: 3600
• x: 3600

• ...

  [12960000 values with dtype=float32]

• Coordinates: (4)
  • band
    ()
    int64
    1

    array(1)

  • x
    (x)
    float64
    -105.0 -105.0 ... -104.0 -104.0

    array([-104.999861, -104.999583, -104.999306, ..., -104.000694, -104.000417,
           -104.000139])

  • y
    (y)
    float64
    41.0 41.0 41.0 ... 40.0 40.0 40.0

    array([40.999861, 40.999583, 40.999306, ..., 40.000694, 40.000417, 40.000139])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -105.0 0.00027777777777515666 0.0 41.0 0.0 -0.00027777777777515666

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-104.99986111111112, -104.99958333333333, -104.99930555555557,
            -104.9990277777778, -104.99875000000002, -104.99847222222223,
           -104.99819444444447,  -104.9979166666667, -104.99763888888891,
           -104.99736111111113,
           ...
            -104.0026388888983, -104.00236111112054, -104.00208333334275,
           -104.00180555556497,  -104.0015277777872, -104.00125000000943,
           -104.00097222223165, -104.00069444445387,  -104.0004166666761,
           -104.00013888889833],
          dtype='float64', name='x', length=3600))

  • y
    PandasIndex

    PandasIndex(Index([40.999861111111116,  40.99958333333334, 40.999305555555566,
            40.99902777777779, 40.998750000000015,  40.99847222222224,
           40.998194444444465,  40.99791666666669, 40.997638888888915,
            40.99736111111114,
           ...
             40.0026388888983,  40.00236111112053,  40.00208333334275,
            40.00180555556498,   40.0015277777872,  40.00125000000943,
            40.00097222223165,  40.00069444445388,   40.0004166666761,
            40.00013888889833],
          dtype='float64', name='y', length=3600))

• Attributes: (1)

  AREA_OR_POINT :

  Area

# Define the bounds of the Pawnee and Comanche National Grasslands, in order to identify which tiles to download
p_bounds_min_lon, p_bounds_min_lat, p_bounds_max_lon, p_bounds_max_lat = (
    pawnee_gdf.total_bounds
)
c_bounds_min_lon, c_bounds_min_lat, c_bounds_max_lon, c_bounds_max_lat = (
    comanche_gdf.total_bounds
)

# Create a list of soil urls needed for the Pawnee gdf
p_soil_url_list=[]
for min_lon in range (floor(p_bounds_min_lon), ceil(p_bounds_max_lon)):
    for min_lat in range (floor(p_bounds_min_lat), ceil(p_bounds_max_lat)):
        soil_url = p_soil_url_template.format(
            min_lat = min_lat, max_lat = min_lat+1,
            min_lon=min_lon, max_lon = min_lon+1
        )
    p_soil_url_list.append(soil_url)
p_soil_url_list

['http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ph/mean/30_60/lat4142_lon-105-104.tif',
 'http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ph/mean/30_60/lat4142_lon-104-103.tif']

# Create a list of soil urls needed for the Comanche gdf
c_soil_url_list=[]
for min_lon in range (floor(c_bounds_min_lon), ceil(c_bounds_max_lon)):
    for min_lat in range (floor(c_bounds_min_lat), ceil(c_bounds_max_lat)):
        soil_url = c_soil_url_template.format(
            min_lat = min_lat, max_lat = min_lat+1,
            min_lon=min_lon, max_lon = min_lon+1
        )
    c_soil_url_list.append(soil_url)
c_soil_url_list

['http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ph/mean/30_60/lat3738_lon-105-104.tif',
 'http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ph/mean/30_60/lat3738_lon-104-103.tif',
 'http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ph/mean/30_60/lat3738_lon-103-102.tif']

%store p_soil_url_list
%store c_soil_url_list

Stored 'p_soil_url_list' (list)
Stored 'c_soil_url_list' (list)

# Some useful code that I don't want to delete yet

# p_bounds_min_lon, p_bounds_min_lat, p_bounds_max_lon, p_bounds_max_lat = (
#     pawnee_gdf.total_bounds
# )
# c_bounds_min_lon, c_bounds_min_lat, c_bounds_max_lon, c_bounds_max_lat = (
#     comanche_gdf.total_bounds
# )
# for min_lon in range (floor(p_bounds_min_lon), ceil(p_bounds_max_lon)):
#     print(min_lon)

First step is to set up a conditional download for the soil variable rasters. I need to:

• make sure they're going to the correct directory
• define the download URL by copying the link address and pasting it into my code (there may be a way to automate this for each of the 4 .tiff files... they each only differ by one or two characters in the online file name)
• create a conditional code that downloads the file(s) into the correct directory, only if they have not already been downloaded.
• load the data
• print the data, to make sure the code ran successfully.

# <Note to self>: This was Elsa's suggestion for opening and 
#   plotting iterated soil rasters. More useful code I don't want to delete yet.

# def soil_download(soil_url, pH, latitue, longitude):

#     # Define info for soil variable downloads
#     variable_name=[]
#     soil_url = (
#         "http://hydrology.cee.duke.edu/POLARIS/"
#         f"PROPERTIES/v1.0/silt/mean/100_200/lat4041_lon{variable_name}.tif"
#     )
#     soil_dir = os.path.join(data_dir, 'soils')
#     os.makedirs(soil_dir, exist_ok=True)
#     soil_path = os.path.join(soil_dir, 'soil.tiff')

#     # Only download once
#     if not os.path.exists (soil_path):
#         soil_raster = rio.open_rasterio(soil_url)
#         soil_raster.rio.to_raster(soil_path)

#     # Load from file
#     soil_raster = rio.open_rasterio(soil_path)

#     # Check the data
#     soil_raster.plot()

Calculate at least one derived topographic** variable** (slope or aspect) to use in your model. You probably will wish to use the xarray-spatial library, which is available in the latest earth-analytics-python environment (but will need to be installed/updated if you are working on your own machine). Note that calculated slope may not be correct if you are using a CRS with units of degrees; you should re-project into a projected coordinate system with units of meters, such as the appropriate UTM Zone.

Harmonize your data - make sure that the grids for each of your layers match up. Check out the ds.rio.reproject_match() method from rioxarray.

Build your model. You can use any model you wish, so long as you explain your choice. However, if you are not sure what to do, we recommend building a fuzzy logic model (see below).

Present your results in at least one figure for each grassland/climate scenario combination.

Make a composite plot of soil pH rasters

This section will combine the lists of raster tiles from the previous section, and plot them all together.

# Import tools
import os
import pathlib #finds the home folder
import pandas as pd
import geopandas as gpd
import matplotlib
import matplotlib.pyplot as plt
import hvplot.pandas
import zipfile
import numpy
import xarray as xr
import rioxarray as rxr
import cartopy as ccrs
import earthaccess
import glob
from glob import glob
from math import floor, ceil

# load stored grassland_gdf's 
%store -r pawnee_gdf comanche_gdf
%store -r p_soil_url_list c_soil_url_list

# Combine the lists
soil_urls = p_soil_url_list + c_soil_url_list

# Set up subplots (adjust rows and columns for layout)
fig, axes = plt.subplots(nrows=1, ncols=len(soil_urls), figsize=(20, 5), constrained_layout=True)

# Loop through each raster, open it, and plot in a subplot
cbar_mappable = None  # To store the QuadMesh object for the colorbar
for i, soil_url in enumerate(soil_urls):
    soil_da = rxr.open_rasterio(soil_url, mask_and_scale=True).squeeze()

    # Plot the raster on the corresponding subplot
    quadmesh = soil_da.plot(ax=axes[i], add_colorbar=False)
    axes[i].set_title(f"Raster {i + 1}")  # Add a title to each subplot

    # Store the QuadMesh object for the colorbar
    if cbar_mappable is None:
        cbar_mappable = quadmesh

# Add a global colorbar
fig.colorbar(cbar_mappable, ax=axes, orientation="horizontal", fraction=0.02, pad=0.1).set_label("Value")

plt.show()

Figure 4: A combined plot of rasters showing soil pH data for the Pawnee and Comanche National Grasslands. Rasters 1 and 2 cover the Pawnee grassland, and rasters 3, 4, and 5 cover the Comanche grassland.

Download and plot elevation data

In this section, I will download and plot elevation data rasters (Digital Elevation Models) for the Pawnee and Comanche National Grasslands.

import os
import re
import pathlib
from glob import glob

import matplotlib.pyplot as plt
import earthaccess
import xrspatial
import geopandas as gpd
import rioxarray as rxr
import rioxarray.merge as rxrmerge

%store -r comanche_gdf pawnee_gdf c_soil_url_list p_soil_url_list

# build project and elevation directories

data_dir = os.path.join(
    pathlib.Path.home(),
    'earth-analytics',
    'data'
)
project_dir = os.path.join(data_dir, 'habitat_suitability')
elevation_dir = os.path.join(data_dir, 'srtm')

os.makedirs(elevation_dir, exist_ok=True)

# login to earthaccess
earthaccess.login(strategy="interactive", persist=True)

<earthaccess.auth.Auth at 0x1a399daec50>

# search for the appropriate DEM

datasets = earthaccess.search_datasets(keyword='SRTM DEM', count=11)
for dataset in datasets:
    print(dataset['umm']['ShortName'], dataset['umm']['EntryTitle'])

NASADEM_SHHP NASADEM SRTM-only Height and Height Precision Mosaic Global 1 arc second V001
NASADEM_SIM NASADEM SRTM Image Mosaic Global 1 arc second V001
NASADEM_SSP NASADEM SRTM Subswath Global 1 arc second V001
C_Pools_Fluxes_CONUS_1837 CMS: Terrestrial Carbon Stocks, Emissions, and Fluxes for Conterminous US, 2001-2016
SRTMGL1 NASA Shuttle Radar Topography Mission Global 1 arc second V003
GEDI01_B GEDI L1B Geolocated Waveform Data Global Footprint Level V002
GEDI02_B GEDI L2B Canopy Cover and Vertical Profile Metrics Data Global Footprint Level V002
NASADEM_HGT NASADEM Merged DEM Global 1 arc second V001
SRTMGL3 NASA Shuttle Radar Topography Mission Global 3 arc second V003
SRTMGL1_NC NASA Shuttle Radar Topography Mission Global 1 arc second NetCDF V003
SRTMGL30 NASA Shuttle Radar Topography Mission Global 30 arc second V002

# Define a pattern to identify DEM tiles associated with Comanche National Grassland.
srtm_c_pattern = os.path.join(elevation_dir, 'N36*hgt.zip' and 'N37*hgt.zip')
bounds_c = tuple(comanche_gdf.total_bounds)
buffer = 0.25
xmin, ymin, xmax, ymax = bounds_c
bounds_buffer = (xmin-buffer, ymin-buffer, xmax+buffer, ymax+buffer)
if not glob(srtm_c_pattern):
    srtm_c_results = earthaccess.search_data(
        short_name = "SRTMGL1",
        bounding_box=bounds_buffer    
    )
    srtm_c_results = earthaccess.download(srtm_c_results, elevation_dir)

# Print the DEM for Comanche natl. Grassland
srtm_c_da_list=[]
for srtm_c_path in glob(srtm_c_pattern):
    tile_da = rxr.open_rasterio(srtm_c_path, mask_and_scale=True).squeeze()
    cropped_da = tile_da.rio.clip_box(*bounds_buffer)
    srtm_c_da_list.append(cropped_da)
    
srtm_c_da = rxrmerge.merge_arrays(srtm_c_da_list)
srtm_c_da.plot(cmap='terrain')
comanche_gdf.boundary.plot(ax=plt.gca(), color='black')

<Axes: title={'center': 'band = 1, spatial_ref = 0'}, xlabel='x', ylabel='y'>

Figure 5: Digital Elevation Model for the Comanche National Grassland, with the boundary for the grassland plotted on top.

%store srtm_c_da

Stored 'srtm_c_da' (DataArray)

# define a pattern which identifies the DEM's belongning to Pawnee National Grasslands

srtm_p_pattern = os.path.join(elevation_dir, 'N*hgt.zip')
bounds_p = tuple(pawnee_gdf.total_bounds)
if not glob(srtm_p_pattern):
    srtm_p_results = earthaccess.search_data(
        short_name = "SRTMGL1",
        bounding_box=bounds_p    
    )
    srtm_p_results = earthaccess.download(srtm_p_results, elevation_dir)

srtm_p_pattern

'C:\\Users\\moenc\\earth-analytics\\data\\srtm\\N*hgt.zip'

bounds_p

(np.float64(-104.79144209999998),
 np.float64(40.60956304000001),
 np.float64(-103.57328565),
 np.float64(41.00184675000003))

# Print the DEM for Pawnee natl. Grassland
srtm_p_da_list=[]
for srtm_p_path in glob(srtm_p_pattern):
    tile_da = rxr.open_rasterio(srtm_p_path, mask_and_scale=True).squeeze()
    try:
        cropped_da = tile_da.rio.clip_box(*bounds_p)
    except: 
        continue
    srtm_p_da_list.append(cropped_da)
    
srtm_p_da = rxrmerge.merge_arrays(srtm_p_da_list)
srtm_p_da.plot(cmap='terrain')
pawnee_gdf.boundary.plot(ax=plt.gca(), color='black')

<Axes: title={'center': 'band = 1, spatial_ref = 0'}, xlabel='x', ylabel='y'>

Figure 6: The same plot as figure 5, but for the Pawnee National Grassland. An srtm DEM (elevation) raster overlaid by the Pawnee National Grasslands boundary.

# reproject the Comanche DEM into utm 13 N crs
utm_13n_epsg = 32613
srtm_c_proj_da = srtm_c_da.rio.reproject(utm_13n_epsg)
# srtm_p_proj_da = srtm_da.to_crs()
srtm_c_proj_da.plot()

<matplotlib.collections.QuadMesh at 0x1ba838dd950>

# Reproject so units are in meters
utm13_epsg = 32613
srtm_c_proj_da = srtm_c_da.rio.reproject(utm13_epsg)
comanche_proj_gdf = comanche_gdf.to_crs(utm_13n_epsg)
bounds_proj = tuple(comanche_proj_da.total_bounds)

# Calculate slope
slope_full_da = xrspatial.slope(srtm_c_proj_da)
# slope_da = slope_full_da.rio.clip_box(*bounds_proj)
slope_da = slope_full_da.rio.clip(comanche_proj_gdf.geometry)

# Plot slope, with Comanche bounds overlay
slope_c_da.plot(cmap='terrain')
comanche_proj_gdf.boundary.plot(ax=plt.gca(), color='white', linewidth=.5)
plt.show()

Figure 7: This is a plot of the slope in meters, derived from the DEM and overlaid with the boundary of the Comanche National Grassland.

Pseudocode for creating a slope raster of the Pawnee National Grassland:

• reproject the Pawnee SRTM and Comanche bounds into UTM 13 N
• use xrspatial to calculate a slope raster from the DEM
• plot the slope raster, overlaid by the Pawnee boundary.

import os
import pathlib
import geopandas as gpd

import pandas as pd
import rioxarray as rxr

# Import packages
import numpy as np
import netCDF4
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import seaborn as sns

%store -r comanche_gdf

model_name = 'BNU-ESM'
variable_name = 'pr'
var_long_name = 'precipitation'
dir_path = 'http://thredds.northwestknowledge.net:8080/thredds/dodsC/'

# build a directory for the macav2 data

comanche_path = os.path.join(
    pathlib.Path.home(),
    'earth-analytics',
    'data',
    'habitat_suitability'
)
os.makedirs(comanche_path, exist_ok=True)
# comanche_gdf.to_file('habitat_suitability')
# comanche_gdf = gpd.read_file(comanche_path)

# Define the file path for the GeoDataFrame
comanche_file = os.path.join(comanche_path, 'comanche_data.shp')

# Save the GeoDataFrame to the file
comanche_gdf.to_file(comanche_file)

# Read the GeoDataFrame back from the file
comanche_gdf = gpd.read_file(comanche_file)

def convert_longitude(longitude):
        """Convert longitude from 0-360 to -180-180"""
        return (longitude - 360) if longitude > 180 else longitude

maca_da_list = []
for site_name, site_gdf in {'comanche': comanche_gdf}.items():
        for variable in ["pr"]:
                for start_year in [2096]:
                        end_year = start_year + 3
                        maca_url = (
                                "http://thredds.northwestknowledge.net:8080/thredds"
                                f"/dodsC/MACAV2/BNU-ESM/macav2metdata_{variable}_BNU-ESM_"
                                f"r1i1p1_rcp85_{start_year}_{end_year}_CONUS_monthly.nc"
                        # "http://thredds/fileServer/MACAV2/BNU-ESM/macav2metdata"
                        # f"_{variable}_BNU-ESM_r1i1p1_rcp45_{start_year}"
                        # f"_{end_year}_CONUS_monthly.nc"
                        )
                        maca_da = xr.open_dataset(maca_url).squeeze().precipitation
                        bounds = site_gdf.to_crs(maca_da.rio.crs).total_bounds
                        maca_da = maca_da.assign_coords(
                                lon=("lon", [convert_longitude(l) for l in maca_da.lon.values])
                        )
                        maca_da = maca_da.rio.set_spatial_dims(x_dim='lon', y_dim='lat')
                        maca_da.rio.clip_box(*bounds)
                        maca_da_list.append(dict(
                                site_name=site_name,
                                variable=variable,
                                start_year=start_year,
                                da=maca_da))
pd.DataFrame(maca_da_list)

	site_name	variable	start_year	da
0	comanche	pr	2096	[[[<xarray.DataArray 'precipitation' ()> Size:...

maca_df = pd.DataFrame(maca_da_list)
maca_df

	site_name	variable	start_year	da
0	comanche	pr	2096	[[[<xarray.DataArray 'precipitation' ()> Size:...

%store maca_df

Stored 'maca_df' (DataFrame)

maca_da.rio.clip_box(*comanche_gdf.to_crs(maca_da.rio.crs).total_bounds)
maca_da.rio.clip_box(*bounds)

<xarray.DataArray 'precipitation' (time: 48, lat: 23, lon: 43)> Size: 190kB
[47472 values with dtype=float32]
Coordinates:
  * lat      (lat) float64 184B 36.98 37.02 37.06 37.1 ... 37.81 37.85 37.9
  * time     (time) object 384B 2096-01-15 00:00:00 ... 2099-12-15 00:00:00
  * lon      (lon) float64 344B -104.1 -104.0 -104.0 ... -102.4 -102.4 -102.3
    crs      int64 8B 0
Attributes:
    long_name:      Monthly Precipitation Amount
    units:          mm
    standard_name:  precipitation
    cell_methods:   time: sum(interval: 24 hours): sum over days
    comments:       Total monthly precipitation at surface: includes both liq...
    _ChunkSizes:    [  8 115 276]

xarray.DataArray

'precipitation'

• time: 48
• lat: 23
• lon: 43

• ...

  [47472 values with dtype=float32]

• Coordinates: (4)
  • lat
    (lat)
    float64
    36.98 37.02 37.06 ... 37.85 37.9

    long_name :

    latitude

    standard_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    description :

    Latitude of the center of the grid cell

    array([36.979553, 37.021221, 37.062885, 37.104553, 37.146217, 37.187885,
           37.229549, 37.271217, 37.312881, 37.354549, 37.396214, 37.437881,
           37.479546, 37.521214, 37.562878, 37.604546, 37.64621 , 37.687878,
           37.729542, 37.77121 , 37.812874, 37.854542, 37.896206])

  • time
    (time)
    object
    2096-01-15 00:00:00 ... 2099-12-...

    description :

    days since 1900-01-01

    array([cftime.DatetimeNoLeap(2096, 1, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 2, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 3, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 4, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 5, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 6, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 7, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 8, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 9, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 10, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 11, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2096, 12, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 1, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 2, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 3, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 4, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 5, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 6, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 7, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 8, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 9, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 10, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 11, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2097, 12, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 1, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 2, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 3, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 4, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 5, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 6, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 7, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 8, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 9, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 10, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 11, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2098, 12, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 1, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 2, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 3, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 4, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 5, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 6, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 7, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 8, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 9, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 10, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 11, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2099, 12, 15, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • lon
    (lon)
    float64
    -104.1 -104.0 ... -102.4 -102.3

    axis :

    X

    long_name :

    longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([-104.064163, -104.022491, -103.980835, -103.939148, -103.897491,
           -103.855835, -103.814148, -103.772491, -103.730835, -103.689178,
           -103.647491, -103.605835, -103.564178, -103.522522, -103.480835,
           -103.439178, -103.397522, -103.355835, -103.314178, -103.272522,
           -103.230835, -103.189178, -103.147522, -103.105835, -103.064178,
           -103.022522, -102.980835, -102.939178, -102.897522, -102.855835,
           -102.814178, -102.772522, -102.730835, -102.689178, -102.647522,
           -102.605835, -102.564178, -102.522522, -102.480835, -102.439178,
           -102.397522, -102.355865, -102.314178])

  • crs
    ()
    int64
    0

    crs_wkt :

    GEOGCS["undefined",DATUM["undefined",SPHEROID["undefined",6378137,298.257223563]],PRIMEM["undefined",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    undefined

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    undefined

    geographic_crs_name :

    undefined

    horizontal_datum_name :

    undefined

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["undefined",DATUM["undefined",SPHEROID["undefined",6378137,298.257223563]],PRIMEM["undefined",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]

    GeoTransform :

    -104.0849963596889 0.04166630336216518 0.0 36.958720207214355 0.0 0.04166603088378906

    array(0)

• Indexes: (3)
  • lat
    PandasIndex

    PandasIndex(Index([ 36.97955322265625,  37.02122116088867,  37.06288528442383,
            37.10455322265625, 37.146217346191406,  37.18788528442383,
           37.229549407958984, 37.271217346191406,  37.31288146972656,
           37.354549407958984,  37.39621353149414,  37.43788146972656,
            37.47954559326172,  37.52121353149414,   37.5628776550293,
            37.60454559326172, 37.646209716796875,   37.6878776550293,
            37.72954177856445, 37.771209716796875,  37.81287384033203,
            37.85454177856445,  37.89620590209961],
          dtype='float64', name='lat'))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([2096-01-15 00:00:00, 2096-02-15 00:00:00, 2096-03-15 00:00:00,
                 2096-04-15 00:00:00, 2096-05-15 00:00:00, 2096-06-15 00:00:00,
                 2096-07-15 00:00:00, 2096-08-15 00:00:00, 2096-09-15 00:00:00,
                 2096-10-15 00:00:00, 2096-11-15 00:00:00, 2096-12-15 00:00:00,
                 2097-01-15 00:00:00, 2097-02-15 00:00:00, 2097-03-15 00:00:00,
                 2097-04-15 00:00:00, 2097-05-15 00:00:00, 2097-06-15 00:00:00,
                 2097-07-15 00:00:00, 2097-08-15 00:00:00, 2097-09-15 00:00:00,
                 2097-10-15 00:00:00, 2097-11-15 00:00:00, 2097-12-15 00:00:00,
                 2098-01-15 00:00:00, 2098-02-15 00:00:00, 2098-03-15 00:00:00,
                 2098-04-15 00:00:00, 2098-05-15 00:00:00, 2098-06-15 00:00:00,
                 2098-07-15 00:00:00, 2098-08-15 00:00:00, 2098-09-15 00:00:00,
                 2098-10-15 00:00:00, 2098-11-15 00:00:00, 2098-12-15 00:00:00,
                 2099-01-15 00:00:00, 2099-02-15 00:00:00, 2099-03-15 00:00:00,
                 2099-04-15 00:00:00, 2099-05-15 00:00:00, 2099-06-15 00:00:00,
                 2099-07-15 00:00:00, 2099-08-15 00:00:00, 2099-09-15 00:00:00,
                 2099-10-15 00:00:00, 2099-11-15 00:00:00, 2099-12-15 00:00:00],
                dtype='object', length=48, calendar='noleap', freq=None))

  • lon
    PandasIndex

    PandasIndex(Index([-104.06416320800781, -104.02249145507812,  -103.9808349609375,
           -103.93914794921875, -103.89749145507812,  -103.8558349609375,
           -103.81414794921875, -103.77249145507812,  -103.7308349609375,
           -103.68917846679688, -103.64749145507812,  -103.6058349609375,
           -103.56417846679688, -103.52252197265625,  -103.4808349609375,
           -103.43917846679688, -103.39752197265625,  -103.3558349609375,
           -103.31417846679688, -103.27252197265625,  -103.2308349609375,
           -103.18917846679688, -103.14752197265625,  -103.1058349609375,
           -103.06417846679688, -103.02252197265625,  -102.9808349609375,
           -102.93917846679688, -102.89752197265625,  -102.8558349609375,
           -102.81417846679688, -102.77252197265625,  -102.7308349609375,
           -102.68917846679688, -102.64752197265625,  -102.6058349609375,
           -102.56417846679688, -102.52252197265625,  -102.4808349609375,
           -102.43917846679688, -102.39752197265625, -102.35586547851562,
           -102.31417846679688],
          dtype='float64', name='lon'))

• Attributes: (6)

  long_name :

  Monthly Precipitation Amount

  units :

  mm

  standard_name :

  precipitation

  cell_methods :

  time: sum(interval: 24 hours): sum over days

  comments :

  Total monthly precipitation at surface: includes both liquid and solid phases from all types of clouds (both large-scale and convective)

  _ChunkSizes :

  [ 8 115 276]

Pseudocode for testing additional variables:

• replace "variable", "start_year", and "end_year" with the desired variables.
• to test a different climate scenario, replace "rcp85" with "rcp45" or another desired scenario in the maca download URL.
• re-run and possibly edit subsequent notebooks to produce new suitability rasters.

# This code didn't work. It doesn't seem to be necessary so I've commented it out, but I retained it in case it ends up being important.

# def convert_longitude(longitude):
#         return (longitude - 360) if longitude > 180 else longitude
# # example useage:
# longitudes_0_360 = [0, 90, 180, 270, 360]
# longitude_neg180_180 = [convert_longitude(lon) for lon in longitude]

# print(longitude_neg180_180)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[24], line 5
      3 # example useage:
      4 longitudes_0_360 = [0, 90, 180, 270, 360]
----> 5 longitude_neg180_180 = [convert_longitude(lon) for lon in longitude]
      7 print(longitude_neg180_180)

NameError: name 'longitude' is not defined

# Imports
import os
import rioxarray as rxr
import xarray as xr
import numpy as np

%store -r c_soil_url_list

def calculate_suitability_score(raster, optimal_value, tolerance_range):
    """
    Calculate a fuzzy suitability score (0–1) for each raster cell based on proximity to the optimal value.

    Args:
        raster (xarray.DataArray): Input raster layer.
        optimal_value (float): The optimal value for the variable.
        tolerance_range (float): The range within which values are considered suitable.

    Returns:
        xarray.DataArray: A raster of suitability scores (0–1).
    """
    # Calculate suitability scores using a fuzzy Gaussian function
    suitability = np.exp(-((raster - optimal_value) ** 2) / (2 * tolerance_range ** 2))
    return suitability

raster = c_soil_url_list
optimal_value = 6.4
tolerance_range = 3.2

def build_habitat_suitability_model(
        input_rasters, optimal_values, tolerance_ranges, output_dir, threshold=None):
    """
    Build a habitat suitability model by combining fuzzy suitability scores for each variable.

    Args:
        input_rasters (list): List of paths to input raster files representing environmental variables.
        optimal_values (list): List of optimal values for each variable.
        tolerance_ranges (list): List of tolerance ranges for each variable.
        output_dir (str): Directory to save the combined suitability raster.
        threshold (float, optional): Threshold for highlighting highly suitable areas (default is None).

    Returns:
        str: Path to the final combined suitability raster.
    """
    input_rasters = [c_soil_url_list]
    optimal_values = [(pH>=6.0 AND pH<=6.8), (pr>=25 AND pr<=31), (temp>=85 AND temp<=95), 
                        (soil_type = (sand OR clay))]
    tolerance_ranges = [(pH>=4.8 AND pH<=8.0), (pr>=11 AND pr<=45), (temp>=60)]
    output_dir = os.path.join('earth-analytics', 'data', 
                            'habitat_suitability', 'suitability_raster')
    threshold = 0.85

  Cell In[20], line 17
    optimal_values = [(pH>=6.0 AND pH<=6.8), (pr>=25 AND pr<=31), (temp>=85 AND temp<=95),
                       ^
SyntaxError: invalid syntax. Perhaps you forgot a comma?

os.makedirs(output_dir, exist_ok=True)

    # Load and calculate suitability scores for each raster
    suitability_layers = []
    for raster_path, optimal_value, tolerance_range in zip(input_rasters, optimal_values, tolerance_ranges):
        raster = rxr.open_rasterio(raster_path, masked=True).squeeze()
        suitability_layer = calculate_suitability_score(raster, optimal_value, tolerance_range)
        suitability_layers.append(suitability_layer)

    # Combine suitability scores by multiplying across all layers
    combined_suitability = suitability_layers[0]
    for layer in suitability_layers[1:]:
        combined_suitability *= layer

    # Apply a threshold if provided
    if threshold is not None:
        combined_suitability = xr.where(combined_suitability >= threshold, combined_suitability, 0)

    # Save the combined suitability raster
    output_file = os.path.join(output_dir, "combined_suitability.tif")
    combined_suitability.rio.to_raster(output_file)
    print(f"Combined suitability raster saved to: {output_file}")

    return output_file

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 1
----> 1 output_dir = os.makedirs(output_dir, exist_ok=True)
      3 # Load and calculate suitability scores for each raster
      4 suitability_layers = []

NameError: name 'output_dir' is not defined

# Example usage
if __name__ == "__main__":
    # Paths to input raster files (e.g., temperature, precipitation, soil pH)
    input_rasters = [
        "path_to_temperature_raster/temperature.tif",
        "path_to_precipitation_raster/precipitation.tif",
        "path_to_soil_ph_raster/soil_ph.tif"
    ]

    # Optimal values for Sorghastrum nutans for each variable
    optimal_values = [25.0, 1000.0, 6.5]  # Example: temperature in °C, precipitation in mm, soil pH

    # Tolerance ranges for each variable
    tolerance_ranges = [5.0, 200.0, 0.5]  # Example: acceptable deviation for each variable

    # Output directory to save the combined suitability raster
    output_dir = "path_to_output_directory"

    # Optional threshold to highlight highly suitable areas (e.g., 0.8)
    threshold = 0.8

    # Build the habitat suitability model
    combined_suitability_file = build_habitat_suitability_model(
        input_rasters, optimal_values, tolerance_ranges, output_dir, threshold
    )

    print("Habitat suitability model created:", combined_suitability_file)

import os
import pathlib
import geopandas as gpd

import pandas as pd
import rioxarray as rxr

# Import packages
import numpy as np
import netCDF4
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import seaborn as sns

# Load stored objects

%store -r srtm_c_da maca_df comanche_gdf

srtm_c_da.plot()

<matplotlib.collections.QuadMesh at 0x18d22852210>

maca_df.drop('da', axis='columns')

	site_name	variable	start_year
0	comanche	pr	2096

(maca_df[(maca_df.start_year==2096) 
         & (maca_df.variable=='pr')]
         .da.values.item()
).plot.hist()

(array([2.293443e+07, 2.565790e+05, 2.710100e+04, 5.818000e+03,
        1.643000e+03, 4.360000e+02, 1.210000e+02, 4.500000e+01,
        1.000000e+01, 9.000000e+00]),
 array([   0.        ,  297.32266235,  594.64532471,  891.96801758,
        1189.29064941, 1486.61328125, 1783.93603516, 2081.25854492,
        2378.58129883, 2675.90405273, 2973.2265625 ]),
 <BarContainer object of 10 artists>)

Figure 8: This was a histogram plot of the maca precipitation data that was used to check if the maca data was sourced correctly. Note: it is NOT clipped to the study area, so this plot should not be used to analyze Comanche grassland precipitation values.

# define "maca_2096_original_da" (in this case, maca precipitation data projected forward to the year 2096)

maca_2096_original_da = (
    maca_df
    [(maca_df.start_year==2096) & (maca_df.variable=='pr')]
    .da.values.item()
    .groupby('time.year')
    .sum()
    .min('year')
    .rio.write_crs(4326)
    .rio.set_spatial_dims(y_dim='lat', x_dim='lon')
)

maca_2096_original_da

<xarray.DataArray 'precipitation' (lat: 585, lon: 1386)> Size: 3MB
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)
Coordinates:
  * lat          (lat) float64 5kB 25.06 25.1 25.15 25.19 ... 49.31 49.35 49.4
    crs          int32 4B 1
  * lon          (lon) float64 11kB -124.8 -124.7 -124.7 ... -67.11 -67.06
    spatial_ref  int64 8B 0

xarray.DataArray

'precipitation'

• lat: 585
• lon: 1386

• 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

  array([[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)

• Coordinates: (4)
  • lat
    (lat)
    float64
    25.06 25.1 25.15 ... 49.35 49.4

    long_name :

    latitude

    standard_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    description :

    Latitude of the center of the grid cell

    array([25.063078, 25.104744, 25.14641 , ..., 49.312691, 49.354359, 49.396023])

  • crs
    ()
    int32
    1

    grid_mapping_name :

    latitude_longitude

    longitude_of_prime_meridian :

    0.0

    semi_major_axis :

    6378137.0

    inverse_flattening :

    298.257223563

    array(1, dtype=int32)

  • lon
    (lon)
    float64
    -124.8 -124.7 ... -67.11 -67.06

    array([-124.772156, -124.730499, -124.688843, ...,  -67.148071,  -67.106415,
            -67.064758])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    array(0)

• Indexes: (2)
  • lat
    PandasIndex

    PandasIndex(Index([25.063077926635742,  25.10474395751953,  25.14640998840332,
            25.18807601928711,   25.2297420501709, 25.271408081054688,
           25.313074111938477, 25.354740142822266, 25.396406173706055,
           25.438072204589844,
           ...
            49.02103042602539,  49.06269454956055,  49.10436248779297,
           49.146026611328125,  49.18769454956055,   49.2293586730957,
           49.271026611328125,  49.31269073486328,   49.3543586730957,
            49.39602279663086],
          dtype='float64', name='lat', length=585))

  • lon
    PandasIndex

    PandasIndex(Index([-124.77215576171875, -124.73049926757812,  -124.6888427734375,
           -124.64715576171875, -124.60549926757812,  -124.5638427734375,
           -124.52217102050781, -124.48049926757812,  -124.4388427734375,
           -124.39717102050781,
           ...
            -67.43975830078125,   -67.3980712890625,  -67.35641479492188,
            -67.31475830078125,   -67.2730712890625,  -67.23141479492188,
            -67.18975830078125,   -67.1480712890625,  -67.10641479492188,
            -67.06475830078125],
          dtype='float64', name='lon', length=1386))

• Attributes: (0)

maca_2096_original_da.plot()

<matplotlib.collections.QuadMesh at 0x18d2298b090>

Figure 9: Another plot used to check the maca data download. Note: the shape of the entire United States is visible, indicating that the data is not clipped to the grassland's extent!

comanche_gdf.total_bounds

array([-104.05931301,   36.99460905, -102.31361457,   37.91380663])

# reproject the comanche_gdf to the same crs as maca_2096_original_da

comanche_gdf.to_crs(maca_2096_original_da.rio.crs)

	NATIONALGR	GRASSLANDN	GIS_ACRES	SHAPE_AREA	SHAPE_LEN	geometry
16	295522010328	Comanche National Grassland	444413.904	0.183064	26.658022	MULTIPOLYGON (((-104.02263 37.69225, -104.0225...

# Clip the maca raster to the bounds of the comanche gdf.

maca_2096_cropped = maca_2096_original_da.rio.clip_box(
    *comanche_gdf.to_crs(maca_2096_original_da.rio.crs).total_bounds
    )
maca_2096_cropped.plot()

<matplotlib.collections.QuadMesh at 0x18d2288e610>

Figure 10: This is the clipped maca precipitation data. The bounds of the plot have now been clipped to the extent of the Comanche National Grasslands boundaries. This is evident from the longitude and latitude values on the axis of the plot, which match the Comanche bounds.

# Harmonize the maca data with the srtm data by reprojecting the maca raster into the same crs as srtm_c_da

maca_2096_da = maca_2096_original_da.rio.reproject_match(srtm_c_da)
maca_2096_da.plot()

<matplotlib.collections.QuadMesh at 0x18d2b339290>

Figure 11: the reprojected maca precipitation data has been harmonized to match the coordinate reference system of the srtm DEM (Digital Elevation Model).

# Input precipitation tolerance value from research. 
# For subsequent runs with different species, update this to match the new tolerance value.

maca_suitable = maca_2096_da > 280
maca_suitable.plot()

<matplotlib.collections.QuadMesh at 0x18d2b3e8e50>

Figure 12: This raster is derived from the maca precipitation data, projected forward to the year 2096. The raster gives a value of 1 (yellow) in areas that will be suitable for sorghastrum nutans in the year 2096 given the rcp85 scenario, which is a more extreme climate change scenario, and a 0 (purple) in areas that will be unsuitable. According to the model, only a small patch in the Southwest corner will be habitable by Sorghastrum nutans.

# Input elevation tolerance value from research
# For subsequent runs with different species, update this value to match the new species tolerance.

srtm_suitable_da = srtm_c_da < 1700
srtm_suitable_da.plot()

<matplotlib.collections.QuadMesh at 0x18d2278a750>

Figure 13: Elevation tolerance of Sorghastrum nutans. This raster contains a 1 (yellow) in areas that have suitable elevation for Sorghastrum nutans to grow, and a 0 (purple) for elevations that are not suitable.

# find areas with suitable precip and elevation for Sorghastrum nutans, by multiplying the suitability rasters together to create a combined raster.

suitable_2096 = maca_suitable * srtm_suitable_da
suitable_2096.plot()

<matplotlib.collections.QuadMesh at 0x18d2d081190>

Figure 14: This is the combined suitability raster for Sorghastrum nutans, obtained by multiplying the precipitation suitability raster by the elevation suitability raster. The teal square means that in an rcp85 climate scenario, nowhere in the square between latitude 37 and 38, and longitude -102.25 and -104.25 is habitable by Sorghastrum nutans.

Pseudocode to create suitability rasters from the POLARIS soil pH data:

• Reproject the soil suitability rasters into UTM 13 N
• Clip soil suitability rasters to the study area bounds
• Set the tolerance ranges. Suitable habitat for nutans is between pH 4.8 and 8.0.
• pH_suitable = 8.0 > pH < 4.8
• This could then be combined with elevation or precipitation data by multiplication to create a composite suitability raster.