MIDTERM

Markdown Author: Jessie Bell

Download this Rmd: Top right corner → Code → Download Rmd

1. PNW Map

1.1

Map precipitation using ggplot.

a: Read in pnw_norms.rds

PNW_df <- readRDS("MidtermData/data/pnw_norms.rds")

b: Use map_data to create states of interest

usa <- map_data("state")

idaho <- subset(usa, region=="idaho")

washington <- subset(usa, region=="washington")

oregon <- subset(usa, region=="oregon")

PNW <- rbind(idaho, washington, oregon)

PNW_df$Inches <- PNW_df$prcp

c: Create ggplot

p1 <- ggplot(data=PNW, aes(x=long, y=lat)) + 
  geom_polygon(data=PNW, aes(group=group), color="white", fill="gray") +
  guides(fill=FALSE)

p1+ geom_point(data=PNW_df, aes(x=longitude, y=latitude, color=Inches, size=Inches))+
  scale_color_viridis_b(alpha=0.5, guide="legend")+ #colors
  scale_size_continuous(range=c(1,6))+ #combine 2 legends
  theme(axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank(),
        panel.background=element_blank(), #remove gray grid
    legend.position = c(0.87, 0.7), #legend location
    text = element_text(family = "Corbel"), #corbel font
    plot.title = element_text(hjust = 0.5))+
    labs(title='Total Annual Precipitation (1991-2020)', y=(expression(paste("z ", symbol('\256')))), x="")

Figure 1: Map showing the 30-year averages of precipitation throughout Idaho, Oregon, and Washington State, USA using 435 stations. Annual climate normals were calculated using average precipitation in inches with data from FedData from the ROpenSci package in R.

1.2

Map temperature using tmap.

a: Read in the data (see 1.1a above)

b: Convert data frame into simple features

PNW_sf <- st_as_sf(x = PNW_df,                         
           coords = c("longitude", "latitude"),
           crs = 4326)

c: Create tmap

PNW_sf$"Average Temps (F)" <- PNW_sf$tavg


tmap_mode("view")
tm_shape(PNW_sf) +
  tm_symbols(id="name", col="Average Temps (F)", alpha=0.7, palette = rev(c("red","#ffb3c8", "#f5f5f5","#2698ce","#003666")))

Figure 2: Map showing 30-year average temperatures (°F) throughout Idaho, Oregon, and Washington State, USA using 435 stations. Highest averages occuring near Yakima, WA, Medford and Portland, OR, and at Browlee Dam in Idaho. Click on the points to see where each data point is located, or click the layers icon to change the basemap.

2. Point Pattern Analysis

2.1

Analyze the PPA using various plots, tables, statistics, and include a summary paragraph.

a: Read in hardwood_ppp.rds

hardwoods_ppp <- readRDS("MidtermData/data/hardwood_ppp.rds")

b: Convert ppp → df

# create data frame with coordinates and marks
hardwoods_df <- data.frame(x = hardwoods_ppp$x, y = hardwoods_ppp$y, Trees=hardwoods_ppp$marks)

c: Summarize hardwoods_ppp using the summary function

summary(hardwoods_ppp)

## Marked planar point pattern:  1351 points
## Average intensity 0.001552001 points per square feet
## 
## Coordinates are given to 5 decimal places
## 
## Multitype:
##          frequency proportion    intensity
## blackoak       135 0.09992598 0.0001550852
## hickory        702 0.51961510 0.0008064433
## maple          514 0.38045890 0.0005904727
## 
## Window: rectangle = [0, 933] x [0, 933] feet
## Window area = 870489 square feet
## Unit of length: 1 feet

d: Explore the data with ggplot

#make cute colors
cols <- c("#f8766d", "#00bfc4", "#94bd30")

ggplot()+
  geom_point(hardwoods_df, mapping=aes(x, y, color=Trees, pch=Trees))+
   scale_colour_manual(name = "TREES",
                      labels = c("Black oak", "Hickory", "Maple"),
                      values = cols) +   
  scale_shape_manual(name = "TREES",
                     labels = c("Black oak", "Hickory", "Maple"),
                     values = c(19, 17, 16))+
  theme_light()+
  labs(x="EASTING (FT)", y="NORTHING (FT)")

