spatstat ppp objects from FIA?

[hunter-stanke]

spatstat ppp objects from FIA?

[hunter-stanke]

Follow up from my previous, how does this look for the tree-level point patterns?

library(rFIA)
library(dplyr)
library(sf)

## Using the canned RI subset
data(fiaRI)

## Get spatial information for all of most recent plot visits
plts <- fiaRI %>% 
  clipFIA() %>%
  area(byPlot = TRUE, returnSpatial = TRUE) %>%
  ## Albers equal area (meters)
  st_transform(crs = 5070)

## Get coordinates of center subplot as a dataframe
c.coords <- st_coordinates(plts) %>%
  as.data.frame() %>%
  ## Label the center SUBP appropriately
  mutate(pltID = plts$pltID, 
         SUBP = 1)

## Imagine the non-center SUBP coordinates as the vertices of an 
## equilateral triangle, where sides are of length sqrt(120^2 + 120^2) dt
## We can then compute the x and y coordinates of non-center SUBPs directly
side.length = sqrt(120^2 + 120^2) * 0.3048 # meters

## Compute coordinates of SUBP 2
t.coords <- c.coords %>%
  ## Shift X 120ft up
  mutate(Y = Y + (120 * 0.3048)) %>%
  mutate(SUBP = 2)

## Coordinates of bottom SUBPs
bl.coords <- c.coords %>%
  ## Some trig
  mutate(Y = Y - (side.length / 2),
         X = X - (side.length / 2)) %>% 
  mutate(SUBP = 3)
br.coords <- bl.coords %>%
  ## Some trig
  mutate(X = X + side.length) %>% 
  mutate(SUBP = 4)


## Combine all coordinates
coords <- bind_rows(c.coords, t.coords, bl.coords, br.coords) %>%
  rename(X.PLOT = X, Y.PLOT = Y)


## Use spatial buffers to make polygons subplots
subp <- coords %>%
  st_as_sf(coords = c('X.PLOT', 'Y.PLOT')) %>%
  st_buffer(dist = 24 * 0.3048)


## Now get the offsets in tree location from plot center
tree <- fiaRI %>% 
  clipFIA() %>%
  tpa(grpBy = c(SUBP, TREE, DIST, AZIMUTH), byPlot = TRUE) %>%
  ## Subplot only
  filter(TPA == 6.018046) %>%
  ## Convert degrees to radians
  ## Convert feet to meters
  mutate(rad = AZIMUTH * (pi/180),
         DIST = DIST * 0.3048) %>%
  mutate(x = case_when(AZIMUTH %in% c(0, 180) ~ 0,
                       AZIMUTH %in% c(90, 270) ~ DIST,
                       AZIMUTH < 90 ~ sin(rad) * DIST,
                       AZIMUTH < 180 ~ sin(pi - rad) * DIST,
                       AZIMUTH < 270 ~ -sin(rad - pi) * DIST,
                       AZIMUTH < 360 ~ -sin(2*pi - rad) * DIST),
         y = case_when(AZIMUTH %in% c(0, 180) ~ 0,
                       AZIMUTH %in% c(90, 270) ~ DIST,
                       AZIMUTH < 90 ~ cos(rad) * DIST,
                       AZIMUTH < 180 ~ -cos(pi - rad) * DIST,
                       AZIMUTH < 270 ~ -cos(rad - pi) * DIST,
                       AZIMUTH < 360 ~ cos(2*pi - rad) * DIST)) %>%
  select(pltID, 
         SUBP,
         X.TREE = x, 
         Y.TREE = y)


## Make tree coordinates absolute and make spatial
tree.coords <- coords %>%
  left_join(tree, by = c('pltID', 'SUBP')) %>%
  mutate(X = X.PLOT + X.TREE,
         Y = Y.PLOT + Y.TREE) %>%
  ## Drop non-treed forested plots
  filter(!is.na(X) | !is.na(Y)) %>%
  st_as_sf(coords = c('X', 'Y'))




## Check it out for a single plot
library(ggplot2)
ggplot() +
  geom_sf(data = filter(subp, pltID == '1_44_1_228')) +
  geom_sf(data = filter(tree.coords, pltID == '1_44_1_228')) 

[hunter-stanke]

Sorry for the (very) delayed reply here! I accidentally shut off my email notifications, and I'm just now seeing issues that were opened in August.

I haven't worked with ppp objects before, however I do have some code that will produce spatial polygons of FIA plot boundaries, and I will follow up with some code that will produce a point pattern (sf) of tree locations within plots.

For the plot boundaries, how does this look?

library(rFIA)
library(dplyr)
library(sf)

## Using the canned RI subset
data(fiaRI)

## Get spatial information for all of most recent plot visits
plts <- fiaRI %>% 
  clipFIA() %>%
  area(byPlot = TRUE, returnSpatial = TRUE) %>%
  ## Albers equal area (meters)
  st_transform(crs = 5070)

## Get coordinates of center subplot as a dataframe
c.coords <- st_coordinates(plts) %>%
  as.data.frame() %>%
  ## Label the center subplot appropriately
  mutate(plot = 1:n(), 
         subplot = 1)

## Imagine the non-center subplot coordinates as the vertices of an 
## equilateral triangle, where sides are of length sqrt(120^2 + 120^2) dt
## We can then compute the x and y coordinates of non-center subplots directly
side.length = sqrt(120^2 + 120^2) * 0.3048 # meters

## Compute coordinates of subplot 2
t.coords <- c.coords %>%
  ## Shift X 120ft up
  mutate(Y = Y + (120 * 0.3048)) %>%
  mutate(subplot = 2)

## Coordinates of bottom subplots
bl.coords <- c.coords %>%
  ## Some trig
  mutate(Y = Y - (side.length / 2),
         X = X - (side.length / 2)) %>% 
  mutate(subplot = 3)
br.coords <- bl.coords %>%
  ## Some trig
  mutate(X = X + side.length) %>% 
  mutate(subplot = 4)


## Combine all coordinates
coords <- bind_rows(c.coords, t.coords, bl.coords, br.coords) 


## Use spatial buffers to make polygons delineating micro-, sub-, and macro-plots
# subplots
subp <- coords %>%
  st_as_sf(coords = c('X', 'Y')) %>%
  st_buffer(dist = 24 * 0.3048)

# microplots
micr <- coords %>%
  mutate(X = X + (12 * 0.3048)) %>%
  st_as_sf(coords = c('X', 'Y')) %>%
  st_buffer(dist = 6.8 * 0.3048)

# macroplots
macr <- coords %>%
  st_as_sf(coords = c('X', 'Y')) %>%
  st_buffer(dist = 58.9 * 0.3048)

## Set projection appropriate prior to writing
st_crs(subp) <- 5070
st_crs(micr) <- 5070
st_crs(macr) <- 5070


## Make sure things look right
library(ggplot2)
ggplot() +
  geom_sf(data = filter(macr, plot == 1)) +
  geom_sf(data = filter(subp, plot == 1)) +
  geom_sf(data = filter(micr, plot == 1))

[jgrn307]

Hunter, thanks a ton for looking into this -- stay tuned, I think we developed (independently) a workflow working and after I confirm, I can tweak into a functional form and send to you as a potential add-on to rFIA if you'd like!

[hunter-stanke]

That would be excellent - we're always open to extensions!