<jupyter_start><jupyter_text>Species Distribution Modeling -- Pitch Pine in the US Northeast <NAME> Gtech733 Final Project Draft  Introduction&emsp;Environmental conservation and management depends in large part knowing where certain species of flora or fauna occur. Biologists have been documenting and mapping the distribution of plants and animals for almost 200 years, although field samples still cannot cover the entire landscape. Additionally, anthropogenic and natural forces contribute to habitat change over time, making the full distribution of a species difficult to map and define. Recent advances in computing and statistics allows us to make use of field samples and other predictors for mapping of species habitat, using an approach called species distribution modeling (SDM).    &emsp;Species distribution modeling uses environmental data to predict the distribution of species across geographic space. SDMs originated as correlative models for ecologial niche modeling, which uses the presence of species at a geographic location, as well as environmental data for that location, to build multiple regression models. A species’ distribution is commonly characterized by either occurrence probability or population density. Logistic regression is widely used for SDM, since logistic regression can be used to predict categorical variables, like presence or absence. The resulting occurrence probability map predicts areas of high or low habitat suitability for the species being modeled, and can inform future research of these locations.Species distribution modeling is based on the "fundamental niche" concept, or the range of environments where the species __could__ be found. - The fundamental niche is the range of environments within which a species can maintain a population that reproduces frequently enough that the population does not shrink. - The realized niche is the range of environments where the species is actually found. The realized niche is a subset of the fundamental niche.    &emsp;Using occurrence data from the Global Biodiversity Information Facility (GBIF) and environmental data from the WorldClim database, the purpose of this project is to model areas of predicted suitability for a species of pine tree called the Pitch Pine (Pinus Rigada) in the North East, U.S. Pitch pine is one of the species of pine that grows in the north eastern coastal pine forests, and has been under increasing strain from invasive species, habitat destruction, and other human-induced causes. For example, research shows that lower frequency of fires in pine forests negatively impacts their growth and reproduction rates. It is widely documented that it also puts them at risk of pine beetle outbreaks, which may accelerate the conversion of pitch pine-dominated forests to mixed oak-hardwood forests (Ledig, et al 1974). Modeling areas of pitch pine suitability, or the range of environements where the species could be found, will help inform management and conservation efforts.<jupyter_code>%%html
<img src="./img/realNiche.jpg", width = 600, height = 600><jupyter_output><empty_output><jupyter_text>Regression Modeling in SDM&emsp;Logsitic regression models are widely used in SDM, since the formula allows the relationship between environmental variables and presence/absence data to be modeled. Recent advances in computing and statistics allows for more complex modeling relationships to be implemented.&emsp;Logistic regression is different from linear regression, because it can predict categorical variables, where in simple linear regression the dependent variable must be continuous. In SDM, the dependent variable is binary (only two categories). The absence/presence data is usually coded as 0 for absence and 1 for presence. The logistic regression equation is used to predict the habitat suitability of the landscape, based on envionment variables at those locations where the species is observed.__The Logistic Curve Formula__P = 1 / (1 + e$^$-(a+bx)) Machine Learning Models&emsp;Machine learning models allows for much more complex relationships between the independent and dependent variables. Machine learning in logistic regression modeling uses the same logsitic curve formula equation, however, it includes more complicated equations to relate the independent variables to the dependent variable. This is done through a form of "data mining," where the computer looks for complex associations in the data that usually go unnoticed. In SDM, this usually involves applying the logistic regression formula to localized subsets of data.  1. Loading and Cleaning Data&emsp;The data used for this research comes from the Global Biodiversity Information Facility (GBIF) and the WorldClim database. GBIF aggregates data from various databases on biodiversity, with an emphasis on species occurrences (locations). The WorldClim database has high spatial resolution global weather and climate data, and this study will use the data from 1970 through 2000. Also, the USDA soil survey data will be used to include soil characteristics. The spatial scale for environmental variables used in this study will be 30 seconds (~1 km2).&emsp;The dependent variable used in this model is occurence data for a species of pine tree called pinus rigada, commonly known as pitch pine. This data has geographic coordinates (WGS84) associated with each observation, which is necessary for SDM.__The independent variables used in this study are:__ - BIO1 = Annual Mean Temperature - BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp)) - BIO3 = Isothermality (BIO2/BIO7) (×100) - BIO4 = Temperature Seasonality (standard deviation ×100) - BIO15 = Precipitation Seasonality (Coefficient of Variation) - BIO_ELEV = SRTM elevation data (used to derive bioclimatic variables in worldClim)<jupyter_code># Load all the necessary libraries
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import shapely
import pyproj
import gdal
import os
import rasterio
import sklearn