Figure 3: Showing the spatial arrangement of three tree types (n=1351) within 20 acres. Data were taken from Diggle (1883).

e: Create a density plot using persp

Tree_Density <- density(hardwoods_ppp)

#make pretty continuous colors for density plot
colramp <- colorRampPalette(cols)(100)
m<- persp(Tree_Density, colmap = colramp, theta = 0, phi = 10)

Figure 4: Density plot showing highest density of trees with lower values of x and higher values of y (in the Northwest corner of the plot). There seems to be a high density of trees along the entire eastern plot border as well. This does not tell us if there are interactions between tree types though.

f: Create density plots for each species individually

#split the treetypes
split_pp <- split(hardwoods_ppp)

dens_all <- density(split_pp)
plot(dens_all, col=cols, main="")

Figure 5: Plots of density patterns for each tree type individually. Both black oak and hickory seem to be more dense in the northern part of the plot and maple seems to be more dense in the southern half. Black oak and hickory might show clustering at some point. I imagine that, by looking at these plots, maple will be competitive with the others.

g: Look for interactions between tree types using alltypes function and then plot function

Treetypes <- alltypes(hardwoods_ppp, "L", envelope=T, verbose=F)

plot(Treetypes[1,3], main="Black Oak & Maple")

Figure 6: Plot showing the transformation of Ripley’s K for black oak and maple trees in a 20 acre plot. As predicted, maple and black oak are competing for some resource.

plot(Treetypes[1,2], main="Black Oak & Hickory")

Figure 7: Plot showing the transformation of Ripley K for black oak and hickory trees in a 20 acre plot. As predicted, these trees begin to clump at about 50 ft.

plot(Treetypes[2,3], main="Hickory & Maple")

Figure 8: Plot showing the transformation of Ripley’s K for hickory and maple trees in a 20 acre plot. They are competitive at all distances shown.

Summary

In the context of a 20-acre plot, the point pattern analysis above revealed significant insights into the spatial distribution of trees, particularly maple, black oak, and hickory species. By examining \(\hat{L}\), the transformed Ripley’s K function, it was evident that competition existed between maple and the other two tree types (Fig. 6 & 8). Beyond a distance of approximately 50 feet, a pattern of clumping emerged between black oak and hickory, indicating a tendency for these species to cluster together (Fig. 7). Conversely, prior to the 50-foot threshold, the distribution of black oak and hickory exhibited complete spatial randomness.

2.2

Based on the limited information I’ve given you, how might soil moisture vary in this plot?

Answer

I have very limited tree knowledge, so I had to look this up. Maple trees are considered, mesophytic, meaning that they require a moderate amount of water to grow. Both oak and hickory are supposidly more drought tolerant and require less water, and because of this, I think that moisture levels are highest near the maple trees on the southern end of the plot (Fig. 3 & 5).

Reference

3. Variograms & Anisotropy

3.1

Use the variogram function in the gstat library to plot an empirical variogram of the particulate matter. Describe what you see. Is there autocorrelation in the data? Output: A plot and a sentence or two.

a: Read in pm_observations.rds

pm_sf <- readRDS("MidtermData/data/pm_observations.rds")

pm_df <- as.data.frame(pm_sf, xy=T)

b: Plot the data with ggplot

xyL = data.frame(st_coordinates(st_cast(pm_sf$geometry,"MULTIPOINT")))

xyL$pm <- pm_sf$pm

#break data up 
xyL_df_mutated <- xyL %>%
  mutate(Unitless = cut(xyL$pm, breaks = 5))

ggplot(xyL_df_mutated)+
geom_point(aes(x=X, y=Y,color=Unitless, size=Unitless), alpha=0.7)+
  scale_color_discrete(labels=c('0.00', '0.25', '0.50', "0.75", "1.00"))+
  scale_size_discrete(labels=c('0.00', '0.25', '0.50', "0.75", "1.00"))+
  labs(x="EASTING (M)", y="NORTHING (M)")

