load("C:\\Users\\mmann1123\\Desktop\\CrapFolder\\data.RData")
source("C:\\Users\\mmann1123\\Desktop\\Share\\Scripts\\fire_plotter.R")
library(raster, quietly = T)
## raster 2.0-12 (1-September-2012)
library(fields, quietly = T)
## Warning: package 'fields' was built under R version 2.15.2
## Spam version 0.29-2 (2012-08-17) is loaded. Type 'help( Spam)' or 'demo(
## spam)' for a short introduction and overview of this package. Help for
## individual functions is also obtained by adding the suffix '.spam' to the
## function name, e.g. 'help( chol.spam)'.
## Attaching package: 'spam'
## The following object(s) are masked from 'package:base':
##
## backsolve, forwardsolve
library(ggplot2, quietly = T)
## Warning: package 'ggplot2' was built under R version 2.15.2
library(maptools, quietly = T)
## Checking rgeos availability: TRUE
proj = " +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
fire_perm2 = readShapePoly("C:\\Users\\mmann1123\\Desktop\\Share\\Fire\\FRAP_FirePerim11_1\\FRAP_Fire111.shp",
proj4string = CRS(proj))
fire_perm_30 = fire_perm2[fire_perm2$YEAR_ %in% seq(1971, 2000), ]
fire_perm_25 = fire_perm2[fire_perm2$YEAR_ %in% seq(1976, 2000), ]
data1 = data.frame(acres = fire_perm_30$GIS_ACRES, year = 30)
data2 = data.frame(acres = fire_perm_25$GIS_ACRES, year = 25)
data = rbind(data1, data2)
data$year = as.factor(data$year)
ggplot(data, aes(x = acres, fill = year)) + geom_density(alpha = 0.3) + xlim(c(0,
1000)) + ggtitle("Density plot of acreage of fire polygons by 25 or 30y")
## Warning: Removed 1015 rows containing non-finite values (stat_density).
## Warning: Removed 1149 rows containing non-finite values (stat_density).
aetAave_30 = raster("B:\\AggregateEnric\\Summary_Files\\aet1971_2000_Aave.tif")
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.1, released 2012/05/15 Path to GDAL shared
## files: C:/Users/mmann1123/Documents/R/win-library/2.15/rgdal/gdal Loaded
## PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470] Path to
## PROJ.4 shared files:
## C:/Users/mmann1123/Documents/R/win-library/2.15/rgdal/proj
aetAsd_30 = raster("B:\\AggregateEnric\\Summary_Files\\aet1971_2000_Asd.tif")
aetAave_25 = raster("B:\\Aggregated1080\\Summary_Files\\aet1976_2000_Aave.tif")
aetAsd_25 = raster("B:\\Aggregated1080\\Summary_Files\\aet1976_2000_Asd.tif")
aetAave_30p = na.omit(getValues(aetAave_30))
aetAave_25p = na.omit(getValues(aetAave_25))
aetAsd_30p = na.omit(getValues(aetAsd_30))
aetAsd_25p = na.omit(getValues(aetAsd_25))
data1 = data.frame(Aaet = aetAave_30p, Aaetsd = aetAsd_30p, year = 30)
data2 = data.frame(Aaet = aetAave_25p, Aaetsd = aetAsd_25p, year = 25)
data = rbind(data1, data2)
data$year = as.factor(data$year)
ggplot(data, aes(x = Aaet, fill = year)) + geom_histogram(position = "dodge") +
ggtitle("Mean Aaet")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
## this.
ggplot(data, aes(x = Aaetsd, fill = year)) + geom_histogram(position = "dodge") +
ggtitle("SDev Aaet")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
## this.
pptAave_30 = raster("B:\\AggregateEnric\\Summary_Files\\ppt1971_2000_Aave.tif")
pptAsd_30 = raster("B:\\AggregateEnric\\Summary_Files\\ppt1971_2000_Asd.tif")
pptAave_25 = raster("B:\\Aggregated1080\\Summary_Files\\ppt1976_2000_Aave.tif")
pptAsd_25 = raster("B:\\Aggregated1080\\Summary_Files\\ppt1976_2000_Asd.tif")
pptAave_30p = na.omit(getValues(pptAave_30))
pptAave_25p = na.omit(getValues(pptAave_25))
pptAsd_30p = na.omit(getValues(pptAsd_30))
pptAsd_25p = na.omit(getValues(pptAsd_25))
data1 = data.frame(Appt = pptAave_30p, Apptsd = pptAsd_30p, year = 30)
data2 = data.frame(Appt = pptAave_25p, Apptsd = pptAsd_25p, year = 25)
data = rbind(data1, data2)
data$year = as.factor(data$year)
ggplot(data, aes(x = Appt, fill = year)) + geom_histogram(, position = "dodge") +
ggtitle("Mean Appt")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
## this.
ggplot(data, aes(x = Apptsd, fill = year)) + geom_histogram(position = "dodge") +
ggtitle("SDev Appt")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
## this.
Annual - all stats calculated for each year and the annual results averaged for 25yr period (mike)
Month - all stats calculated on the average month of the 25yr period (Meg & Bioclim)
pptAave_25 = raster("B:\\Aggregated1080\\Summary_Files\\ppt1975_1999_Aave.tif")
pptAsd_25 = raster("B:\\Aggregated1080\\Summary_Files\\ppt1975_1999_Asd.tif")
pptMave_25 = raster("B:\\Aggregated1080\\Summary_Files\\ppt1975_1999_Mave.tif")
pptMsd_25 = raster("B:\\Aggregated1080\\Summary_Files\\ppt1975_1999_Msd.tif")
pptAave_25p = na.omit(getValues(pptAave_25))
pptMave_25p = na.omit(getValues(pptMave_25))
pptAsd_25p = na.omit(getValues(pptAsd_25))
pptMsd_25p = na.omit(getValues(pptMsd_25))
dataA = data.frame(ppt = pptAave_25p, pptsd = pptAsd_25p, type = "Ann")
dataM = data.frame(ppt = pptMave_25p, pptsd = pptMsd_25p, type = "Mon")
dataAM = rbind(dataA, dataM)
dataAM$type = as.factor(dataAM$type)
ggplot(dataAM, aes(x = ppt, fill = type)) + geom_histogram(position = "dodge") +
ggtitle("Annual vs Month Mean for PPT")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
## this.
ggplot(dataAM, aes(x = pptsd, fill = type)) + geom_histogram(position = "dodge") +
ggtitle("Annual vs Month SDev for PPT")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
## this.
par(mfrow = c(3, 1))
plot(pptAsd_25, main = "annual SD")
plot(pptMsd_25, main = "monthly SD")
plot(pptAsd_25 - pptMsd_25, main = "annual minus monthly SD")
I view the outputs of flint and flint as being too noisy even on a 25-30 norm. Actual resource gradient better represented with smoothed
normal = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$aetAave, fun = sum)
normal2 = normal/mean(getValues(normal)[!is.na(getValues(normal))], rm.na = T)
logged = log(normal + 1)
logged2 = log(normal2 + 1)
plot(normal, main = "Annual Aet Ave", col = bpy.colors(40))
plot(logged, main = "log(Annual Aet Ave)", col = bpy.colors(40))
plot(normal2, main = "Demeaned", col = bpy.colors(40))
plot(logged2, main = "log(Demeaned)", col = bpy.colors(40))
I plan on using log(de-meaned data), maintains spatial distribution and minimizes noise
fire = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$fire_poly, fun = sum)
logged = raster("C:/Users/mmann1123/Desktop/Share/Fire/Fire_All_Script_Outputs\\fire_count_7599_1.tif")
notlogged = raster("C:/Users/mmann1123/Desktop/Share/Fire/Fire_All_Script_Outputs\\fire_count_7599_2.tif")
plot(fire, main = "Fire Count")
plot(logged, main = "Fire Count - logged model")
plot(notlogged, main = "Fire Count - not logged model")
library(tiff)
## Warning: package 'tiff' was built under R version 2.15.2
meg = readTIFF("C:\\Users\\mmann1123\\Desktop\\Share\\Fire\\Max Discussions\\BayDelta_FprobMeg_1080m.tif")
## image 2500 x 2500 x 0, tiles 0 x 0, bps = 8, spp = 3 (output 3), config = 1, colormap = no
mike = readTIFF("C:\\Users\\mmann1123\\Desktop\\Share\\Fire\\Max Discussions\\BayDelta_FprobMike_1080m.tif")
## image 2500 x 2500 x 0, tiles 0 x 0, bps = 8, spp = 3 (output 3), config = 1, colormap = no
mike <- as.raster(mike[, , 1:3])
meg <- as.raster(meg[, , 1:3])
plot(1:2, type = "n")
rasterImage(meg, 1, 1, 2, 2, interpolate = FALSE)
plot(1:2, type = "n")
rasterImage(mike, 1, 1, 2, 2, interpolate = FALSE)
sd = standard deviation pixel i
coef variation = standard deviation i / mean i
mm = Count of less than 1SDev below mean events for PPT/CWD pixel i
mp = Count of greater than 1SDev above mean events for PPT/CWD pixel i
pm = Count of mp events followed by mm events for pixel i
pptsd = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$pptAsd, fun = sum)
pptcv = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$pptAcov, fun = sum)
pptmmppt = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$mmppt, fun = sum)
pptmpppt = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$mpppt, fun = sum)
pptpmppt = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$pmppt, fun = sum)
plot(pptsd, main = "Std Dev PPT")
plot(pptcv, main = "Cof Var PPT")
plot(pptmmppt, main = "Count of < 1SDev Event PPT (MM)")
plot(pptmpppt, main = "Count of > 1SDev Event PPT (MP)")
plot(pptpmppt, main = "Count of > 1SDev followed by < 1SDev Event PPT")
cwdsd = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$cwdAsd, fun = sum)
cwdcv = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$cwdAcov, fun = sum)
cwdmmppt = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$mmcwd, fun = sum)
cwdmpppt = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$mpcwd, fun = sum)
cwdpmppt = rasterize(data1975_1999[, c("x", "y")], raster("B:\\Aggregated1080\\Summary_Files\\OtherVariables\\Elev.tif"),
data1975_1999$pmcwd, fun = sum)
plot(cwdsd, main = "Std Dev CWD")
plot(cwdcv, main = "Cof Var CWD")
plot(cwdmmppt, main = "Count of < 1SDev Event CWD")
plot(cwdmpppt, main = "Count of > 1SDev Event CWD")
plot(cwdpmppt, main = "Count of < 1SDev followed by > 1SDev Event CWD")
Count process: Number of fire in 25yr period including zero
Zero inflatio: affect on log odds of observing an excess zero
inner = fire_plotter(zip_in = multi_sample_zip[[1]], data_in = data1975_1999,
back_cast_data = data1950_1974)
## [1] "returns: currentplot,backcastplot, outaic,current_resp_fire,mfri"
her specification is thrown off by my estimation method
(av. coefficient vs. average outcome)
basically her coefficients aren't stable over runs
so outlier pulls av. model off track
print(multi_sample_zip[[1]])
##
## Call:
## zeroinfl(formula = equation, data = data_in_z[data_in_z$train ==
## T, ], dist = "poisson")
##
## Count model coefficients (poisson with log link):
## (Intercept) cwdAsum tmxAave I(tmxAave^2)
## -6.07234 0.00134 0.29821 -0.00646
##
## Zero-inflation model coefficients (binomial with logit link):
## (Intercept) aetAsd I(aetAsd^2) petAsd pptAsd
## 6.703568 -0.067837 0.000266 -0.016065 -0.134370
## Incorp
## 0.009280
print(paste("AIC: ", inner[[3]]))
## [1] "AIC: 4166.27"
# par(mfrow=c(2,1)) plot(inner[[1]]) plot(inner[[2]])
image.plot(inner[[5]], main = "current Mean Fire return", zlim = c(0, 250))
inner = fire_plotter(zip_in = multi_sample_zip[[2]], data_in = data1975_1999,
back_cast_data = data1950_1974)
## [1] "returns: currentplot,backcastplot, outaic,current_resp_fire,mfri"
print(multi_sample_zip[[2]])
##
## Call:
## zeroinfl(formula = equation, data = data_in_z[data_in_z$train ==
## T, ], dist = "poisson")
##
## Count model coefficients (poisson with log link):
## (Intercept) log(cwdAave + 1) I(log(cwdAave + 1)^2)
## -20.086 -3.155 0.722
## log(pptAave + 1) I(log(pptAave + 1)^2) den5x5
## 7.863 -0.776 0.524
## I(den5x5^2)
## -0.147
##
## Zero-inflation model coefficients (binomial with logit link):
## (Intercept) PubLand
## -18.26 -1.95
## log(den + 1) LightMarc
## 1.21 8.47
## log(aetAave + 1) LightMarc:log(aetAave + 1)
## 5.30 -2.53
print(paste("AIC: ", inner[[3]]))
## [1] "AIC: 4017.85"
par(mfrow = c(2, 1))
plot(inner[[1]])
plot(inner[[2]])
image.plot(inner[[5]], main = "current Mean Fire return", zlim = c(0, 250))
plot(inner[[4]], main = "current fire count")
inner = fire_plotter(zip_in = multi_sample_zip[[3]], data_in = data1975_1999,
back_cast_data = data1950_1974)
## [1] "returns: currentplot,backcastplot, outaic,current_resp_fire,mfri"
print(multi_sample_zip[[3]])
##
## Call:
## zeroinfl(formula = equation, data = data_in_z[data_in_z$train ==
## T, ], dist = "poisson")
##
## Count model coefficients (poisson with log link):
## (Intercept) log(cwdAave + 1) I(log(cwdAave + 1)^2)
## -19.9214 -2.9805 0.6894
## log(pptAave + 1) I(log(pptAave + 1)^2) den5x5
## 7.6955 -0.7762 0.5867
## I(den5x5^2) Slope
## -0.1578 0.0242
##
## Zero-inflation model coefficients (binomial with logit link):
## (Intercept) PubLand
## -18.71 -1.96
## log(den + 1) LightMarc
## 1.26 7.54
## log(aetAave + 1) LightMarc:log(aetAave + 1)
## 5.38 -2.30
print(paste("AIC: ", inner[[3]]))
## [1] "AIC: 4012.38"
par(mfrow = c(2, 1))
plot(inner[[1]])
plot(inner[[2]])
image.plot(inner[[5]], main = "current Mean Fire return", zlim = c(0, 250))
plot(inner[[4]], main = "current fire count")
inner = fire_plotter(zip_in = multi_sample_zip[[4]], data_in = data1975_1999,
back_cast_data = data1950_1974)
## [1] "returns: currentplot,backcastplot, outaic,current_resp_fire,mfri"
print(multi_sample_zip[[4]])
##
## Call:
## zeroinfl(formula = equation, data = data_in_z[data_in_z$train ==
## T, ], dist = "poisson")
##
## Count model coefficients (poisson with log link):
## (Intercept) log(cwdAave + 1) log(aetAave + 1)
## -1.9172 1.2282 -13.6761
## I(log(aetAave + 1)^2) I(log(aetAave + 1)^3) den5x5
## 5.9932 -0.7241 0.4247
## I(den5x5^2) Slope log(pmcwd + 1)
## -0.1404 0.0379 0.6997
##
## Zero-inflation model coefficients (binomial with logit link):
## (Intercept) PubLand
## -25.78 -3.09
## log(den + 1) LightMarc
## 1.21 12.46
## log(aetAave + 1) LightMarc:log(aetAave + 1)
## 7.30 -3.57
print(paste("AIC: ", inner[[3]]))
## [1] "AIC: 3995.39"
par(mfrow = c(2, 1))
plot(inner[[1]])
plot(inner[[2]])
image.plot(inner[[5]], main = "current Mean Fire return", zlim = c(0, 250))
plot(inner[[4]], main = "current fire count")
inner = fire_plotter(zip_in = multi_sample_zip[[5]], data_in = data1975_1999,
back_cast_data = data1950_1974)
## [1] "returns: currentplot,backcastplot, outaic,current_resp_fire,mfri"
print(multi_sample_zip[[5]])
##
## Call:
## zeroinfl(formula = equation, data = data_in_z[data_in_z$train ==
## T, ], dist = "poisson")
##
## Count model coefficients (poisson with log link):
## (Intercept) log(cwdAave + 1) log(aetAave + 1)
## -1.7405 1.2766 -16.8653
## I(log(aetAave + 1)^2) I(log(aetAave + 1)^3) den5x5
## 7.9019 -1.0192 0.4106
## I(den5x5^2) Slope log(pmcwd + 1)
## -0.1375 0.0381 0.6798
##
## Zero-inflation model coefficients (binomial with logit link):
## (Intercept) log(den + 1)
## -9.626 1.161
## PubLand log(aetAave + 1)
## -5.397 2.764
## log(aetAave + 1):LightMarc
## -0.167
print(paste("AIC: ", inner[[3]]))
## [1] "AIC: 4023.26"
par(mfrow = c(2, 1))
plot(inner[[1]])
plot(inner[[2]])
image.plot(inner[[5]], main = "current Mean Fire return", zlim = c(0, 250))
plot(inner[[4]], main = "current fire count")