%matplotlib inline
# for ineractive plots, use:
# %matplotlib qt 
## First, load and clip the images to the aoi shapefile

# set input and output paths for clips
inputPath = './data/bioclim/'
outputPath = './data/output/'

# Get a list of all the raster files in bioclim
bioList = [file for file in os.listdir(inputPath)]

# Load the shapefile AOI
aoi = './data/aoi.shp'

bioList
## For each bioclim raster, clip to AOI with the GDAL warp function, and name it properly
for bio in bioList:
    print(outputPath + bio[10:-4]+'_clip'+bio[-4:]) #file names
    options = gdal.WarpOptions(cutlineDSName=aoi,cropToCutline=True)
    outBioVar = gdal.Warp(srcDSOrSrcDSTab=inputPath + bio,
                        destNameOrDestDS=outputPath + bio[10:-4]+'_clip'+bio[-4:],
                        options=options)
    outBioVar= None
## Clip the pitch pine points using the aoi shapefile
# Read in the aoi shapefile and pitch pine occurence points into geopandas dataframes
aoi_gdf = gpd.read_file(aoi)

pitch_pine_fp = './ppOccurences.shp'
pp_gdf = gpd.read_file(pitch_pine_fp)

print('AOI crs: ' + aoi_gdf.crs['init'])
print('Pitch Pine layer crs: ' + pp_gdf.crs['init'])

pp_gdf
## Add a new column called "class" indicating presence or absence. Make all equal to 1 (since this is presence data)
pp_gdf.insert(2, "class", 1, True)
pp_gdf
## Clip the pp_gdf with the aoi shapefile
# Fix error before draft
# use:
# pp_gdf.clip(aoi_gdf)

# Geopandas clip is not working, may be due to older Geopandas version
# After updating geopandas, this function still does not work
pp_clip = gpd.clip(pp_gdf, aoi_gdf)
## Load in clipped pitch pine occurence points
fp = './ppOccurences_clip.shp'
pp_clipped = gpd.read_file(fp)
pp_clipped
# Rename the "type" column to "class", indicating the presence/absence of an observation
# this is due to "type" being a reserved word in python, and can cause complications if we have a column or variable
# with the same name

pp_clipped.rename(columns={'type':'class'}, inplace=True)
pp_clipped
# Get a list of all the raster files in bioclim

bioClim_path = './data/output'
bioClim_list = [file for file in os.listdir(bioClim_path)]

bioClim_list
## Create a function to generate random points within the aoi (for psuedo-absence points necessary in SDM)

import random
from shapely.geometry import Point

def generate_random(number, polygon):
    list_of_points = []
    minx, miny, maxx, maxy = polygon.bounds
    counter = 0
    while counter < number:
        pnt = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
        if polygon.contains(pnt):
            list_of_points.append(pnt)
            counter += 1
    return list_of_points

# Use the list of points to create a data frame
# add a column called class, and give all those values 0 (indicating absence)
# and append to the pitch pine GeoDataFrame

### Extracting raster values using the clipped points shapefile
### Draft 1

# # Read used pp_clipped points to extract raster values to points
# pp_clipped.index = range(len(pp_clipped))
# for i in range(len(pp_clipped)):
#     coords = [(x,y) for x, y in zip(pp_clipped.geometry[i][0].x, pp_clipped.geometry[i][0].y)]


