Marten Space-time
Marten (Martes americana), U.P., Michigan
1 Get County Boundaries
Read county boundaries shapefile.
Counties = readOGR(dsn = "C:/Users/humph173/Documents/Michigan_State/SLP_Beam_Diam/Counties_v17a",
layer = "Counties_v18a",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\Michigan_State\SLP_Beam_Diam\Counties_v17a", layer: "Counties_v18a"
## with 83 features
## It has 15 fields
## Integer64 fields read as strings: OBJECTID FIPSNUMUP = subset(Counties, PENINSULA == "upper")
UP.union = gUnaryUnion(UP)
#Create a rasterized version
Ras = raster(res = 0.01, ext = extent(UP.union),
crs = proj4string(UP.union))
Domain.r = rasterize(UP.union, Ras,
field = 0,
background = NA)
#Creat a point grid version
Grd.pnt = rasterToPoints(Domain.r, spatial = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2]) %>%
select(-layer)2 Load Marten Harvest Data
Second data set (with age and sex)
Mart.df = read.csv("./Marten_harvests_2012_2016.csv",
header = TRUE, sep=",")
#Table with PLSS abd coordinates for Michigan
PLSS = read.csv("C:/Users/humph173/Documents/Michigan_State/PLSS/MI_PLSS.csv",
header = TRUE, sep=",")
Mart.df$Long = with(PLSS,
Long[match(
Mart.df$TRS_Taken,
TRS)])
Mart.df$Lat = with(PLSS,
Lat[match(
Mart.df$TRS_Taken,
TRS)])
Mart.df = Mart.df %>%
mutate(TakeDate = as.POSIXct(Date_Taken_Date,
tz = "America/New_York",
format = "%m/%d/%Y"),
Year = year(TakeDate),
Month = month(TakeDate),
Age = Lab_Age,
Sex = Lab_Gender,
Spp = "Marten",
OBS = 1,
Source = "MDNR") %>%
filter(is.na(Year) == FALSE,
is.na(Long) == FALSE,
is.na(Lat) == FALSE) %>%
select(TakeDate, Year, Month, Age, Sex, Spp, OBS, Year, Long, Lat, Source)
dim(Mart.df)## [1] 1608 10#Count records per year
Mart.df %>%
group_by(Year) %>%
summarise(Count = length(Year))## # A tibble: 6 x 2
## Year Count
## <dbl> <int>
## 1 2012 364
## 2 2013 361
## 3 2014 394
## 4 2015 335
## 5 2016 150
## 6 2017 4Mart.df = Mart.df %>% filter(Year < 2017)
Mart.df %>%
group_by(Year) %>%
summarise(Count = length(Year))## # A tibble: 5 x 2
## Year Count
## <dbl> <int>
## 1 2012 364
## 2 2013 361
## 3 2014 394
## 4 2015 335
## 5 2016 150dim(Mart.df)## [1] 1604 10Convert to Points
Marten.pnt = SpatialPointsDataFrame(Mart.df[,c("Long", "Lat")], Mart.df)
proj4string(Marten.pnt) = proj4string(UP)
#Keep only those in the UP
Marten.pnt$UP = is.na(over(Marten.pnt, UP.union))
Marten.pnt = subset(Marten.pnt, UP == FALSE)Quick Plot
UP_map = fortify(UP)## Regions defined for each Polygonstmp.df = Marten.pnt@data %>%
select(Long, Lat, Year)
ggplot(UP_map, aes(long,lat, group=group)) +
geom_polygon(fill = "tan", col="black") +
geom_point(data=tmp.df,
aes(Long, Lat,
group=NULL, fill=NULL),
col = "red", size=1) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Marten Harvest") +
guides(color=guide_legend(title = "Year")) +
facet_wrap(~Year, ncol = 3) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_blank(),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_text(size=20, face="bold", hjust = 0.5))
Change Projection
nProj = "+proj=utm +zone=16 +datum=NAD83 +units=km +no_defs +ellps=GRS80 +towgs84=0,0,0"
UP.prj = spTransform(UP, nProj)
UP.union.prj = spTransform(UP.union, nProj)
Marten.pnt.prj = spTransform(Marten.pnt, nProj)
Marten.pnt.prj@data = Marten.pnt.prj@data %>%
mutate(Longp = Marten.pnt.prj@coords[,1],
Latp = Marten.pnt.prj@coords[,2])3 Construct Mesh
max.edge = 2
bound.outer = 30
bdry = inla.sp2segment(UP.union.prj)
mesh = inla.mesh.2d(boundary = bdry,
loc=cbind(Marten.pnt.prj),
max.edge = c(1, 12)*max.edge,
cutoff = 5,
min.angle = 23,
offset = c(max.edge, bound.outer))
plot(mesh, lwd=0.5, main=" ",
rgl=F, draw.vertices = TRUE,
vertex.color = "blue")
mesh$n## [1] 1860Get Mesh Nodes
dd = as.data.frame(cbind(mesh$loc[,1],
mesh$loc[,2]))
names(dd) = c("Longp", "Latp")
dd = dd %>%
mutate(OBS = 0,
Source = "Mesh")
##Make a copy for each year
Year.lvls = levels(factor(Marten.pnt.prj$Year))
for(i in 1:length(Year.lvls)){
tmp.dd = dd %>%
mutate(Year = Year.lvls[i])
if(i==1){dd.df = tmp.dd}
else{dd.df = rbind(dd.df, tmp.dd)}}
Month.lvls = levels(factor(Marten.pnt.prj$Month))
Age.lvls = levels(factor(Marten.pnt.prj$Age))
dd.length = dim(dd.df)[1]
#Random values for background points
dd.df$Month = sample(Month.lvls, dd.length, replace = T)
dd.df$Sex = sample(c("F","M"), dd.length, replace = T)
dd.df$Age = sample(Age.lvls, dd.length, replace = T)
levels(factor(dd.df$Year))## [1] "2012" "2013" "2014" "2015" "2016"head(dd.df)## Longp Latp OBS Source Year Month Sex Age
## 1 522.7310 5050.779 0 Mesh 2012 6 M 1.5
## 2 700.1876 5066.286 0 Mesh 2012 10 M 3.5
## 3 489.1956 5057.737 0 Mesh 2012 6 M 0.5
## 4 539.7490 5065.947 0 Mesh 2012 10 F 3.5
## 5 512.3788 5064.130 0 Mesh 2012 12 M 8.5
## 6 526.6298 5064.886 0 Mesh 2012 6 F 8.5Combine Marten and Nodes Also converting to points.
Mod.pnts = Marten.pnt.prj@data %>%
select(Longp, Latp, OBS, Source, Year, Month, Sex, Age)
Mod.pnts = rbind(Mod.pnts, dd.df)
Mod.pnts = SpatialPointsDataFrame(Mod.pnts[,c("Longp", "Latp")], Mod.pnts)
proj4string(Mod.pnts) = proj4string(Marten.pnt.prj)4 Monthly PRISM Data
This will download PRISM data to specified folder and then convert it to a raster. Matching climate to month of harvest.
Max Temperature (TMAX)
library(prism)
range(Marten.pnt$Year)## [1] 2012 2016range(Marten.pnt$Month)## [1] 1 12options(prism.path = "./PRISM/TMAX") #Where to save the data; seperate folder for temp
get_prism_monthlys(type = 'tmax', #climate variable wanted, options: "ppt", "tmean", "tmin", "tmax"
years=2012:2016, #Years wanted
mon = 1:12, #Months wanted
keepZip = FALSE) #Don't keep zip file##
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|=== | 5%
|
|==== | 7%
|
|===== | 8%
|
|====== | 10%
|
|======== | 12%
|
|========= | 13%
|
|========== | 15%
|
|=========== | 17%
|
|============ | 18%
|
|============= | 20%
|
|============== | 22%
|
|=============== | 23%
|
|================ | 25%
|
|================= | 27%
|
|================== | 28%
|
|==================== | 30%
|
|===================== | 32%
|
|====================== | 33%
|
|======================= | 35%
|
|======================== | 37%
|
|========================= | 38%
|
|========================== | 40%
|
|=========================== | 42%
|
|============================ | 43%
|
|============================= | 45%
|
|============================== | 47%
|
|=============================== | 48%
|
|================================ | 50%
|
|================================== | 52%
|
|=================================== | 53%
|
|==================================== | 55%
|
|===================================== | 57%
|
|====================================== | 58%
|
|======================================= | 60%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================== | 65%
|
|=========================================== | 67%
|
|============================================ | 68%
|
|============================================== | 70%
|
|=============================================== | 72%
|
|================================================ | 73%
|
|================================================= | 75%
|
|================================================== | 77%
|
|=================================================== | 78%
|
|==================================================== | 80%
|
|===================================================== | 82%
|
|====================================================== | 83%
|
|======================================================= | 85%
|
|======================================================== | 87%
|
|========================================================= | 88%
|
|========================================================== | 90%
|
|============================================================ | 92%
|
|============================================================= | 93%
|
|============================================================== | 95%
|
|=============================================================== | 97%
|
|================================================================ | 98%
|
|=================================================================| 100%GetFiles = ls_prism_data(name=TRUE) #File names
#Count files
dim(GetFiles)## [1] 60 2#View (2 columns)
head(GetFiles)## files
## 1 PRISM_tmax_stable_4kmM2_201201_bil
## 2 PRISM_tmax_stable_4kmM2_201202_bil
## 3 PRISM_tmax_stable_4kmM2_201203_bil
## 4 PRISM_tmax_stable_4kmM2_201204_bil
## 5 PRISM_tmax_stable_4kmM2_201205_bil
## 6 PRISM_tmax_stable_4kmM2_201206_bil
## product_name
## 1 Jan 2012 - 4km resolution - Maximum temperature
## 2 Feb 2012 - 4km resolution - Maximum temperature
## 3 Mar 2012 - 4km resolution - Maximum temperature
## 4 Apr 2012 - 4km resolution - Maximum temperature
## 5 May 2012 - 4km resolution - Maximum temperature
## 6 Jun 2012 - 4km resolution - Maximum temperatureMyMonths = substr(GetFiles$product_name, 1, 3) #Pull month abbreviation
MyYear = substr(GetFiles$product_name, 6, 9) #Pull year
#stack rasters
TMAX.stk = prism_stack(ls_prism_data())
#Simplified names
names(TMAX.stk) = paste(MyMonths, MyYear, sep="_")
#Plot one
plot(TMAX.stk[[1]])
##Create TMAX Dataframe
TMAX.df = as.data.frame(
extract(TMAX.stk, #Finding values for each location
spTransform(Mod.pnts,
proj4string(TMAX.stk)),
method = "simple"))
TMAX.df$RowID = 1:nrow(TMAX.df)
dim(TMAX.df)## [1] 10884 61head(TMAX.df)## Jan_2012 Feb_2012 Mar_2012 Apr_2012 May_2012 Jun_2012 Jul_2012 Aug_2012
## 1 -1.49 0.97 9.06 9.30 19.04 22.39 26.54 24.23
## 2 -2.35 0.13 8.77 9.64 19.36 23.37 27.19 24.30
## 3 -1.53 0.37 8.91 11.39 20.33 23.79 27.50 25.20
## 4 -0.83 0.77 9.64 11.68 20.07 23.67 27.16 24.80
## 5 -1.28 0.54 9.18 11.28 20.04 23.88 27.58 24.85
## 6 -0.83 0.77 9.64 11.68 20.07 23.67 27.16 24.80
## Sep_2012 Oct_2012 Nov_2012 Dec_2012 Jan_2013 Feb_2013 Mar_2013 Apr_2013
## 1 18.68 11.67 5.04 0.64 -2.94 -3.65 -0.57 4.73
## 2 18.39 11.08 4.07 -0.24 -4.02 -4.92 -1.55 4.40
## 3 19.05 12.44 5.17 0.27 -2.65 -3.35 1.11 5.96
## 4 18.76 12.36 5.22 0.55 -2.15 -3.65 1.01 6.30
## 5 18.88 12.24 4.89 0.57 -2.55 -3.74 0.53 5.77
## 6 18.76 12.36 5.22 0.55 -2.15 -3.65 1.01 6.30
## May_2013 Jun_2013 Jul_2013 Aug_2013 Sep_2013 Oct_2013 Nov_2013 Dec_2013
## 1 13.99 19.89 23.40 23.37 19.21 12.45 3.90 -5.40
## 2 14.69 20.74 24.33 24.03 19.53 11.94 2.74 -6.62
## 3 18.63 22.32 25.31 24.08 19.25 12.90 4.00 -6.06
## 4 18.71 22.17 24.87 23.83 19.03 12.51 3.99 -6.01
## 5 18.34 22.47 25.28 23.95 19.24 12.40 3.87 -5.95
## 6 18.71 22.17 24.87 23.83 19.03 12.51 3.99 -6.01
## Jan_2014 Feb_2014 Mar_2014 Apr_2014 May_2014 Jun_2014 Jul_2014 Aug_2014
## 1 -8.30 -7.90 -3.45 5.91 14.83 20.46 21.86 21.59
## 2 -9.87 -8.87 -3.81 5.77 15.43 21.66 22.62 22.14
## 3 -8.28 -7.32 -2.84 6.91 15.97 23.19 22.68 23.28
## 4 -7.95 -7.40 -2.17 7.27 15.99 23.23 22.21 23.38
## 5 -8.29 -7.69 -3.29 6.96 15.63 23.45 22.25 23.13
## 6 -7.95 -7.40 -2.17 7.27 15.99 23.23 22.21 23.38
## Sep_2014 Oct_2014 Nov_2014 Dec_2014 Jan_2015 Feb_2015 Mar_2015 Apr_2015
## 1 19.25 10.34 0.68 -0.65 -5.62 -9.88 1.28 7.72
## 2 19.10 9.92 -0.48 -1.49 -6.63 -10.88 0.72 8.14
## 3 19.25 11.03 1.07 -0.67 -5.89 -9.78 0.49 8.75
## 4 19.06 10.78 1.07 -0.49 -5.86 -9.25 0.80 9.03
## 5 19.20 10.58 0.78 -0.62 -6.40 -9.97 0.03 8.56
## 6 19.06 10.78 1.07 -0.49 -5.86 -9.25 0.80 9.03
## May_2015 Jun_2015 Jul_2015 Aug_2015 Sep_2015 Oct_2015 Nov_2015 Dec_2015
## 1 16.75 19.03 24.16 22.87 22.69 12.52 8.55 3.02
## 2 17.07 20.29 25.09 23.21 22.79 12.07 7.62 1.94
## 3 18.46 21.79 25.88 23.99 22.84 12.40 8.34 3.56
## 4 17.67 21.61 25.03 23.39 22.42 12.43 8.16 3.37
## 5 17.90 21.63 25.54 23.53 22.78 12.36 8.10 3.44
## 6 17.67 21.61 25.03 23.39 22.42 12.43 8.16 3.37
## Jan_2016 Feb_2016 Mar_2016 Apr_2016 May_2016 Jun_2016 Jul_2016 Aug_2016
## 1 -3.35 -1.71 3.52 5.86 16.304 20.617 24.497 24.651
## 2 -4.25 -2.73 3.15 6.10 16.970 21.632 25.465 24.825
## 3 -3.29 -2.04 3.99 7.47 18.981 22.435 26.062 25.962
## 4 -3.07 -1.89 4.52 7.79 18.145 22.229 25.515 25.661
## 5 -3.44 -2.10 3.93 7.38 18.384 22.234 25.651 25.699
## 6 -3.07 -1.89 4.52 7.79 18.145 22.229 25.515 25.661
## Sep_2016 Oct_2016 Nov_2016 Dec_2016 RowID
## 1 20.748 13.486 9.781 -1.924 1
## 2 20.265 12.883 9.033 -2.622 2
## 3 21.933 14.196 9.043 -0.968 3
## 4 21.794 14.396 9.033 -0.696 4
## 5 21.686 14.100 9.056 -0.992 5
## 6 21.794 14.396 9.033 -0.696 6Total Precipitation (PPT)
options(prism.path = "./PRISM/PPT") #Where to save the data; seperate folder for precip
get_prism_monthlys(type = 'ppt', #climate variable wanted, options: "ppt", "tmean", "tmin", "tmax"
years=2012:2016, #Years wanted
mon = 1:12, #Months wanted
keepZip = FALSE) #Don't keep zip file##
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|=== | 5%
|
|==== | 7%
|
|===== | 8%
|
|====== | 10%
|
|======== | 12%
|
|========= | 13%
|
|========== | 15%
|
|=========== | 17%
|
|============ | 18%
|
|============= | 20%
|
|============== | 22%
|
|=============== | 23%
|
|================ | 25%
|
|================= | 27%
|
|================== | 28%
|
|==================== | 30%
|
|===================== | 32%
|
|====================== | 33%
|
|======================= | 35%
|
|======================== | 37%
|
|========================= | 38%
|
|========================== | 40%
|
|=========================== | 42%
|
|============================ | 43%
|
|============================= | 45%
|
|============================== | 47%
|
|=============================== | 48%
|
|================================ | 50%
|
|================================== | 52%
|
|=================================== | 53%
|
|==================================== | 55%
|
|===================================== | 57%
|
|====================================== | 58%
|
|======================================= | 60%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================== | 65%
|
|=========================================== | 67%
|
|============================================ | 68%
|
|============================================== | 70%
|
|=============================================== | 72%
|
|================================================ | 73%
|
|================================================= | 75%
|
|================================================== | 77%
|
|=================================================== | 78%
|
|==================================================== | 80%
|
|===================================================== | 82%
|
|====================================================== | 83%
|
|======================================================= | 85%
|
|======================================================== | 87%
|
|========================================================= | 88%
|
|========================================================== | 90%
|
|============================================================ | 92%
|
|============================================================= | 93%
|
|============================================================== | 95%
|
|=============================================================== | 97%
|
|================================================================ | 98%
|
|=================================================================| 100%GetFiles = ls_prism_data(name=TRUE) #File names
#Count files
dim(GetFiles)## [1] 60 2#View (2 columns)
head(GetFiles)## files
## 1 PRISM_ppt_stable_4kmM3_201201_bil
## 2 PRISM_ppt_stable_4kmM3_201202_bil
## 3 PRISM_ppt_stable_4kmM3_201203_bil
## 4 PRISM_ppt_stable_4kmM3_201204_bil
## 5 PRISM_ppt_stable_4kmM3_201205_bil
## 6 PRISM_ppt_stable_4kmM3_201206_bil
## product_name
## 1 Jan 2012 - 4km resolution - Precipitation
## 2 Feb 2012 - 4km resolution - Precipitation
## 3 Mar 2012 - 4km resolution - Precipitation
## 4 Apr 2012 - 4km resolution - Precipitation
## 5 May 2012 - 4km resolution - Precipitation
## 6 Jun 2012 - 4km resolution - PrecipitationMyMonths = substr(GetFiles$product_name, 1, 3) #Pull month abbreviation
MyYear = substr(GetFiles$product_name, 6, 9) #Pull year
#stack rasters
PPT.stk = prism_stack(ls_prism_data())
#Simplified names
names(PPT.stk) = paste(MyMonths, MyYear, sep="_")
#Plot one
plot(PPT.stk[[1]])
##Create TMAX Dataframe
PPT.df = as.data.frame(
extract(PPT.stk, #Finding values for each location
spTransform(Mod.pnts,
proj4string(PPT.stk)),
method = "simple"))
PPT.df$RowID = 1:nrow(PPT.df)
dim(PPT.df)## [1] 10884 61head(PPT.df)## Jan_2012 Feb_2012 Mar_2012 Apr_2012 May_2012 Jun_2012 Jul_2012 Aug_2012
## 1 66.56 14.89 41.69 35.40 57.51 105.18 93.65 61.38
## 2 55.05 16.56 52.06 37.44 63.10 89.98 90.66 70.90
## 3 59.03 16.51 62.16 25.55 30.48 138.97 68.66 74.36
## 4 72.24 29.98 71.89 29.42 43.95 168.68 80.88 49.68
## 5 72.38 28.28 63.01 29.18 30.70 143.20 68.14 62.76
## 6 72.24 29.98 71.89 29.42 43.95 168.68 80.88 49.68
## Sep_2012 Oct_2012 Nov_2012 Dec_2012 Jan_2013 Feb_2013 Mar_2013 Apr_2013
## 1 170.36 123.26 60.21 36.82 69.77 81.39 55.91 98.84
## 2 114.68 118.60 46.24 38.89 51.97 73.76 55.42 102.43
## 3 65.95 143.01 37.57 66.00 63.67 56.80 51.09 79.74
## 4 65.78 117.59 26.94 59.51 61.84 72.68 65.73 64.83
## 5 96.75 137.25 44.44 68.73 78.12 79.07 68.06 85.00
## 6 65.78 117.59 26.94 59.51 61.84 72.68 65.73 64.83
## May_2013 Jun_2013 Jul_2013 Aug_2013 Sep_2013 Oct_2013 Nov_2013 Dec_2013
## 1 60.76 86.07 146.69 88.30 51.57 83.37 125.68 53.66
## 2 61.94 124.17 130.42 56.72 49.55 81.63 104.18 58.54
## 3 57.94 88.95 123.07 124.01 112.55 83.16 189.34 61.44
## 4 49.55 57.14 164.69 102.63 91.40 75.74 169.94 61.50
## 5 60.09 77.65 120.21 129.17 104.91 69.55 163.15 58.66
## 6 49.55 57.14 164.69 102.63 91.40 75.74 169.94 61.50
## Jan_2014 Feb_2014 Mar_2014 Apr_2014 May_2014 Jun_2014 Jul_2014 Aug_2014
## 1 61.38 34.25 53.77 91.24 72.10 97.76 113.57 72.26
## 2 62.32 33.40 54.98 86.45 70.50 79.02 102.50 102.59
## 3 60.53 27.35 29.78 115.61 68.62 78.39 61.12 87.80
## 4 55.49 22.87 32.49 95.68 78.67 68.95 50.14 67.41
## 5 71.92 22.27 33.05 96.90 90.54 62.62 63.14 83.61
## 6 55.49 22.87 32.49 95.68 78.67 68.95 50.14 67.41
## Sep_2014 Oct_2014 Nov_2014 Dec_2014 Jan_2015 Feb_2015 Mar_2015 Apr_2015
## 1 108.23 175.05 140.96 60.68 42.45 45.24 32.19 66.88
## 2 135.00 166.09 119.86 59.93 27.91 34.68 34.36 71.04
## 3 99.50 211.95 111.07 57.55 46.10 28.00 38.94 47.69
## 4 88.62 232.13 123.18 64.71 44.46 25.94 35.94 44.77
## 5 94.17 187.67 130.72 73.91 58.56 36.80 42.97 47.61
## 6 88.62 232.13 123.18 64.71 44.46 25.94 35.94 44.77
## May_2015 Jun_2015 Jul_2015 Aug_2015 Sep_2015 Oct_2015 Nov_2015 Dec_2015
## 1 182.09 63.06 54.41 54.44 60.39 94.48 102.72 127.972
## 2 144.77 71.62 57.45 48.49 52.29 68.64 88.36 128.732
## 3 82.52 43.86 66.56 69.25 72.38 84.69 128.65 144.784
## 4 108.17 57.58 58.50 54.90 91.95 92.50 136.87 149.570
## 5 101.68 47.55 122.67 41.89 81.39 80.59 121.07 133.616
## 6 108.17 57.58 58.50 54.90 91.95 92.50 136.87 149.570
## Jan_2016 Feb_2016 Mar_2016 Apr_2016 May_2016 Jun_2016 Jul_2016 Aug_2016
## 1 102.199 69.320 65.95 77.087 69.235 71.408 110.315 89.843
## 2 65.292 64.755 82.40 90.865 87.082 72.342 79.503 110.897
## 3 46.178 49.873 85.36 91.666 28.937 34.329 80.746 121.108
## 4 44.824 50.466 76.20 89.356 32.561 54.835 101.914 86.358
## 5 75.469 59.471 68.83 85.030 27.396 67.655 107.944 81.389
## 6 44.824 50.466 76.20 89.356 32.561 54.835 101.914 86.358
## Sep_2016 Oct_2016 Nov_2016 Dec_2016 RowID
## 1 131.943 171.324 49.199 89.384 1
## 2 143.622 134.064 47.813 62.812 2
## 3 97.688 100.867 57.749 98.152 3
## 4 126.854 92.184 55.122 72.663 4
## 5 164.984 88.242 47.563 84.035 5
## 6 126.854 92.184 55.122 72.663 6Match Climate to Marten by Month
Mod.pnts$Month_ab = month.abb[as.integer(Mod.pnts$Month)] #Month abbreviation
Mod.pnts$YrMonth = paste(Mod.pnts$Month_ab, Mod.pnts$Year, sep = "_") #label matches climate columns
Mod.pnts$RowID = 1:nrow(Mod.pnts@data)
TMAX.mlt = melt(TMAX.df, "RowID")
TMAX.mlt$YrMonth = as.character(TMAX.mlt$variable)
#For marten where we know the month
tmp.frm = merge(Mod.pnts@data, TMAX.mlt,
by.x = c("RowID", "YrMonth"),
by.y = c("RowID", "YrMonth"),
all.x = TRUE, all.y = FALSE)
#For background pounts, we dont know the month
#pick random month
Annual.TMAX = as.data.frame(
tmp.frm %>%
group_by(Year) %>%
summarise(Average = mean(value, na.rm=T)))
tmp.frm$Average = with(Annual.TMAX,
Average[match(
tmp.frm$Year,
Year)])
#If marten use the month value, if not, use average
Mod.pnts$TMAX = ifelse(Mod.pnts$OBS == 1, tmp.frm$value, tmp.frm$Average)
Mod.pnts$TMAX[is.na(Mod.pnts$TMAX)] = mean(Mod.pnts$TMAX, na.rm=T)
#################################
####AS ABOVE FOR PPT
PPT.mlt = melt(PPT.df, "RowID")
PPT.mlt$YrMonth = PPT.mlt$variable
#For marten where we know the month
tmp.frm = merge(Mod.pnts@data, PPT.mlt,
by.x = c("RowID", "YrMonth"),
by.y = c("RowID", "YrMonth"),
all.x = TRUE, all.y = FALSE)
#For background pounts, we dont know the month
#Find annual averages
Annual.PPT = as.data.frame(
tmp.frm %>%
group_by(Year) %>%
summarise(Average = mean(value, na.rm=T)))
tmp.frm$Average = with(Annual.PPT,
Average[match(
tmp.frm$Year,
Year)])
#If marten use the month value, if not, use average
Mod.pnts$PPT = ifelse(Mod.pnts$OBS == 1, tmp.frm$value, tmp.frm$Average)
Mod.pnts$PPT[is.na(Mod.pnts$PPT)] = mean(Mod.pnts$PPT, na.rm=T)5 Distance to Roads
suppressMessages(library(spatstat))Rds = readOGR(dsn = "C:/Users/humph173/Documents/Michigan_State/StateRds",
layer = "Roads",
stringsAsFactors = FALSE)
Rds = spTransform(Rds, nProj) #convert projection
Road.psp = as.psp(Rds) #format roadsM.pp = as.ppp(Mod.pnts) #format marten points
#Distance from each point to nearest road
Mod.pnts$nRoad = nncross(M.pp, Road.psp)[,"dist"]
range(Mod.pnts$nRoad)## [1] 9.272237e-08 9.910993e+016 LandFire
LandFire are factor-type attributes, not continuous.
Veg.dist = raster("C:/Users/humph173/Documents/LandFire/US_VegDIST2014_10262017/Grid/us_vdist2014")
#Extract
Mod.pnts$vDist = extract(Veg.dist,
spTransform(Mod.pnts,
proj4string(Veg.dist), method = "simple"))
#Negative values indicate missing
Mod.pnts$vDist[Mod.pnts$vDist < 0] = NA
Canopy.Cvr = raster("C:/Users/humph173/Documents/LandFire/US_140CC_20180620/Grid/us_140cc")
Mod.pnts$CC = extract(Canopy.Cvr,
spTransform(Mod.pnts,
proj4string(Canopy.Cvr), method = "simple"))
Mod.pnts$CC[Mod.pnts$CC < 0] = NA
Bulk.dens = raster("C:/Users/humph173/Documents/LandFire/US_140CBD_20180620/Grid/us_140cbd")
Mod.pnts$CBD = extract(Bulk.dens,
spTransform(Mod.pnts,
proj4string(Bulk.dens), method = "simple"))
Mod.pnts$CBD[Mod.pnts$CBD < 0] = NA
Veg.type = raster("C:/Users/humph173/Documents/LandFire/US_140EVC_20180618/Grid/us_140evc")
Mod.pnts$EVC = extract(Veg.type,
spTransform(Mod.pnts,
proj4string(Veg.type), method = "simple"))
Mod.pnts$EVC[Mod.pnts$EVC < 0] = NA7 Data Scaling
Tmax.sc = scale(Mod.pnts$TMAX, scale=T, center=T)
Mod.pnts$TMAX.sc = as.numeric(Tmax.sc)
Ppt.sc = scale(Mod.pnts$PPT, scale=T, center=T)
Mod.pnts$PPT.sc = as.numeric(Ppt.sc)
Mod.pnts$nRoad10 = Mod.pnts$nRoad/108 Construct Model
Select Time Steps This model using an autoregressive term for temporal correlation.
k = length(levels(factor(Mod.pnts$Year)))
#Change years to integers
Mod.pnts$Step = as.integer(as.factor(Mod.pnts$Year))
range(Mod.pnts$Step)## [1] 1 5Spatial Temporal Indexing Relating points to the mesh.
locs = cbind(Mod.pnts$Longp, Mod.pnts$Latp)
A.mart = inla.spde.make.A(mesh,
alpha = 2,
loc=locs,
group = Mod.pnts$Step)Set Spatial Prior
FE.df = Mod.pnts@data
spde0 = inla.spde2.pcmatern(mesh,
alpha = 2,
prior.range=c(500, 0.01),
prior.sigma=c(1, 0.01),
constr = TRUE)
field1 = inla.spde.make.index("field1",
spde0$n.spde,
n.group=k) Construct Data Stack Organizing data as a list object instead of a dataframe for model fitting.
FE.lst = list(c(field1, #Spatial index (created in step above)
list(intercept1 = 1)), #intercept
list(TMAX = FE.df[,"TMAX.sc"], #Max Temp covariate
PPT = FE.df[,"PPT.sc"], #Precip covariate
nRoad = FE.df[,"nRoad10"], #nearest road
Age = as.numeric(FE.df[,"Age"]), #nearest road
vDist = FE.df[,"vDist"], #Veg disturbance
CC = FE.df[,"CC"], #canopy cover
CBD = FE.df[,"CBD"], #bulk density
EVC = FE.df[,"EVC"], #cover type
Year = FE.df[,"Year"],
Sex = FE.df[,"Sex"],
Month = FE.df[,"Month"]))
Stack1 = inla.stack(data = list(PA = FE.df$OBS), #PA = either a 1 (marten) or zero (mesh point/node)
A = list(A.mart, 1), #created 2 steps above, loaction info
effects = FE.lst,
tag = "mart.0") #just an arbitrary name for the file
#Save data to run in HPCC
#save(list=c("Stack1", "spde0"), file="./HPC/Mart_060319.RData") #File for HPCCNon Spatial Model
h.spec = list(theta=list(prior='pccor1', param=c(0, 0.9)))
YearPrior = list(prec = list(prior="pc.prec", param = c(1, 0.0001)))
iid.prior = list(theta=list(prior = "normal", param=c(0, 5))) #LandFire prior
Frm0 = PA ~ -1 + intercept1 +
#f(field1,
# model=spde,
# group = field1.group,
# control.group=list(model="ar1"),
# hyper=h.spec) +
f(Month,
model="iid",
hyper = iid.prior) +
f(vDist,
model="iid",
hyper = iid.prior) +
f(CC,
model="iid",
hyper = iid.prior) +
f(CBD,
model="iid",
hyper = iid.prior) +
f(EVC,
model="iid",
hyper = iid.prior) +
f(Year,
model="iid",
hyper = iid.prior) +
f(Sex,
model="iid",
hyper = iid.prior) +
TMAX + PPT + Age + nRoad
Model0 = inla(Frm0,
data = inla.stack.data(Stack1, spde=spde0),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.01,
prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(Stack1),
compute = TRUE,
link = 1),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model0)##
## Call:
## c("inla(formula = Frm0, family = \"binomial\", data = inla.stack.data(Stack1, ", " spde = spde0), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack1), ", " compute = TRUE, link = 1), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.fixed = list(prec = 0.01, prec.intercept = 0.001))" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.0671 53.2901 1.9070 57.2642
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -2.6905 1.1288 -4.9074 -2.6903 -0.4767 -2.6899 0
## TMAX -0.1440 0.0355 -0.2134 -0.1441 -0.0740 -0.1443 0
## PPT 0.3087 0.0365 0.2378 0.3084 0.3810 0.3079 0
## Age -0.5464 0.0214 -0.5892 -0.5462 -0.5052 -0.5456 0
## nRoad -0.3468 0.0592 -0.4653 -0.3460 -0.2329 -0.3444 0
##
## Random effects:
## Name Model
## Month IID model
## vDist IID model
## CC IID model
## CBD IID model
## EVC IID model
## Year IID model
## Sex IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for Month 0.4577 0.1395 0.2398 0.4394 0.783 0.4047
## Precision for vDist 1.2738 0.5641 0.5068 1.1646 2.678 0.9745
## Precision for CC 0.9449 0.4638 0.3408 0.8480 2.117 0.6840
## Precision for CBD 1.1039 0.5209 0.4174 0.9970 2.412 0.8160
## Precision for EVC 0.6798 0.1981 0.3691 0.6538 1.142 0.6047
## Precision for Year 1.3629 0.5858 0.5515 1.2540 2.813 1.0603
## Precision for Sex 0.6469 0.2684 0.2730 0.5975 1.309 0.5099
##
## Expected number of effective parameters(std dev): 54.79(0.00)
## Number of equivalent replicates : 198.64
##
## Deviance Information Criterion (DIC) ...............: 2719.72
## Deviance Information Criterion (DIC, saturated) ....: 2719.71
## Effective number of parameters .....................: 53.57
##
## Watanabe-Akaike information criterion (WAIC) ...: 2717.36
## Effective number of parameters .................: 48.69
##
## Marginal log-Likelihood: -1440.30
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedSpatial-temporal Model
h.spec = list(theta=list(prior='pccor1', param=c(0, 0.9)))
YearPrior = list(prec = list(prior="pc.prec", param = c(1, 0.0001)))
iid.prior = list(theta=list(prior = "normal", param=c(0, 5))) #LandFire prior
Frm1 = PA ~ -1 + intercept1 +
f(field1,
model=spde,
group = field1.group,
control.group=list(model="ar1"),
hyper=h.spec) +
f(Month,
model="iid",
hyper = iid.prior) +
f(vDist,
model="iid",
hyper = iid.prior) +
f(CC,
model="iid",
hyper = iid.prior) +
f(CBD,
model="iid",
hyper = iid.prior) +
f(EVC,
model="iid",
hyper = iid.prior) +
f(Year,
model="iid",
hyper = iid.prior) +
f(Sex,
model="iid",
hyper = iid.prior) +
TMAX + PPT + Age + nRoad
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(5.019970, 1.272142, 7.275585, 2.745718)
Model1 = inla(Frm1,
data = inla.stack.data(Stack1, spde=spde0),
family = "binomial",
verbose = TRUE,
control.fixed = list(prec = 0.01,
prec.intercept = 0.001),
control.predictor = list(
A = inla.stack.A(Stack1),
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = theta0),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) summary(Model1)##
## Call:
## c("inla(formula = Frm1, family = \"binomial\", data = inla.stack.data(Stack1, ", " spde = spde0), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack1), ", " compute = TRUE, link = 1), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.fixed = list(prec = 0.01, prec.intercept = 0.001))" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.4699 1850.2694 3.0107 1855.7500
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 -4.3941 1.1845 -6.7207 -4.3938 -2.0712 -4.3931 0
## TMAX -0.1349 0.0397 -0.2125 -0.1351 -0.0565 -0.1354 0
## PPT 0.2996 0.0402 0.2214 0.2993 0.3792 0.2988 0
## Age -0.5574 0.0236 -0.6046 -0.5571 -0.5120 -0.5565 0
## nRoad -0.5794 0.1287 -0.8346 -0.5785 -0.3292 -0.5768 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
## Month IID model
## vDist IID model
## CC IID model
## CBD IID model
## EVC IID model
## Year IID model
## Sex IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field1 205.2608 45.7521 132.0238 199.4768 310.8649
## Stdev for field1 2.3852 0.4245 1.6738 2.3428 3.3373
## GroupRho for field1 0.9939 0.0064 0.9770 0.9959 0.9995
## Precision for Month 0.4467 0.1355 0.2359 0.4288 0.7634
## Precision for vDist 1.2002 0.5346 0.4741 1.0969 2.5301
## Precision for CC 1.1068 0.5239 0.4105 1.0016 2.4198
## Precision for CBD 1.1108 0.5054 0.4319 1.0108 2.3747
## Precision for EVC 0.7918 0.2392 0.4267 0.7568 1.3588
## Precision for Year 1.3878 0.6008 0.5516 1.2785 2.8672
## Precision for Sex 0.7079 0.2998 0.2917 0.6527 1.4472
## mode
## Range for field1 188.0589
## Stdev for field1 2.2570
## GroupRho for field1 0.9986
## Precision for Month 0.3948
## Precision for vDist 0.9163
## Precision for CC 0.8187
## Precision for CBD 0.8381
## Precision for EVC 0.6916
## Precision for Year 1.0798
## Precision for Sex 0.5543
##
## Expected number of effective parameters(std dev): 122.53(0.00)
## Number of equivalent replicates : 88.83
##
## Deviance Information Criterion (DIC) ...............: 2204.96
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 115.52
##
## Watanabe-Akaike information criterion (WAIC) ...: 2198.51
## Effective number of parameters .................: 102.30
##
## Marginal log-Likelihood: -1236.00
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed9 Results
Year Effect
mic.df = as.data.frame(Model1$summary.random$Year)
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Year Effect") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))
Month Effect
mic.df = as.data.frame(Model1$summary.random$Month)
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Month Effect") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))
Sex Effect
mic.df = as.data.frame(Model1$summary.random$Sex)[2:3,]
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(x=ID, y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Sex Effect") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))
Vegetation Disturbance Effect
mic.df = as.data.frame(Model1$summary.random$vDist)
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(x=factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Vegetation Disturbanve Class") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))
Canopy Cover Effect
mic.df = as.data.frame(Model1$summary.random$CC)
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(x=factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Canopy Cover Class") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))
Bulk Density Effect
mic.df = as.data.frame(Model1$summary.random$CBD)
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(x=factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Bulk Density Class") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))
Vegetative Cover Effect
mic.df = as.data.frame(Model1$summary.random$EVC)
names(mic.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
ggplot(mic.df, aes(x=factor(ID), y=Mean)) +
geom_point(size=2, pch=19, col = "red") +
geom_linerange(aes(ymin=Q025, ymax=Q975), colour="black") +
geom_hline(yintercept = 0,
linetype = "dotted",
colour = "red",
size = 1) +
theme_classic() +
xlab("Cover Type Class") +
ylab("Effect (logit)") +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=14),
axis.text.y = element_text(face="bold", size=16),
axis.text.x = element_text(face="bold", size=10,
vjust=1, hjust=1, angle=45),
strip.text = element_text(face="bold", size = 20))