# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# %matplotlib inline
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
from statsmodels.regression.linear_model import RegressionResults


# +
# turns a ascii raster file into a list of values

def raster_to_list(filename):
    data_list = []
    with open(filename, 'r') as fn:   
        for line in fn:
            try:
                float(line[0])  # see if the line represents data
                line_list = line.rstrip().split(' ')  # reformat line string
                for value in line_list:
                    data_list.append(float(value))
            except:
                continue 
    return data_list


# +
# convert the relevant rasters

model_output = raster_to_list('ABM_Output_Mean.asc')
observations = raster_to_list('cov_density.asc')
parcels = raster_to_list('parcels_1910/parcels_1910.asc')
water_distance = raster_to_list('DistToWater/disttowater.asc')

# +
# turn the raster data into a table.  Each relevant pixel is a row

data_tuples = list(zip(range(len(model_list)), model_output, observations, parcels, water_distance))
df = pd.DataFrame(data_tuples, columns=['ID', 'model', 'observation', 'parcel', 'water_distance'])

# only include pixels that were undeveloped at the start of the model run and that aren't water
# we know those won't have covenants
df = df.loc[(df['parcel'] == 0) & (df['water_distance'] > 0)]

# check the head of the dataframe to confirm it looks correcnt
df.head(20)
# -

# display a scatter plot of the relationship between the modeled output and observation data
# If the model is good, there should be a positive correlation
plt.scatter(df['model'], df['observation'])

# +
# create a linear regression model of the relationship between the predicted and observed values

# data from dataframe
X = df['model'].values.reshape(-1, 1)
y = df['observation'].values.reshape(-1, 1)

# stats model needs a constant to create an intercept
X = sm.add_constant(X)

model = sm.OLS(y, X)
results = model.fit()
print(results.summary())
# -

residualsdf = pd.DataFrame({'ID': df['ID'], 'Residual': results.resid})
residualsdf.head(20)

# +
# turn the residuals table into an ascii raster 
# Import into ArcGIS to check if they are spatially autocorrelated

with open('residuals.asc', 'w') as residuals:
    residuals.write(
        "ncols         41\nnrows         70\nxllcorner     473979.5318\nyllcorner     4970867.4064\ncellsize      250\nNODATA_value  -9999\n"
        )
    for i in range(len(model_output)):
        try:
            residuals.write(str(residualsdf.loc[i]['Residual']))
        except:
            residuals.write('0')
        if i % 41 == 40:
            residuals.write('\n')
        else:
            residuals.write(' ')

StarcoderdataJupyterScriptsDedupFiltered