# for (bio, i) in zip(bioClim_list, range(len(bioClim_list))):
#     # Open the raster and store metadata
#     src = rasterio.open(bio)
#     # add new column for bioclim raster values
#     pp_clipped.insert(-1, "raster %d" % (i), [x[0] for x in src.sample(coords)]) 
#     # Sample the raster at every point location and store values in DataFrame
# #     pp_clipped['Raster Value'] = [x[0] for x in src.sample(coords)]
### Extracting raster values using the clipped points shapefile
## HELP needed:
## Unable to extract pixel values to points using rasterio sample function()

## Extracting raster values using the original point shapefile layer

# Read in pp_gdf points extract coordinates to a list
pp_gdf.index = range(len(pp_gdf))
coords = [(x,y) for x, y in zip(pp_gdf.geometry.x, pp_gdf.geometry.y)]


for (bio, i) in zip(bioClim_list, range(len(bioClim_list))):
    # Open the raster and store metadata
    src = rasterio.open(bioClim_path + "/" + bio)
    # add new column for bioclim raster values, and include raster values sampled at coordinate points
    pp_gdf.insert(-1, "raster %d" % (i), [x[0] for x in src.sample(coords)])
## Loading the clean pitch pine data (with psuedo absence points) generated in qgis

# Load in absence/presence pitch pine data with extracted raster values
pp_fp = './ppPts_rasterValues_clean.shp'
pp_data = gpd.read_file(pp_fp)

pp_data
# Drop NA values
pp_df = pp_data[pp_data['bio11'].notna()]
pp_df = pp_df[pp_df['bio_elev'].notna()]
pp_df.shape

pp_df
# Select only variables needed for analysis (review literature)
pp_clean = pp_df[['decimalLon', 'decimalLat', 'geometry', 'bio1', 'bio2', 'bio3', 'bio4', 'bio15', 'bio_elev', 'class']]
                 
pp_clean.shape
pp_clean<jupyter_output><empty_output><jupyter_text>2. Exploratory Data Analysis<jupyter_code># Use a multi pairs plot for visual check on collinearity between independent variables

import seaborn as sns; sns.set(style="ticks", color_codes=True)

g = sns.pairplot(pp_clean.iloc[:,3:-1], kind='reg')
# Visualize the clipped rasters that are used as the explanatory (independent) variables
from PIL import Image
# First, load the data
# set input path to the clipped rasters
bioInput = './data/output/'

# Get a list of all the raster files in bioclim
bioVarList = [file for file in os.listdir(bioInput)]

# Method 1: use a for loop to loop through biolist rasters and stack the variables used in analysis

# Method 2: Open each individual raster file used in analysis and stack them together
bio1_pil = Image.open(bioInput + 'bio_1_clip.tif')
bio2_pil = Image.open(bioInput + 'bio_2_clip.tif')
bio3_pil = Image.open(bioInput + 'bio_3_clip.tif')
bio4_pil = Image.open(bioInput + 'bio_4_clip.tif')
bio15_pil = Image.open(bioInput + 'bio_15_clip.tif')
bio_elev_pil = Image.open(bioInput + 'elev_clip.tif')

# Convert each band into a numpy array
bio1_arr = np.array(bio1_pil)
bio2_arr = np.array(bio2_pil)
bio3_arr = np.array(bio3_pil)
bio4_arr = np.array(bio4_pil)
bio15_arr = np.array(bio15_pil)
bio_elev_arr = np.array(bio_elev_pil)
# Set NA values to the bioclim rasters (all NA values are -3.4e+38 for bioclim variables, while elevation is -32768)

# set the values to NA for bioclim variables
bio1_arr[bio1_arr < 0] = np.nan
bio2_arr[bio2_arr < 0] = np.nan
bio3_arr[bio3_arr < 0] = np.nan
bio4_arr[bio4_arr < 0] = np.nan
bio15_arr[bio15_arr < 0] = np.nan

# set na values for the elevation raster
bio_elev_arr = bio_elev_arr.astype('float')
bio_elev_arr[bio_elev_arr == (-32768)] = np.nan


# Display image variables
# use the following magic line to render interactive plots to a window output
%matplotlib qt

# Dont forget to switch back to normal magic afterwards!
# %matplotlib inline

fig=plt.figure(figsize=(24,12))

# display bio1
a = fig.add_subplot(2,3,1)
bio1_plot = plt.imshow(bio1_arr, cmap="BuPu")
a.set_title('Bio1')

