deepSSF Landscape Predictions

Author

Affiliation

Scott Forrest

Queensland University of Technology, CSIRO

Abstract

As trained convolution filters can be applied to images (i.e. rasters) of any size, it is possible to generate habitat selection surfaces of any size by using the trained filters from the habitat selection subnetwork. These have been trained to extract the salient features of the input data, and produce higher probability values where an animal is likely to take its next step, which also depends on the time of day and the day of the year. Additionally, these habitat selection filters have been conditioned on the movement process of the animal, as only local covariates that are ‘available’ to the animal are used in the training data, which was the reason for the initial development of step selection functions (Fortin et al. 2005; Rhodes et al. 2005).

When applied to a larger landscape area, the resulting probability surface denotes only habitat selection, and not the expected distribution of animals, as of course when this surface is generated it ignores the movement process, and assumes that an animal could access any point in the landscape at the next step. This is an assumption of resource selection functions (RSFs), and has been called a ‘naive’ estimation by Signer, Fieberg, and Avgar (2017), because it ignores the movement dynamics. Predicting over broader areas should also be used with caution as covariate values that were not encountered during model training may produce inaccurate or misleading results. Acknowledging these assumptions, we include examples of landscape-scale habitat selection to suggest they have some utility for diagnosing model performance and understaning the influence of certain covariates. Primarily, the resultant habitat selection map provides a visual representation of what the habitat selection subnetwork has learned from the environmental covariates, and how this changes with respect to the other covariates such as the hour and day of the year. This can be used as a diagnostic tool for further model development, or as an end in itself, as it highlights the features present in the environmental covariates that the model considered to be influential in determining the next-step’s location. From these maps, we can also directly plot the correlation between the covariate values and the habitat selection prediction probability, which represents the marginal contribution of the covariate values.

Import packages

# If using Google Colab, uncomment the following line
# !pip install rasterio

import sys
print(sys.version)  # Print Python version in use

import numpy as np                                      # Array operations
import matplotlib.pyplot as plt                         # Plotting library
import torch                                            # Main PyTorch library
from datetime import datetime                           # Date/time utilities
import os                                               # Operating system 
import glob                                             # Pattern matching
import imageio.v2 as imageio                            # Image manipulation - for creating GIFs
from IPython.display import Image, display              # For plotting GIFs
import pandas as pd                                     # Data manipulation
import rasterio                                         # Geospatial raster data

import deepSSF_model as deepSSF_model    # Import the .py file containing the deepSSF model                                          

# Get today's date
today_date = datetime.today().strftime('%Y-%m-%d')

# To save the images, create a directory
base_dir = f'outputs'
os.makedirs(base_dir, exist_ok=True)
print(f"Base output directory: {base_dir}")

3.12.9 | packaged by Anaconda, Inc. | (main, Feb  6 2025, 12:55:12) [Clang 14.0.6 ]
Using mps device
Set default tensor type to float32
Base output directory: outputs

If using Google Colab, uncomment the following lines

The file directories will also need to be changed to match the location of the files in your Google Drive.

# from google.colab import drive
# drive.mount('/content/drive')

Set the device (accelerator - cuda for NVIDIA GPU or mps for Mac)

# Set the device to be used (GPU or CPU)
device = torch.accelerator.current_accelerator().type if torch.accelerator.is_available() else "cpu"
print(f"Using {device} device")

if torch.backends.mps.is_available():
    # Set default tensor type for PyTorch
    torch.set_default_dtype(torch.float32)
    print('Set default tensor type to float32')

Using mps device
Set default tensor type to float32

Use the directory where the model is stored to save the outputs

# Output directory for saving plots
base_id_dir = f'{base_dir}/model_training/djelk_derived_covs_CNN_4L8F_move_2_2025-05-26/landscape_predictions'
os.makedirs(base_id_dir, exist_ok=True)
print(f"Output directory: {base_id_dir}")

Output directory: outputs/model_training/djelk_derived_covs_CNN_4L8F_move_2_2025-05-26/landscape_predictions

Import the GPS data

We only use this for selecting a spatial extent for the area we want to predict over.

# select the id of that data that the model was trained on
buffalo_id = 2005
n_samples = 10297 # 2005 has 10297 samples

# Specify the path to your CSV file
# csv_file_path = f'../buffalo_local_data_id/buffalo_temporal_cont_{buffalo_id}_data_df_lag_1hr_n{n_samples}.csv'
csv_file_path = '/Users/scottforrest/deepSSF/buffalo_all_steps_with_paths_n103558_steps.csv'

# Read the CSV file into a DataFrame
buffalo_df = pd.read_csv(csv_file_path)

# How many samples?
print(f"Number of samples in the DataFrame: {len(buffalo_df)}")

# Display the first few rows of the DataFrame
print(buffalo_df.head())

Number of samples in the DataFrame: 103558
     id           x1_           x2_           y1_           y2_         sl_  \
0  2005  41969.310875  41921.521939 -1.435671e+06 -1.435654e+06   50.674891   
1  2005  41921.521939  41779.439594 -1.435654e+06 -1.435601e+06  151.845214   
2  2005  41779.439594  41841.203272 -1.435601e+06 -1.435635e+06   70.659861   
3  2005  41841.203272  41655.463332 -1.435635e+06 -1.435604e+06  188.309703   
4  2005  41655.463332  41618.651923 -1.435604e+06 -1.435608e+06   37.077972   

   direction_p       ta_                        t1_  \
0     2.802478  1.367942  2018-07-25 11:04:23+10:00   
1     2.781049 -0.021429  2018-07-25 12:04:39+10:00   
2    -0.507220  2.994917  2018-07-25 13:04:17+10:00   
3     2.976198 -2.799767  2018-07-25 14:04:39+10:00   
4    -3.021610  0.285377  2018-07-25 15:04:27+10:00   

                         t2_  ...   bearing  bearing_sin  bearing_cos  \
0  2018-07-25 12:04:39+10:00  ...  2.802478     0.332652    -0.943050   
1  2018-07-25 13:04:17+10:00  ...  2.781049     0.352783    -0.935705   
2  2018-07-25 14:04:39+10:00  ... -0.507220    -0.485749     0.874098   
3  2018-07-25 15:04:27+10:00  ...  2.976198     0.164641    -0.986354   
4  2018-07-25 16:04:24+10:00  ... -3.021610    -0.119695    -0.992811   

                t1_int               t2_int  \
0  1532480663000000000  1532484279000000000   
1  1532484279000000000  1532487857000000000   
2  1532487857000000000  1532491479000000000   
3  1532491479000000000  1532495067000000000   
4  1532495067000000000  1532498664000000000   

                                           ndvi_path  \