Figure 9: Particulate matter (unitless) for 100 \(km^2\) site.

c: Create a variogram & summary sentence (see figure 10 caption)

pmVar <- variogram(pm~1, pm_sf, cloud=F)
plot(pmVar,pch=20,cex=1, col=c("red", "orange", "pink", "steelblue"),
     ylab=expression(Semivariance~(gamma)),
     xlab="Distance (m)", main = "Particulate Matter")

Figure 10: Variogram showing isotropy. Points are autocorrelated with a separation distance of 0-1000 m, after which there is no longer any obvious autocorrelation.

3.2

Repeat but use the alpha argument to look at directionality in the data. If you set alpha to c(0,45,90,135) you’ll be looking at axes from south to north, southwest to northwest, west to east, and southeast to northwest. Describe what you see. How has your understanding changed? Output: A plot and a sentence or two.

a: Create a directional variogram & summary sentence (see figure 11 caption)

pmVar <- variogram(pm~1, pm_sf, alpha=c(0,45,90,135))
plot(pmVar,pch=20,cex=1,col=c("red", "orange", "pink", "steelblue"),
     ylab=expression(Semivariance~(gamma)),
     xlab="Distance (m)", main = "Particulate Matter Directionality")

Figure 11: Directional variogram showing anisotropy, with all directions except 45 looking relatively similar. The 0 graph shows S to N, 45 shows SW to NW, 90 shows W to E, and 135 shows SE to NE. There seems to be spatial autocorrelation from SW to NW in the 45 plot at distances of up to 1000 m between points. After 1000 m between points there is no spatial autocorrelation.

4. Moran’s I & Scale

4.1

Make a map of the canopy height data. Show each plot (i.e., HARV 037, HARV 043, HARV 046, and HARV 047).

a: read in the data

canopyheight <- readRDS("MidtermData/data/harvSubsetCHM.rds")

b: Summarize the data deets & convert

#summary(canopyheight)

#crs(canopyheight)

xyL2 <- data.frame(st_coordinates(st_cast(canopyheight$geometry),"MULTIPOINT",  crs = st_crs(32618)))

xyL2$covertype <- canopyheight$covertype
xyL2$chm <- canopyheight$chm
xyL2$plotID <- canopyheight$plotID

c: Bring in CHM data

CHM_Harv <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")

d: Convert .tif into df

CHM_Harv_df <- as.data.frame(CHM_Harv, xy=T)

e: Plot map using ggplot

a <- ggplot() +
  geom_raster(data = CHM_Harv_df , aes(x = x, y = y, fill = HARV_chmCrop)) +
    scale_fill_viridis_c(name="M")+
  geom_sf(data = canopyheight$geometry, color=alpha("white",0.5))+
  coord_sf()+
  labs(x="EASTING (M)", y="NORTHING (M)", title="Canopy Height", subtitle = "UTM Zone 18N")+
    xlim(731400, 732000)+
  ylim(4713150, 4713800)+
  theme(axis.text.x=element_text(angle=45, hjust=1), axis.text.y=element_text(angle=45, hjust = 1))+
  annotate("text", x = 731510, y=4713355, label = "46", size=2.2)+
  annotate("text", x = 731932.5, y=4713232.5, label = "47", size=2.2)+
  annotate("text", x = 731840, y=4713715, label = "43", size=2.2)+
    annotate("text", x = 731635, y=4713300, label = "37", size=2.2)

f: Try to make what Andy made (and fail because points just are not rasters)

Harv_043 <- subset(canopyheight, plotID=="HARV_043")


b <- ggplot() +
  geom_sf(data=Harv_043$geometry, aes(color=Harv_043$chm), size=5)+
  scale_color_viridis_c(name="m")+
   xlim(c(731800, 731890))+
  ylim(c(4713675, 4713755))

g: Create polygons around point sites and clip CHM_Harv to polygons