# display bio2
a = fig.add_subplot(2,3,2)
nirrg_plot = plt.imshow(bio2_arr, cmap = 'YlOrRd')
a.set_title('Bio2')

# display bio3
a = fig.add_subplot(2,3,3)
nirrg_plot = plt.imshow(bio3_arr, cmap='winter')
a.set_title('Bio3')

# display bio4
a = fig.add_subplot(2,3,4)
nirrg_plot = plt.imshow(bio4_arr, cmap='cool')
a.set_title('Bio4')

# display bio15
a = fig.add_subplot(2,3,5)
nirrg_plot = plt.imshow(bio15_arr, cmap='Blues')
a.set_title('Bio15')

# display bio_elev
a = fig.add_subplot(2,3,6)
nirrg_plot = plt.imshow(bio_elev_arr, cmap="gist_earth")
a.set_title('Bio_Elev')

plt.show()
bio1_arr[800,1000]<jupyter_output><empty_output><jupyter_text>3. Training the Model  Regression Modeling in SDM&emsp;Logsitic regression models are widely used in SDM, since the formula allows the relationship between environmental variables and presence/absence data to be modeled. Recent advances in computing and statistics allows for more complex modeling relationships to be implemented.&emsp;Logistic regression is different from linear regression, because it can predict categorical variables, where in simple linear regression the dependent variable must be continuous. In SDM, the dependent variable is binary (only two categories). The absence/presence data is usually coded as 0 for absence and 1 for presence. The logistic regression equation is used to predict the habitat suitability of the landscape, based on envionment variables at those locations where the species is observed.__The Logistic Curve Formula__P = 1 / (1 + e$^$-(a+bx)) Machine Learning Models&emsp;Machine learning models allows for much more complex relationships between the independent and dependent variables. Machine learning in logistic regression modeling uses the same logsitic curve formula equation, however, it includes more complicated equations to relate the independent variables to the dependent variable. This is done through a form of "data mining," where the computer looks for complex associations in the data that usually go unnoticed. In SDM, this usually involves applying the logistic regression formula to localized subsets of data.<jupyter_code>%%html
<img src="./img/log_reg_sdm.png">
%%html
<img src="./img/log_reg_basic.png">
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

# Fit a logistic regression model using all bioclim variables

# First, split into test train data sets using the train_test_split function from scikit_learn
X = pp_clean.iloc[:,3:-1]
Y = pp_clean['class']

x_train, x_test, y_train, y_test = train_test_split(X, Y, random_state=4)
x_train.shape


# Plot the decision boundary. For that, we will assign a color to each
# point in the mesh [x_min, x_max]x[y_min, y_max].
x_min, x_max = x_train[:, 0].min() - .5, x_train[:, 0].max() + .5
y_min, y_max = x_train[:, 1].min() - .5, x_train[:, 1].max() + .5
h = .02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.figure(1, figsize=(4, 3))
plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)

# Plot also the training points
plt.scatter(X[:, 0], X[:, 1], c=Y, edgecolors='k', cmap=plt.cm.Paired)
plt.xlabel('Sepal length')
plt.ylabel('Sepal width')

plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xticks(())
plt.yticks(())

plt.show()<jupyter_output><empty_output><jupyter_text>Train the model, then predict using the x_testlogreg = LogisticRegression()logreg.fit(x_train, y_train) logreg.get_params() store the predicted values in a variabley_pred = logreg.predict(x_test)<jupyter_code># Calculate the accuracy of the model using the score method on a logistic_regression class
acc = logreg.score(x_test, y_test)
print(acc * 100)

# OR calculate the accuracy of the model using the accuracy score function from sklearn.metrics
accuracy = sklearn.metrics.accuracy_score(y_test, y_pred)
accuracy_percentage = 100 * accuracy
accuracy_percentage
# Calculate probability estimates using the predict_proba method (this will be used for the mapping portion)
prob_pred = logreg.predict_proba(x_test)
prob_pred[:,1]<jupyter_output><empty_output><jupyter_text>Predicting probability of presence<jupyter_code># Stack all rasters into one numpy array to use for analysis
bio_stack = np.dstack((bio1_arr, bio2_arr, bio3_arr, bio4_arr, bio15_arr, bio_elev_arr))