0  /Users/scottforrest/deepSSF/ndvi/ndvi_id2005_t...   
1  /Users/scottforrest/deepSSF/ndvi/ndvi_id2005_t...   
2  /Users/scottforrest/deepSSF/ndvi/ndvi_id2005_t...   
3  /Users/scottforrest/deepSSF/ndvi/ndvi_id2005_t...   
4  /Users/scottforrest/deepSSF/ndvi/ndvi_id2005_t...   

                                         canopy_path  \
0  /Users/scottforrest/deepSSF/canopy/cnpy_id2005...   
1  /Users/scottforrest/deepSSF/canopy/cnpy_id2005...   
2  /Users/scottforrest/deepSSF/canopy/cnpy_id2005...   
3  /Users/scottforrest/deepSSF/canopy/cnpy_id2005...   
4  /Users/scottforrest/deepSSF/canopy/cnpy_id2005...   

                                          herby_path  \
0  /Users/scottforrest/deepSSF/herby/herb_id2005_...   
1  /Users/scottforrest/deepSSF/herby/herb_id2005_...   
2  /Users/scottforrest/deepSSF/herby/herb_id2005_...   
3  /Users/scottforrest/deepSSF/herby/herb_id2005_...   
4  /Users/scottforrest/deepSSF/herby/herb_id2005_...   

                                          slope_path  \
0  /Users/scottforrest/deepSSF/slope/slop_id2005_...   
1  /Users/scottforrest/deepSSF/slope/slop_id2005_...   
2  /Users/scottforrest/deepSSF/slope/slop_id2005_...   
3  /Users/scottforrest/deepSSF/slope/slop_id2005_...   
4  /Users/scottforrest/deepSSF/slope/slop_id2005_...   

                                         target_path  
0  /Users/scottforrest/deepSSF/target/targ_id2005...  
1  /Users/scottforrest/deepSSF/target/targ_id2005...  
2  /Users/scottforrest/deepSSF/target/targ_id2005...  
3  /Users/scottforrest/deepSSF/target/targ_id2005...  
4  /Users/scottforrest/deepSSF/target/targ_id2005...  

[5 rows x 35 columns]

Importing spatial data

Instead of importing the stacks of local layers (one for each step), here can import the spatial covariates at any extent, which we can predict the habitat selection over. We use an extent that covers all of the observed locations, which refer to as the ‘landscape’.

NDVI

We have NDVI layers for every month of the year for 2018 and 2019, so we will just import a few examples for illustration.

file_path_ndvi_aug2018 = '../mapping/cropped rasters/ndvi_aug_2018.tif'
file_path_ndvi_feb2019 = '../mapping/cropped rasters/ndvi_feb_2019.tif'

# read the raster file
with rasterio.open(file_path_ndvi_aug2018) as src:
    # Read the raster band as separate variable
    ndvi_aug2018 = src.read(1)
    # Get the metadata of the raster
    ndvi_meta_aug2018 = src.meta
    # Get the transform of the raster for transforming pixels to the raster CRS
    raster_transform = src.transform

with rasterio.open(file_path_ndvi_feb2019) as src:
    # Read the raster band as separate variable
    ndvi_feb2019 = src.read(1)
    # Get the metadata of the raster
    ndvi_meta_feb2019 = src.meta

Prepare the NDVI data

There are a few things we need to do to prepare the landscape layers.

First, we need to ensure that there are no NA values in the data. For NDVI we will replace any NA values with -1 (which denotes water), as in our case that is typically why they were set to NA.

Secondly, the model expects the covariates to on the same scale as the training data. We will therefore scale the NDVI data using the same max and min scaling parameters as the training data. To get these, there are some min and max print statements in the deepSSF_train.ipynb script. When we plot the NDVI data below we will see that the values will no longer range from 0 to 1, which is because there are values in the landscape layers that are outside of the range of the training data.

# Check the coordinate reference system
print("NDVI metadata:")
print(ndvi_meta_aug2018)
print("\n")

# Have a look at the affine transformation parameters that are used to convert pixel 
# coordinates to geographic coordinates and vice versa
print("Affine transformation parameters:")
print(raster_transform)
print("\n")

# Check the shape (row, columns) of the raster
print("Shape of the raster:")
print(ndvi_aug2018.shape)

# Replace NaNs in the original array with -1, which represents water
ndvi_aug2018 = np.nan_to_num(ndvi_aug2018, nan=-1.0)
ndvi_feb2019 = np.nan_to_num(ndvi_feb2019, nan=-1.0)

# Create a water mask:
# For pixels where ndvi equals -1.0 (representing water), assign -np.inf.
# For all other pixels, assign 1.
water_mask = np.where(ndvi_aug2018 == -1.0, -np.inf, 1)

# from the stack of local layers (training data)
ndvi_max = 0.8220
ndvi_min = -0.2772

# Convert the numpy arrays to PyTorch tensors
ndvi_aug2018_tens = torch.from_numpy(ndvi_aug2018)
ndvi_feb2019_tens = torch.from_numpy(ndvi_feb2019)

# Normalising the data
ndvi_aug2018_norm = (ndvi_aug2018_tens - ndvi_min) / (ndvi_max - ndvi_min)
ndvi_feb2019_norm = (ndvi_feb2019_tens - ndvi_min) / (ndvi_max - ndvi_min)

plt.imshow(ndvi_aug2018_norm.numpy())
plt.colorbar()
plt.title('NDVI August 2018')
plt.show()

plt.imshow(ndvi_feb2019_norm.numpy())
plt.colorbar()
plt.title('NDVI February 2019')
plt.show()