#adds polygon to outline each plot
polygons_Harv <- canopyheight %>%
  st_as_sf(coords = c("lon", "lat")) %>%
  group_by(plotID) %>%
  summarize(geometry = st_union(geometry)) %>%
  st_convex_hull()

#convert polygons to shape file (because it is easier to work with):
#st_write(polygons_Harv, "data/Harvplots.shp", driver="Esri Shapefile")

#read in new shape file
polygons_Harv_shp  <- read_sf('data/Harvplots.shp')

#clip the raster for canopy height (the .tif) to the boudary of your new shape file
masked <- mask(CHM_Harv, polygons_Harv_shp) #creates raster layer 

#note that mask function WILL NOT KNIT if you have cache=T in your chunk. Took me a long time to get this to knit.

#make dataframe for ggplot
masked <- as.data.frame(masked, xy=T)

h: Create plots using ggplot geom_raster function with limits

#build cute circles for your final map
H_043 <- ggplot() +
  geom_raster(data = na.omit(masked), aes(x, y, fill=HARV_chmCrop), show.legend = F)+
  scale_fill_viridis_c(na.value = "white")+
  xlim(c(731800, 731890))+
  ylim(c(4713675, 4713755))+
  labs(x="", y="")+
  theme(panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    panel.grid.major.x = element_blank(), 
    axis.ticks.x = element_blank(),
               axis.text.x = element_blank(), 
    axis.ticks.y = element_blank(),
               axis.text.y = element_blank(), 
    plot.title = element_text(hjust = 0.5, face="bold", size = 15))

H_047 <- ggplot() +
  geom_raster(data = na.omit(masked), aes(x, y, fill=HARV_chmCrop), show.legend = F)+
  scale_fill_viridis_c(na.value = "white")+
  xlim(c(731890, 731975))+
  ylim(c(4713190, 4713275))+
  labs(x="", y="")+
  theme(panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    panel.grid.major.x = element_blank(), 
    axis.ticks.x = element_blank(),
               axis.text.x = element_blank(), 
    axis.ticks.y = element_blank(),
               axis.text.y = element_blank(), 
    plot.title = element_text(hjust = 0.5, face="bold", size = 15)) 



H_037 <- ggplot() +
  geom_raster(data = na.omit(masked), aes(x, y, fill=HARV_chmCrop), show.legend = F)+
  scale_fill_viridis_c(na.value = "white")+
  xlim(c(731890, 731975))+
  ylim(c(4713190, 4713275))+
  labs(x="", y="")+
  theme(panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    panel.grid.major.x = element_blank(), 
    axis.ticks.x = element_blank(),
               axis.text.x = element_blank(), 
    axis.ticks.y = element_blank(),
               axis.text.y = element_blank(), 
    plot.title = element_text(hjust = 0.5, face="bold", size = 15)) 


H_046 <- ggplot() +
  geom_raster(data = na.omit(masked), aes(x, y, fill=HARV_chmCrop), show.legend = F)+
  scale_fill_viridis_c(na.value = "white")+
  xlim(c(731472, 731555))+
  ylim(c(4713315, 4713395))+
  labs(x="", y="")+
  theme(panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    panel.grid.major.x = element_blank(), 
    axis.ticks.x = element_blank(),
               axis.text.x = element_blank(), 
    axis.ticks.y = element_blank(),
               axis.text.y = element_blank(), 
    plot.title = element_text(hjust = 0.5, face="bold", size = 15)) 

i: Plots

a

Figure 12: Showing all plots included in the autocorrelation analysis of the Harvard tree sites. Here you can see that plots 46 and 47 are closest in proximity to one another, which would make me suspect they have similar, if not identical correlograms for ther spatial autocorrelation. The legend is showing canopy height in meters. The canopy height is highest in the southern part of the map and is lowest in the northeast end of the plot, with a few more low canopy sites on the western side of the plot.

ggarrange(H_037, H_043, H_046, H_047, labels = c("A. 37", "B. 43", "C. 46", "D. 47"),
                    ncol = 2, nrow = 2)