bio_stack.shape

# Create an empty target raster
result = np.zeros(bio_stack[:, :, 0].shape)
result.shape
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
%matplotlib qt
# Visualize probability on a raster grid, the same shape and size as the aoi

# What is the precision and recall of the results, comparing the reference Y to the predicted Y:
p = precision_score(y_test, y_pred, average=None)
r = recall_score(y_test, y_pred, average=None)

# Note that, with average set to None, precision and recall are calculated for each class
print("Precision: {0}, Recall: {1}".format(p, r))

# Now lets use the classifier we developed

# Create an empty target raster
result = np.zeros(bio_stack[:, :, 0].shape)

# Some housekeeping to restructure the data so it fits with scikit-learn's expected parameters. Scikit learn
# doesn't know about rasters, so we have to get the data out of the grid format.
a = bio_stack.reshape((6, result.shape[0] * result.shape[1]))
b = np.swapaxes(a, 0, 1)
## Opened raster tif files using PIL Image and the converted to numpy array, where the format is [row, col, array]


# Predict everything at once
c = logreg.predict(b)

# Put the results back into raster grid format
d = c.reshape(result.shape[0], result.shape[1])

# Now we can take a look at the results
plt.imshow(d)
plt.show(block=True)
print("All done!")
# initialize an array to add labeled pixels
result = np.zeros(bio_stack[:, :, 0].shape)
for j in range(bio_stack.shape[0]):
    print(j)
    for k in range(bio_stack.shape[1]):
        # for all band values in a pixel
        b_array = bio_stack[j,k,:]
        # add band values to list to pass to classifier prediction
        b_vals = [val for val in b_array]
#         print(b_vals)
        pred = logreg.predict([b_vals])
        
        # add label value to cell
        result[j,k] = pred[0]
plt.imshow(result)<jupyter_output><empty_output><jupyter_text>Bioclimatic variables coded as:BIO1 = Annual Mean TemperatureBIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))BIO3 = Isothermality (BIO2/BIO7) (×100)BIO4 = Temperature Seasonality (standard deviation ×100)BIO5 = Max Temperature of Warmest MonthBIO6 = Min Temperature of Coldest MonthBIO7 = Temperature Annual Range (BIO5-BIO6)BIO8 = Mean Temperature of Wettest QuarterBIO9 = Mean Temperature of Driest QuarterBIO10 = Mean Temperature of Warmest QuarterBIO11 = Mean Temperature of Coldest QuarterBIO12 = Annual PrecipitationBIO13 = Precipitation of Wettest MonthBIO14 = Precipitation of Driest MonthBIO15 = Precipitation Seasonality (Coefficient of Variation)BIO16 = Precipitation of Wettest QuarterBIO17 = Precipitation of Driest QuarterBIO18 = Precipitation of Warmest QuarterBIO19 = Precipitation of Coldest Quarterbio_elev = SRTM derived digital elevation model (30s)<jupyter_code>from scipy.stats import logistic
%matplotlib inline

fig, ax = plt.subplots(1, 1)
# Calculate the first few moments
mean, var, skew, kurt = logistic.stats(moments="mvsk")
# Display the probability density function (pdf):

x = np.linspace(logistic.ppf(0.01),
                 logistic.ppf(0.99), 100)
ax.plot(x, logistic.pdf(x),
        'r-', lw=5, alpha=0.6, label='logistic pdf')
plt.imshow(bio4_arr)
plt.plot(pp_clean.geometry)
image_array.shape
y_pred
c<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08090739053282817, 'fraction_numerical': 0.014148002493885186, 'mean_word_length': 4.983357245337159, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '340334', 'license': 'MIT', 'n_tokens_mistral': 6578, 'n_tokens_neox': 6195, 'n_words': 2513, 'path': 'SDM_PitchPineNE/.ipynb_checkpoints/SDM_PitchPine-checkpoint.ipynb', 'repo_name': 'arojas314/Geocomp', 'size': '20860'}