NDVI metadata:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_wkt('PROJCS["GDA94 / Geoscience Australia Lambert",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",134],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3112"]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}


Affine transformation parameters:
| 25.00, 0.00, 0.00|
| 0.00,-25.00,-1406000.00|
| 0.00, 0.00, 1.00|


Shape of the raster:
(2280, 2400)

Canopy cover

Canopy cover is just a single static layer.

# Path to the canopy cover raster file
file_path = '../mapping/cropped rasters/canopy_cover.tif'

# read the raster file
with rasterio.open(file_path) as src:
    # Read the raster band as separate variable
    canopy_landscape = src.read(1)
    # Get the metadata of the raster
    canopy_meta = src.meta

Prepare the canopy cover data

As with the NDVI data, we need to ensure that there are no NA values in the data.

As the canopy cover values in the landscape layer are within the same range as the training data, we see that the values range from 0 to 1.

# Check the canopy metadata:
print("Canopy metadata:")
print(canopy_meta)
print("\n")

# Check the shape (rows, columns) of the canopy raster:
print("Shape of canopy raster:")
print(canopy_landscape.shape)
print("\n")

# Check for NA values in the canopy raster:
print("Number of NA values in the canopy raster:")
print(np.isnan(canopy_landscape).sum())

# Define the maximum and minimum canopy values from the stack of local layers:
canopy_max = 82.5000
canopy_min = 0.0

# Convert the canopy data from a NumPy array to a PyTorch tensor:
canopy_landscape_tens = torch.from_numpy(canopy_landscape)

# Normalise the canopy data:
canopy_landscape_norm = (canopy_landscape_tens - canopy_min) / (canopy_max - canopy_min)

# Visualise the normalised canopy cover:
plt.imshow(canopy_landscape_norm.numpy())
plt.colorbar()
plt.title('Canopy Cover')
plt.show()

Canopy metadata:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.3999999521443642e+38, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_wkt('PROJCS["GDA94 / Geoscience Australia Lambert",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",134],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3112"]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}


Shape of canopy raster:
(2280, 2400)


Number of NA values in the canopy raster:
0

Herbaceous vegetation

# Path to the herbaceous vegetation raster file
file_path = '../mapping/cropped rasters/veg_herby.tif'

# read the raster file
with rasterio.open(file_path) as src:
    # Read the raster band as separate variable
    herby_landscape = src.read(1)
    # Get the metadata of the raster
    herby_meta = src.meta

# Check the herbaceous metadata:
print("Herbaceous metadata:")
print(herby_meta)
print("\n")

# Check the shape (rows, columns) of the herbaceous raster:
print("Shape of herbaceous raster:")
print(herby_landscape.shape)
print("\n")

# Check for NA values in the herby raster:
print("Number of NA values in the herbaceous vegetation raster:")
print(np.isnan(herby_landscape).sum())

# Define the maximum and minimum herbaceous values from the stack of local layers:
herby_max = 1.0
herby_min = 0.0

# Convert the herbaceous data from a NumPy array to a PyTorch tensor:
herby_landscape_tens = torch.from_numpy(herby_landscape)

# Normalize the herbaceous data:
herby_landscape_norm = (herby_landscape_tens - herby_min) / (herby_max - herby_min)

# Visualize the normalised herbaceous cover:
plt.imshow(herby_landscape_norm.numpy())
plt.colorbar()
plt.show()

Herbaceous metadata:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.3999999521443642e+38, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_wkt('PROJCS["GDA94 / Geoscience Australia Lambert",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",134],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3112"]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}


Shape of herbaceous raster:
(2280, 2400)


Number of NA values in the herbaceous vegetation raster:
0

Slope

# Path to the slope raster file
file_path = '../mapping/cropped rasters/slope.tif'

# read the raster file
with rasterio.open(file_path) as src:
    # Read the raster band as separate variable
    slope_landscape = src.read(1)
    # Get the metadata of the raster
    slope_meta = src.meta

# Check the slope metadata:
print("Slope metadata:")
print(slope_meta)
print("\n")

# Check the shape (rows, columns) of the slope landscape raster:
print("Shape of slope landscape raster:")
print(slope_landscape.shape)
print("\n")

# Check for NA values in the slope raster:
print("Number of NA values in the slope raster:")
print(np.isnan(slope_landscape).sum())

# Replace NaNs in the slope array with 0.0 (representing water):
slope_landscape = np.nan_to_num(slope_landscape, nan=0.0)

# Define the maximum and minimum slope values from the stack of local layers:
slope_max = 12.2981
slope_min = 0.0006

# Convert the slope landscape data from a NumPy array to a PyTorch tensor:
slope_landscape_tens = torch.from_numpy(slope_landscape)

# Normalize the slope landscape data:
slope_landscape_norm = (slope_landscape_tens - slope_min) / (slope_max - slope_min)

# Visualize the slope landscape (note: displaying the original tensor, not the normalised data):
plt.imshow(slope_landscape_tens.numpy())
plt.colorbar()
plt.show()

Slope metadata:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_wkt('PROJCS["GDA94 / Geoscience Australia Lambert",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",134],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3112"]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}


Shape of slope landscape raster:
(2280, 2400)


Number of NA values in the slope raster:
9356

Prepare data for running the model

Select a smaller extent of the landscape

To illustrate the approach, we will select a smaller extent of the landscape to predict over, which covers the spatial extent of the training data.

# from the buffalo data
buffer = 1250
min_x = min(buffalo_df['x1_']) - buffer
max_x = max(buffalo_df['x1_']) + buffer
min_y = min(buffalo_df['y1_']) - buffer
max_y = max(buffalo_df['y1_']) + buffer

# custom extent in epsg:3112
# min_x = 28148.969145
# max_x = 47719.496935
# min_y = -1442210.335861
# max_y = -1433133.681746

# Convert geographic coordinates to pixel coordinates
min_px, min_py = ~raster_transform * (min_x, min_y)
print(min_px, min_py)
max_px, max_py = ~raster_transform * (max_x, max_y)
print(max_px, max_py)

# Round pixel coordinates to integers
min_px, max_px, min_py, max_py = int(round(min_px)), \
    int(round(max_px)), \
        int(round(min_py)), \
            int(round(max_py))

# Print the pixel coordinates   
print(f"Min x = {min_px}, Max x = {max_px}, \nMin y = {min_py}, Max y = {max_py}")

193.3215017577893 2126.903702106494
2166.104503472696 170.63023888185853
Min x = 193, Max x = 2166, 
Min y = 2127, Max y = 171

Subset all layers

Here we want to crop the layers to the spatial extent defined above.

For plotting purposes, we convert the normalised landscape layers back to their natural scale. This is only for plotting and are not used for modelling.

We also calculate the minimum and maximum values of the normalised landscape subsets, which we will use to define the colour scale when plotting the layers.

# Initialize a subset array with zeros (or another padding value) with the dimensions defined by the pixel indices.
subset = np.zeros((min_py - max_py, max_px - min_px), dtype=ndvi_aug2018.dtype)

# Extract the subsets of each normalised raster using the defined pixel boundaries.
ndvi_aug2018_subset = ndvi_aug2018_norm[max_py:min_py, min_px:max_px]
ndvi_feb2019_subset = ndvi_feb2019_norm[max_py:min_py, min_px:max_px]
canopy_subset = canopy_landscape_norm[max_py:min_py, min_px:max_px]
herby_subset = herby_landscape_norm[max_py:min_py, min_px:max_px]
slope_subset = slope_landscape_norm[max_py:min_py, min_px:max_px]

# Convert normalised NDVI values back to their natural scale by reversing the normalisation.
ndvi_aug2018_subset_natural = ndvi_aug2018_subset * (ndvi_max - ndvi_min) + ndvi_min
ndvi_feb2019_subset_natural = ndvi_feb2019_subset * (ndvi_max - ndvi_min) + ndvi_min

# Convert normalised canopy, herbaceous, and slope values back to their natural scales.
canopy_subset_natural = canopy_subset * (canopy_max - canopy_min) + canopy_min
herby_subset_natural = herby_subset * (herby_max - herby_min) + herby_min
slope_subset_natural = slope_subset * (slope_max - slope_min) + slope_min

# Determine the minimum NDVI value across the two layers (natural scale).
ndvi_plot_min = min(ndvi_aug2018_subset_natural.min(), 
                    ndvi_feb2019_subset_natural.min())
print(f"Minimum NDVI value across layers {ndvi_plot_min}")

# Determine the maximum NDVI value across the three layers (natural scale).
ndvi_plot_max = max(ndvi_aug2018_subset_natural.max(), 
                    ndvi_feb2019_subset_natural.max())
print(f"Maximum NDVI value across layers {ndvi_plot_max}")

Minimum NDVI value across layers -1.0
Maximum NDVI value across layers 0.8156132698059082

Plot the landscape subsets

# NDVI August 2018
plt.imshow(ndvi_aug2018_subset_natural, cmap='viridis', vmin=ndvi_plot_min, vmax=ndvi_plot_max)
plt.colorbar(shrink=0.7)
plt.title(f'NDVI August 2018')
plt.savefig(f"{base_id_dir}/id{buffalo_id}_ndvi_aug2018.png", 
            dpi=600, bbox_inches='tight')
plt.show()
plt.close()  # Close the figure to free memory

# NDVI February 2019
plt.imshow(ndvi_feb2019_subset_natural, cmap='viridis', vmin=ndvi_plot_min, vmax=ndvi_plot_max)
plt.colorbar(shrink=0.7)
plt.title(f'NDVI February 2019')
plt.savefig(f"{base_id_dir}/id{buffalo_id}_ndvi_feb2019.png", 
            dpi=600, bbox_inches='tight')
plt.show()
plt.close()  # Close the figure to free memory

# Canopy cover
plt.imshow(canopy_subset_natural, cmap='viridis')
plt.colorbar(shrink=0.7)
plt.title(f'Canopy cover')
plt.savefig(f"{base_id_dir}/id{buffalo_id}_canopy_cover.png", 
            dpi=600, bbox_inches='tight')
plt.show()
plt.close()  # Close the figure to free memory

# Herbaceous vegetation
plt.imshow(herby_subset_natural, cmap='viridis')
plt.colorbar(shrink=0.7)
plt.title(f'Herbaceous vegetation')
plt.savefig(f"{base_id_dir}/id{buffalo_id}_herby_veg.png", 
            dpi=600, bbox_inches='tight')
plt.show()
plt.close()  # Close the figure to free memory

# Slope
plt.imshow(slope_subset_natural, cmap='viridis')
plt.colorbar(shrink=0.7)
plt.title(f'Slope')
plt.savefig(f"{base_id_dir}/id{buffalo_id}_slope.png", 
            dpi=600, bbox_inches='tight')
plt.show()
plt.close()  # Close the figure to free memory

Load the model

Define the parameters for the model

Here we enter the specific parameter values and hyperparameters for the model. These are the values that will be used to instantiate the model.

n_max_pool_layers = 2 # used to determine the number of inputs entering the fully connected block - needs to be manually changed if the number of max pooling layers is changed
n_spatial_inputs = 4 # number of spatial inputs layers 
n_scalar_inputs = 4 # number of scalar inputs 

params_dict = {"batch_size": 32, #number of samples in each batch
               "image_dim": 101, #number of pixels along the edge of each local patch/image
               "pixel_size": 25, #number of metres along the edge of a pixel
               "input_channels": n_spatial_inputs + n_scalar_inputs, #number of spatial layers in each image + number of scalar layers that are converted to a grid
               "dim_in_nonspatial_to_grid": n_scalar_inputs, #the number of scalar predictors that are converted to a grid and appended to the spatial features
               "dense_dim_in_nonspatial": n_scalar_inputs, #change this to however many other scalar predictors you have (bearing, velocity etc)
               "kernel_size": 3, #the size of the 2D moving windows / kernels that are being learned
               "stride": 1, #the stride used when applying the kernel.  This reduces the dimension of the output if set to greater than 1
               "kernel_size_mp": 2, #the size of the kernel that is used in max pooling operations
               "stride_mp": 2, #the stride that is used in max pooling operations
               "padding": 1, #the amount of padding to apply to images prior to applying the 2D convolution
               "num_movement_params": 12, #number of parameters used to parameterise the movement kernel
               "dropout": 0.1, #the proportion of nodes that are dropped out in the dropout layers

               # hyperparameters that change the model architecture
               "output_channels": 8, #number of convolution filters to learn
               "output_channels_movement": 8, #number of convolution filters to learn for the movement kernel
               "dense_dim_hidden": 128, #number of nodes in the hidden layers

               # this will be updated below
               "dense_dim_in_all": n_scalar_inputs, #number of inputs entering the fully connected block once the nonspatial features have been concatenated to the spatial features
               "device": device
               }

# Now update the dictionary with calculated values
params_dict["dense_dim_in_all"] = int(((params_dict["image_dim"] - (params_dict["image_dim"] % 2))**2) * (params_dict["output_channels_movement"] / (4**n_max_pool_layers)))

As described in the deepSSF_train.ipynb script, we saved the model definition into a file named deepSSF_model.py. We can instantiate the model by importing the file (which was done when importing other packages) and calling the classes parameter dictionary from that script.

params = deepSSF_model.ModelParams(params_dict)
model = deepSSF_model.ConvJointModel(params).to(device)
print(model)

ConvJointModel(
  (scalar_grid_output): Scalar_to_Grid_Block()
  (conv_habitat): Conv2d_block_spatial(
    (conv2d): Sequential(
      (0): Conv2d(8, 8, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (1): ReLU()
      (2): Conv2d(8, 8, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (3): ReLU()
      (4): Conv2d(8, 8, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (5): ReLU()
      (6): Conv2d(8, 1, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
    )
  )
  (conv_movement): Conv2d_block_toFC(
    (conv2d): Sequential(
      (0): Conv2d(8, 8, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (1): ReLU()
      (2): MaxPool2d(kernel_size=2, stride=2, padding=0, dilation=1, ceil_mode=False)
      (3): Conv2d(8, 8, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (4): ReLU()
      (5): MaxPool2d(kernel_size=2, stride=2, padding=0, dilation=1, ceil_mode=False)
      (6): Flatten(start_dim=1, end_dim=-1)
    )
  )
  (fcn_movement_all): FCN_block_all_movement(
    (ffn): Sequential(
      (0): Linear(in_features=5000, out_features=128, bias=True)
      (1): Dropout(p=0.1, inplace=False)
      (2): ReLU()
      (3): Linear(in_features=128, out_features=128, bias=True)
      (4): Dropout(p=0.1, inplace=False)
      (5): ReLU()
      (6): Linear(in_features=128, out_features=12, bias=True)
    )
  )
  (movement_grid_output): Params_to_Grid_Block_ChV()
)

# # load the model weights
# print(model.state_dict())
model.load_state_dict(torch.load(f'{base_id_dir}/../checkpoint_deepSSF_model.pt', 
                                 map_location=torch.device('cpu'),
                                 weights_only=True))
# print(model.state_dict())
# model.eval()

<All keys matched successfully>

Run habitat selection subnetwork on landscape layers

Select which NDVI layer to use for prediction

Now that we have layers for the full landscape area, as well as the subset defined by the observed locations, we can select which layers to use for prediction. We will use the subset layers, but just comment and uncomment the appropriate lines to use the full landscape layers.

Additionally, to run the model we need to select a single NDVI layer to predict over. For illustration, we will use the NDVI layer from August 2018, and set the input parameter corresponding to the day of the year to be within that month.

# To use the full extent
# ndvi_landscape = ndvi_aug2018_norm
# ndvi_landscape = ndvi_feb2019_norm
# canopy_landscape = canopy_landscape_norm
# herby_landscape = herby_landscape_norm
# slope_landscape = slope_landscape_norm

# To use the subset extent
ndvi_landscape = ndvi_aug2018_subset
# ndvi_landscape = ndvi_feb2019_subset
canopy_landscape = canopy_subset
herby_landscape = herby_subset
slope_landscape = slope_subset

Prepare the scalar covariates

We also need the hour of the day and day of the year to predict with (as the model was trained with these as inputs).

In this example we’ll select a single hour and a single day of the year to predict over, and then below we’ll loop over the hours of the day to illustrate the temporal variation in habitat selection that was learned by the model.

# Hour of the day (hour) 
hour_t1 = 1

# Convert hour-of-day values into sine and cosine components.
hour_t1_sin1 = np.sin(2 * np.pi * hour_t1 / 24)
hour_t1_cos1 = np.cos(2 * np.pi * hour_t1 / 24)
# hour_t1_sin2 = np.sin(4 * np.pi * hour_t1 / 24)
# hour_t1_cos2 = np.cos(4 * np.pi * hour_t1 / 24)

# Day of the year (yday)
yday_t1 = 225

# Convert day-of-year values into sine and cosine components.
yday_t1_sin1 = np.sin(2 * np.pi * yday_t1 / 365.25)
yday_t1_cos1 = np.cos(2 * np.pi * yday_t1 / 365.25)
# yday_t1_sin2 = np.sin(4 * np.pi * yday_t1 / 365.25)
# yday_t1_cos2 = np.cos(4 * np.pi * yday_t1 / 365.25)

# Convert the arrays to PyTorch tensors.
# Hour of the day and its sine/cosine components.
hour_t1_tensor = torch.tensor(hour_t1).float()
hour_t1_sin1_tensor = torch.tensor(hour_t1_sin1).float()
hour_t1_cos1_tensor = torch.tensor(hour_t1_cos1).float()
# hour_t1_sin2_tensor = torch.tensor(hour_t1_sin2).float()
# hour_t1_cos2_tensor = torch.tensor(hour_t1_cos2).float()
# Day of the year and its sine/cosine components.
yday_t1_tensor = torch.tensor(yday_t1).float()
yday_t1_sin1_tensor = torch.tensor(yday_t1_sin1).float()
yday_t1_cos1_tensor = torch.tensor(yday_t1_cos1).float()
# yday_t1_sin2_tensor = torch.tensor(yday_t1_sin2).float()
# yday_t1_cos2_tensor = torch.tensor(yday_t1_cos2).float()

def scalar_to_grid(x, dim_x, dim_y):
    """
    Reshape a scalar tensor to a grid with the given dimensions.
    
    Args:
        x (torch.Tensor): Input scalar tensor of shape [1].
        dim_x (int): Number of rows for the grid.
        dim_y (int): Number of columns for the grid.
        
    Returns:
        torch.Tensor: Tensor expanded to shape [1, 1, dim_x, dim_y].
    """
    # View x as shape [1, 1, 1, 1] and expand to a grid of shape [1, 1, dim_x, dim_y]
    scalar_map = x.view(1, 1, 1, 1).expand(1, 1, dim_x, dim_y)
    return scalar_map

# Stack the scalar tensors along a new dimension (dim=0) to form a tensor with shape [4, 1]
# representing the four covariates: hour sin, hour cos, day-of-year sin, day-of-year cos.
scalar_covariates = torch.stack([
    hour_t1_sin1_tensor, 
    hour_t1_cos1_tensor, 
    # hour_t1_sin2_tensor,
    # hour_t1_cos2_tensor,
    yday_t1_sin1_tensor,
    yday_t1_cos1_tensor
    # yday_t1_sin2_tensor,
    # yday_t1_cos2_tensor
], dim=0).to(device)
print(scalar_covariates.shape)  # Expected shape: [4, 1]

# Convert each scalar covariate into a grid matching the dimensions of the NDVI landscape.
# This creates spatial maps where each grid cell contains the same scalar value.
scalar_grids = torch.cat([
    scalar_to_grid(tensor, ndvi_landscape.shape[0], ndvi_landscape.shape[1]) 
    for tensor in scalar_covariates
], dim=1).to(device)
print(scalar_grids.shape)  # Expected shape: [1, 4, rows, cols]

torch.Size([4])
torch.Size([1, 4, 1956, 1973])

Combine the spatial and scalar (as grid) covariates

# Stack the individual landscape layers (NDVI, canopy, herbaceous, and slope) into a single tensor.
# The resulting tensor has shape [4, rows, cols].
landscape_stack = torch.stack([ndvi_landscape, canopy_landscape, herby_landscape, slope_landscape], dim=0).to(device)

# Add an extra batch dimension to the tensor.
# Now the shape becomes [1, 4, rows, cols], where 1 indicates a single sample.
landscape_stack = landscape_stack.unsqueeze(0)
print(landscape_stack.shape)

# Concatenate the landscape stack with the scalar grids along the channel dimension.
# This merges the spatial features (landscape_stack) and the repeated scalar covariate grids (scalar_grids)
# into a full feature tensor with shape [1, total_channels, rows, cols].
full_stack = torch.cat([landscape_stack, scalar_grids], dim=1)
print(full_stack.shape)

torch.Size([1, 4, 1956, 1973])
torch.Size([1, 8, 1956, 1973])

Run habitat selection subnetwork on the landscape layers

All we need to do to run the habitat selection subnetwork of the model is to pull that component out of the model and run it on the landscape layers.

# Pass the full feature stack through the habitat convolutional layer of the model to get predictions.
landscape_predictions = model.conv_habitat(full_stack)

# Print the shape of the output tensor to verify the dimensions of the predictions.
print(landscape_predictions.shape)

torch.Size([1, 1956, 1973])

Plot the predictions

Create a directory with a folder for the day of the year.

# To save the images, create a directory
base_id_day_dir = f'{base_id_dir}/yday{yday_t1}'
os.makedirs(base_id_day_dir, exist_ok=True)
print(f"Output directory: {base_id_day_dir}")

Output directory: outputs/model_training/djelk_derived_covs_CNN_4L8F_move_2_2025-05-26/landscape_predictions/yday225

As the spatial inputs are padded by the model, there are artifacts at the edges of the predictions. Sometimes this can result in quite different values to the rest of the predictions, which changes the colour scale. To prevent this we remove the outer pixels of the predictions.

# Extract the first sample from the output tensor, detach it from the computational graph,
# move it to the CPU, and convert it to a NumPy array for further processing and visualization.
output_image = landscape_predictions[0].detach().cpu().numpy()

# Create masks for the x and y coordinates, as well as a water mask (unused in the code below),
# with the same shape as the output image.
x_mask = np.ones_like(output_image)
y_mask = np.ones_like(output_image)
water_mask = np.ones_like(output_image)

# Get the dimensions of the output image.
y_dim = output_image.shape[0]
x_dim = output_image.shape[1]

# Define a buffer value to mask out edge cells.
buffer = 3

# Apply the buffer mask to the x-axis: set the first and last 'buffer' columns to -infinity.
x_mask[:, :buffer] = -np.inf
x_mask[:, x_dim - buffer:] = -np.inf

# Apply the buffer mask to the y-axis: set the first and last 'buffer' rows to -infinity.
y_mask[:buffer, :] = -np.inf
y_mask[y_dim - buffer:, :] = -np.inf

# Mask out edge cells in the output image by multiplying with the x and y masks.
# Also mask out water cells by multiplying with the water mask.
output_image = output_image * x_mask * y_mask
output_image = output_image * water_mask

# Plot the masked output image using the 'viridis' colormap.
plt.imshow(output_image, cmap='viridis')
plt.colorbar(shrink=0.7)  # Display the color scale.
plt.title(f'Log-probabilities: Day {yday_t1}, Hour {hour_t1}')
# Define the filename for saving the landscape prediction image
plt.savefig(f"{base_id_day_dir}/id{buffalo_id}_hab_log_prob_yday{yday_t1}_hour{hour_t1}.png", 
            dpi=600, bbox_inches='tight') 
plt.show()
plt.close()  # Close the figure to free up memory.

# Plot the exponential of the output image, which may be used to convert log-probabilities
# back to probability values.
plt.imshow(np.exp(output_image), cmap='viridis')
plt.colorbar(shrink=0.7)  # Display the color scale.
plt.title(f'Probabilities: Day {yday_t1}, Hour {hour_t1}')
# Define the filename for saving the landscape prediction image
plt.savefig(f"{base_id_day_dir}/id{buffalo_id}_hab_prob_yday{yday_t1}_hour{hour_t1}.png", 
            dpi=600, bbox_inches='tight') 
plt.show()
plt.close()  # Close the figure to free up memory.

Extracting prediction values

To assess how the covariates are affecting the predictions, we can extract the prediction values and correlate them with values of the landscape layers. This allows us to assess the marginal effect of each covariate on the predictions, which can look across the day and the year.

# Flatten the tensors and convert to numpy arrays
ndvi_array = ndvi_landscape
ndvi_array = ndvi_array.flatten()

canopy_array = canopy_landscape
canopy_array = canopy_array.flatten()

herby_array = herby_landscape
herby_array = herby_array.flatten()

slope_array = slope_landscape.detach().numpy()
slope_array = slope_array.flatten()

landscape_predictions_array = torch.exp(landscape_predictions).detach().cpu().numpy()
landscape_predictions_array = landscape_predictions_array * x_mask * y_mask
landscape_predictions_array = landscape_predictions_array.flatten()

# Create a single DataFrame with updated column names
df = pd.DataFrame({
    'NDVI': ndvi_array,
    'Canopy_cover': canopy_array,
    'Herbaceous_vegetation': herby_array,
    'Slope': slope_array,
    'Predictions': landscape_predictions_array
})

# Display a slice of the DataFrame
print(f"Number of cells = {len(df)}")
# Select some cells that aren't in the masked regions
print(df[15000:15010])

# Export the DataFrame to a CSV file
df.to_csv(f"{base_id_day_dir}/id{buffalo_id}_habitat_selection_landscape_subset_yday{yday_t1}_hour{hour_t1}.csv", index=False)

Number of cells = 3859188
           NDVI  Canopy_cover  Herbaceous_vegetation     Slope   Predictions
15000  0.686534       0.69697                    0.0  0.083876  4.980798e-07
15001  0.695695       0.69697                    0.0  0.078543  4.213175e-07
15002  0.702972       0.69697                    0.0  0.072087  3.567990e-07
15003  0.714315       0.69697                    0.0  0.064439  3.502799e-07
15004  0.708825       0.69697                    0.0  0.055334  3.278334e-07
15005  0.709845       0.69697                    0.0  0.045923  3.131432e-07
15006  0.714491       0.69697                    0.0  0.039663  2.740802e-07
15007  0.720424       0.69697                    0.0  0.039896  2.297743e-07
15008  0.717843       0.69697                    0.0  0.043171  1.940504e-07
15009  0.707907       0.69697                    0.0  0.043540  1.742542e-07

Plot the correlation between the predictions and the covariates

It is a bit hard to see the effects here with so many points being plotted, so we’ll do some more processing and plotting in R, so we’ll just create and export the dataframes here.

# Plot the data with plt.scatter
plt.figure(figsize=(8, 5))
plt.scatter(df['NDVI'], df['Predictions'], alpha=0.1)
plt.title('NDVI vs Predictions')
plt.xlabel('NDVI')
plt.ylabel('Predictions')
plt.grid(True)
plt.show()

plt.figure(figsize=(8, 5))
plt.scatter(df['Canopy_cover'], df['Predictions'], alpha=0.1)
plt.title('Canopy cover vs Predictions')
plt.xlabel('Canopy cover')
plt.ylabel('Predictions')
plt.grid(True)
plt.show()

plt.figure(figsize=(8, 5))
plt.scatter(df['Herbaceous_vegetation'], df['Predictions'], alpha=0.1)
plt.title('Herbaceous vegetation vs Predictions')
plt.xlabel('Herbaceous vegetation')
plt.ylabel('Predictions')
plt.grid(True)
plt.show()

plt.figure(figsize=(8, 5))
plt.scatter(df['Slope'], df['Predictions'], alpha=0.1)
plt.title('Slope vs Predictions')
plt.xlabel('Slope')
plt.ylabel('Predictions')
plt.grid(True)
plt.show()

Loop over hours

A benefit of the deepSSF approach is that it can represent temporal variation in habitat selection (and movement dynamics) across the day, which interacts with the day of the year.

To illustrate this, we can loop over the hours of the day and predict the habitat selection at each hour. We can also assess the contribution of each covariate to the predictions at each hour, giving us some idea of what the model has learned about the temporal variation in habitat selection, which can help us to further understand our species’ spatial ecology.

# To plot the prediction maps with the same colour scale, 
# we need to determine the minimum and maximum values
# Initialize landscape min and max values
landscape_vmin = float('inf')
landscape_vmax = float('-inf')

Select the day of the year

We will select a single day of the year to predict over, and then loop over the hours of the day.

Ensure that the correct NDVI layer is selected for the day of the year. - Day 225 is in August, so we will use the NDVI layer from August 2018. - Day 45 is in February, so we will use the NDVI layer from February 2019.

# Uncomment one set of covariates corresponding to the NDVI layer you want to use

yday_t1 = 225 # corresponds to the August 2018 NDVI layer
ndvi_landscape = ndvi_aug2018_subset

# yday_t1 = 45 # corresponds to the February 2019 NDVI layer
# ndvi_landscape = ndvi_feb2019_subset

# To use the full extent
# ndvi_landscape = ndvi_aug2018_norm
# ndvi_landscape = ndvi_feb2019_norm

# Convert day-of-year values into sine and cosine components.
yday_t1_sin1 = np.sin(2 * np.pi * yday_t1 / 365.25)
yday_t1_cos1 = np.cos(2 * np.pi * yday_t1 / 365.25)
# yday_t1_sin2 = np.sin(4 * np.pi * yday_t1 / 365.25)
# yday_t1_cos2 = np.cos(4 * np.pi * yday_t1 / 365.25)

# Day of the year and its sine/cosine components.
yday_t1_sin1_tensor = torch.tensor(yday_t1_sin1).float()
yday_t1_cos1_tensor = torch.tensor(yday_t1_cos1).float()
# yday_t1_sin2_tensor = torch.tensor(yday_t1_sin2).float()
# yday_t1_cos2_tensor = torch.tensor(yday_t1_cos2).float()

Update the directory with the correct day of the year.

# To save the probabilities and images, create directories
base_id_day_dir = f'{base_id_dir}/yday{yday_t1}'
os.makedirs(base_id_day_dir, exist_ok=True)
print(f"Output directory: {base_id_day_dir}")

base_id_day_dir_log = f'{base_id_dir}/yday{yday_t1}/log_probabilities'
os.makedirs(base_id_day_dir_log, exist_ok=True)
print(f"Output directory: {base_id_day_dir_log}")

base_id_day_dir_probs = f'{base_id_dir}/yday{yday_t1}/probabilities'
os.makedirs(base_id_day_dir_probs, exist_ok=True)
print(f"Output directory: {base_id_day_dir_probs}")

Output directory: outputs/model_training/djelk_derived_covs_CNN_4L8F_move_2_2025-05-26/landscape_predictions/yday225
Output directory: outputs/model_training/djelk_derived_covs_CNN_4L8F_move_2_2025-05-26/landscape_predictions/yday225/log_probabilities
Output directory: outputs/model_training/djelk_derived_covs_CNN_4L8F_move_2_2025-05-26/landscape_predictions/yday225/probabilities

As we need to run over the full loop to get the landscape_vmin and landscape_vmax values, we will first run the loop over the hours of the day once without plotting (which is much faster), and then in the following code chunk run over the loop again but plot the predictions and save the images.

This loop generates the predictions and saves the probability values as a csv to correlate with covariate values (which we do in R).

# Create a single DataFrame with updated column names
df_hourly = pd.DataFrame({
    'NDVI': ndvi_array,
    'Canopy_cover': canopy_array,
    'Herbaceous_vegetation': herby_array,
    'Slope': slope_array
})

# Define the range of hours you want to loop over
hours = range(1,25) # to start at 1

# As we used sine and cosine terms to represent the hour of the day,
# rather than the hour as integers, we can use a continuous range of hours
# (Uncomment line below)
# hours = np.arange(0,24, 0.1)

for hour_t1 in hours:

    # Convert hour-of-day values into sine and cosine components.
    hour_t1_sin1 = np.sin(2 * np.pi * hour_t1 / 24)
    hour_t1_cos1 = np.cos(2 * np.pi * hour_t1 / 24)

    # Convert numpy objects to PyTorch tensors
    hour_t1_sin1_tensor = torch.tensor(hour_t1_sin1).float()
    hour_t1_cos1_tensor = torch.tensor(hour_t1_cos1).float()

    # Stack tensors column-wise to create a tensor of shape 
    scalar_covariates = torch.stack([hour_t1_sin1_tensor, 
                                     hour_t1_cos1_tensor, 
                                     yday_t1_sin1_tensor,
                                     yday_t1_cos1_tensor], 
                                     dim=0).to(device)
    
    # Convert each scalar covariate into a grid matching the dimensions of the NDVI landscape
    scalar_grids = torch.cat([scalar_to_grid(tensor, 
                                             ndvi_landscape.shape[0], 
                                             ndvi_landscape.shape[1]) 
                                             for tensor in scalar_covariates
                                             ], dim=1).to(device)

    # Stack the spatial (landscape_stack) and scalar covariates (as grids)
    full_stack = torch.cat([landscape_stack, scalar_grids], dim=1)
    # print(full_stack.shape) # Expected shape: [1, n_channels, rows, cols]

    # Run the model
    landscape_predictions = model.conv_habitat(full_stack)
    # print(landscape_predictions.shape)

    # Pull out the prediction
    output_image = landscape_predictions[0].detach().cpu().numpy()

    # Mask out cells on the edges (that affect the colour scale)
    output_image = output_image * x_mask * y_mask

    # Check if output_image is valid before updating landscape min and max
    if output_image.size > 0:
        # Ignore masked values in the calculation
        valid_values = output_image[np.isfinite(output_image)]
        if valid_values.size > 0:
            current_min = valid_values.min()
            current_max = valid_values.max()
            # Update landscape min and max values for scaling the colour map
            landscape_vmin = min(landscape_vmin, current_min)
            landscape_vmax = max(landscape_vmax, current_max)

    predictions_array = np.exp(output_image).flatten()
    df_hourly[f'{hour_t1}'] = predictions_array

# Export the DataFrame to a CSV file
df_hourly.to_csv(f"{base_id_day_dir}/id{buffalo_id}_hourly_habitat_suitability_landscape_subset_yday{yday_t1}.csv", index=False)

print(landscape_vmin, landscape_vmax)

-34.048416 -9.004518

In this chunk we generate the plots (now the landscape_vmin and landscape_vmax values are set) and save them to the output directory.

We can save both the log-probabilities and the probabilities in their natural scale. We’ve saved the images of the log-probabilities, but just uncomment the appropriate lines to save the images of the probabilities as well.

for hour_t1 in hours:

    # Convert hour-of-day values into sine and cosine components.
    hour_t1_sin1 = np.sin(2 * np.pi * hour_t1 / 24)
    hour_t1_cos1 = np.cos(2 * np.pi * hour_t1 / 24)

    # Convert numpy objects to PyTorch tensors
    hour_t1_sin1_tensor = torch.tensor(hour_t1_sin1).float()
    hour_t1_cos1_tensor = torch.tensor(hour_t1_cos1).float()

    # Stack tensors column-wise to create a tensor of shape 
    scalar_covariates = torch.stack([hour_t1_sin1_tensor, 
                                     hour_t1_cos1_tensor, 
                                     yday_t1_sin1_tensor,
                                     yday_t1_cos1_tensor], 
                                     dim=0).to(device)
    
    # Convert each scalar covariate into a grid matching the dimensions of the NDVI landscape
    scalar_grids = torch.cat([scalar_to_grid(tensor, 
                                             ndvi_landscape.shape[0], 
                                             ndvi_landscape.shape[1]) 
                                             for tensor in scalar_covariates
                                             ], dim=1)

    # Stack the spatial (landscape_stack) and scalar covariates (as grids)
    full_stack = torch.cat([landscape_stack, scalar_grids], dim=1)
    # print(full_stack.shape) # Expected shape: [1, n_channels, rows, cols]

    # Run the model
    landscape_predictions = model.conv_habitat(full_stack)
    # print(landscape_predictions.shape)

    # Pull out the prediction
    output_image = landscape_predictions[0].detach().cpu().numpy()

    # Mask out cells on the edges (that affect the colour scale)
    output_image = output_image * x_mask * y_mask

    # Habitat selection log-probabilities
    filename_landscape_preds = f"{base_id_day_dir_log}/id{buffalo_id}_log_landscape_habitat_selection_yday{yday_t1}_hour{hour_t1}.png"
    plt.figure()  # Create a new figure
    plt.imshow(output_image) #, vmin=landscape_vmin, vmax=landscape_vmax)
    plt.colorbar(shrink=0.7)
    plt.title(f'Habitat selection (log): Day {yday_t1}, Hour {hour_t1}')
    plt.savefig(filename_landscape_preds, dpi=500) #, bbox_inches='tight'
    # plt.show()
    plt.close()  # Close the figure to free memory

    # Habitat selection probabilities
    filename_landscape_preds = f"{base_id_day_dir_probs}/id{buffalo_id}_landscape_habitat_selection_yday{yday_t1}_hour{hour_t1}.png"
    plt.figure()  # Create a new figure
    plt.imshow(np.exp(output_image)) #, vmin=np.exp(landscape_vmin), vmax=np.exp(landscape_vmax))
    plt.colorbar(shrink=0.7)
    plt.title(f'Habitat selection: Day {yday_t1}, Hour {hour_t1}')
    plt.savefig(filename_landscape_preds, dpi=500) #, bbox_inches='tight'
    # plt.show()
    plt.close()  # Close the figure to free memory

Make a GIF of the landscape predictions

First, here’s a function to call to make a gif from a given directory.

# Example sorting by the epoch number
def extract_hour(filename):
    # Extract the epoch number from the filename
    # Adjust the extraction based on your naming pattern
    import re
    match = re.search(r'hour(\d+).', filename)
    if match:
        return int(match.group(1))
    return 0

def create_gif(image_folder, output_filename, fps=10):
    """
    Creates a GIF from a sequence of images in a folder.

    Parameters:
    - image_folder: Path to the folder containing images
    - output_filename: Name of the output GIF file
    - duration: Duration of each frame in seconds
    """
    # Get all png files in the specified folder, sorted by name
    images = sorted(glob.glob(os.path.join(image_folder, '*.png')), key=extract_hour)

    # Check if any images were found
    if not images:
        print(f"No images found in {image_folder}")
        return

    # Read all images
    frames = [imageio.imread(image) for image in images]

    # Save as GIF
    imageio.mimsave(output_filename, frames, fps=fps, loop=0)

    display(Image(filename=output_filename))

    print(f"GIF created successfully: {output_filename}")

GIF of log-probabilities

# Path to your images
image_folder =  f'{base_id_day_dir_log}'
# Output GIF filename
output_filename = f'{base_id_day_dir_log}/../../landscape_predictions_log_probabilities_{yday_t1}.gif'
# Create the GIF
create_gif(image_folder, output_filename, fps=5)

<IPython.core.display.Image object>

GIF created successfully: outputs/model_training/djelk_derived_covs_CNN_4L8F_move_2_2025-05-26/landscape_predictions/yday225/log_probabilities/../../landscape_predictions_log_probabilities_225.gif

GIF of probabilities

# Path to your images
image_folder =  f'{base_id_day_dir_probs}'
# Output GIF filename
output_filename = f'{base_id_day_dir_probs}/../../landscape_predictions_probabilities_{yday_t1}.gif'
# Create the GIF
create_gif(image_folder, output_filename, fps=5)

<IPython.core.display.Image object>

GIF created successfully: outputs/model_training/djelk_derived_covs_CNN_4L8F_move_2_2025-05-26/landscape_predictions/yday225/probabilities/../../landscape_predictions_probabilities_225.gif

References

Fortin, Daniel, Hawthorne L Beyer, Mark S Boyce, Douglas W Smith, Thierry Duchesne, and Julie S Mao. 2005. “Wolves Influence Elk Movements: Behavior Shapes a Trophic Cascade in Yellowstone National Park.” Ecology 86 (5): 1320–30. https://doi.org/10.1890/04-0953.

Rhodes, Jonathan R, Clive A McAlpine, Daniel Lunney, and Hugh P Possingham. 2005. “A Spatially Explicit Habitat Selection Model Incorporating Home Range Behavior.” Ecology 86 (5): 1199–1205. https://doi.org/10.1890/04-0912.

Signer, Johannes, John Fieberg, and Tal Avgar. 2017. “Estimating Utilization Distributions from Fitted Step‐selection Functions.” Ecosphere 8 (4): e01771. https://doi.org/10.1002/ecs2.1771.