Figure 13: Showing canopy heights of each plot above, using the same legend. Plots 37 and 43 (A & B) are both evergreen plots made up of white pine (Pinus strobus) and eastern hemlock (Tsuga canadensis). Plots 46 and 47 (C & D) are made up of deciduous broadleaf trees, predominately comprised of red oak (Quercus rubra) and red maple (Acer rubrum).

4.2

Calculate a correlogram for each plot. Describe the behavior of Moran’s I for each plot.

a: Set up the data

cor_037 <- subset(canopyheight, plotID=="HARV_037")
cor_043 <- subset(canopyheight, plotID=="HARV_043")
cor_046 <- subset(canopyheight, plotID=="HARV_046")
cor_047 <- subset(canopyheight, plotID=="HARV_047")

#037
X_037 <- st_coordinates(cor_037)[,1]
Y_037 <- st_coordinates(cor_037)[,2]
CH_037 <- cor_037$chm

#043
X_043 <- st_coordinates(cor_043)[,1]
Y_043 <- st_coordinates(cor_043)[,2]
CH_043 <- cor_043$chm

#046
X_046 <- st_coordinates(cor_046)[,1]
Y_046 <- st_coordinates(cor_046)[,2]
CH_046 <- cor_046$chm

#047
X_047 <- st_coordinates(cor_047)[,1]
Y_047 <- st_coordinates(cor_047)[,2]
CH_047 <- cor_047$chm

b: Calculate max distance

#037
(max(X_037)-min(X_037))/3

## [1] 26

(max(Y_037)-min(Y_037))/3

## [1] 26.33333

c: Plot the correlogram for continuous data & describe the plots (Fig. 14 caption)

# 037
I_037 <- spline.correlog(x=X_037, y=Y_037, z=CH_037, 
                         resamp=100, xmax=40, quiet=TRUE)


# 043
I_043 <- spline.correlog(x=X_043, y=Y_043, z=CH_043, 
                         resamp=100, xmax=40, quiet=TRUE)

# 046
I_046 <- spline.correlog(x=X_046, y=Y_046, z=CH_046, 
                         resamp=100, xmax=40, quiet=TRUE)

# 047
I_047 <- spline.correlog(x=X_047, y=Y_047, z=CH_047, 
                         resamp=100, xmax=40, quiet=TRUE)

par(mfrow=c(2, 2))
plot(I_037, main="A. 37 Moran's I")
plot(I_043, main="B. 43 Moran's I")
plot(I_046, main="A. 46 Moran's I")
plot(I_047, main="A. 47 Moran's I")

Figure 14: Moran’s I spatial correlogram of canopy height for Harvard tree plots. There seems to be strong spatial autocorrelations up to 5 m for plots A & B but this drops off rapidly after 5 m. The spatial autocorrelation is strong for plots C & D at close distances but the correlation drops much less rapidly that those for the evergreen trees (A & B). It is not surprising that the plots with similar tree structure are identical in their correlogram outputs. Plots 37 & 46 (A & C) are very close in proximity to one another and it is surprising that they are not more similar in canopy height, suggesting that tree composition can also change rapidly, speculatively due to elevation and cardinal direction of the plots.

4.3

Speculate about how spatial autocorrelation is linked to the characteristic scale of the two different forest types. Backing up any wild speculation with the literature is encouraged.

Crown shapes and canopy height of both evergreen and deciduous forests likely influence spatial autocorrelation patterns in these plots. Evergreen forests likely show stronger spatial autocorrelation at smaller scales (Fig. 14A & B) because they have denser canopies, while deciduous forests seem to only show spatial autocorrelation at larger scales because they have broader, open canopies (Fig. 14C & D), though there are likley other factors influencing this as well. This speculation is supported by each of the provided paper readings in our midterm, which all suggest that the relationship between organismal characteristics and spatial autocorrelation is very complex and influenced by a range of things including species composition, canopy structure, sex, scale, and environmental conditions.