sf_Boundaries.qmd

---
title: " Split  Boundaries"
author: "Tom Philippi"
date:  "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    highlight: haddock
    theme: journal
  pdf_document:
    highlight: haddock
  html_notebook:
    chunk_output_type: inline
params:
  boundsshp: "d:/datasets/Boundary/nps_boundary.shp"
editor: 
  markdown: 
    wrap: 72
---

# Introduction

The Park Boundaries as delivered by the Lands Division
<https://irma.nps.gov/DataStore/Reference/Profile/2224545?lnv=True> have
several quirks and inconsistencies that make some applications
difficult.

-   In most but not all cases, park / preserve pairs are 2 features, but
    they share the same UNIT_CODE value for the park. \[Craters of the
    Moon is an exception, where the National Monument is CRMO but the
    Preserve is CRMP.\]
-   Conversely, KICA and SEQU are included as separate units, with no
    simple solution for the administrative unit SEKI. And there is
    nothing linking Craters of the Moon National Monument (CRMO) with
    Craters of the Moon National Preserve (CRMP).
-   Several widely dispersed multi-unit parks are single (multipolygon)
    features, so bounding boxes around them are huge, up to 80% of CONUS
    for MAPR )Hanford WA Los Alamos NM & Oakridge TN).
-   Long diagonal units (APPA, BLRI, NATR) also have huge bounding
    boxes. Splitting them into parts supports a series of smaller
    bounding boxes.\
-   The units of the Shape_Area variable are unspecified, and the
    "metadata" <https://irma.nps.gov/DataStore/DownloadFile/733887> are
    an empty shell.

This script takes a shapefile of park boundaries as published by the
Lands Division, and splits multi-part units into multiple features, then
produces additional Unit Code related variables for consistently split
and consistently lumped units, as well as retaining the verbatim
UNIT_CODE value used by Lands.

-   UNIT_CODE: The verbatim UNIT_CODE value from the source Lands
    shapefile.
-   SplitCode: UnitCode for the most granular distinctions (splitting
    NEPE, but not NACE nor NAMA), and giving Preserves of Park/Preserve
    separate unit code values.
-   GroupCode: Consistently lumping Park/Monument & Preserve, also SEKI,
-   UCaction: the action taken on this feature "verbatim", "split",
    "segmented"

This is a reincarnation of work that was done circa 2016-2021 by Tom
Philippi, Lisa Nelson, and Michelle Kinseth. In working with gbif and
other species occurrence data, Tom needed reasonably compact bounding
boxes for each unit. He also needed updated 3km and 30km buffers around
park boundaries as areas of analysis. Tom started with splits for MAPR,
and a stack of bounding boxes for each of the long diagonal units (APPA,
BLRI, NATR). Lisa and Michelle had partially overlapping needs for the
landscape and climate data for the National Protocol. In order to
prevent equivalent analyses from different workflows from producing
slightly different results, we agreed to standardize both buffered and
split units and push that processing upstream into the production of
IMDbounds. Micehelle and Lisa generalized the splitting of units to
include multiple islands (e.g., CHIS, NPSA) and otherwise greatly
improved and extended (as well as thoroughly tested) this concept. It is
documented in several of their reports, including:

> Kinseth M and Nelson L. 2020. Unit boundary processing for National
> Park Service units for use in the national park environmental settings
> monitoring protocol: Fall 2018. Natural Resource Data Series.
> NPS/NRSS/IMD/NRDS—2020/1253. National Park Service. Fort Collins,
> Colorado <https://irma.nps.gov/DataStore/Reference/Profile/2270597>
>
> Kinseth M and Nelson L. 2022. Boundary-derived Areas of Analysis for
> National Park Service Units, Fall 2021. Data Release Report.
> NPS/NRSS/DRR—2022/3. National Park Service.
> https://doi.org/10.36967/2287628
>
> <https://irma.nps.gov/DataStore/Reference/Profile/2287628>

That processing has been discontinued, so this version is more of a
stand-alone script for use in standardized preprocessing.

The main change is the switch from sp spatial objects and tools to sf
objects and tools. However, over the years, the boundary polygons
produced by the Lands Division have changed, and several new units have
been created that are either dispersed units or long trails.

This (2025) version does not create a new variabe UnitCode. Rather, it
creates additional separate variabes for SplitCode and GroupCode.

+----------------+----------------+----------------+
| Case           | SplitCode      | GroupCode      |
+================+================+================+
| normal         | 4 character    | 4 character    |
|                | UNIT_CODE      | UNIT_CODE      |
+----------------+----------------+----------------+
| P              | Park:          | Park UNIT_CODE |
| ark/Monument + | UNIT_CODE      | (DENA)         |
| Preserve       | (DENA)         |                |
|                |                |                |
|                | Preserve: 3    |                |
|                | char of        |                |
|                | UNIT_CODE +    |                |
|                | "P" (DENP)     |                |
+----------------+----------------+----------------+
| Disjunct Parts | UNIT_CODE +    | 4 character    |
|                | name           | UNIT_CODE      |
|                | (MAPR_Hanford) | (MAPR)         |
+----------------+----------------+----------------+
| Diagonal       | UNIT_CODE +    | 4 character    |
|                | digit (APPA_1, | UNIT_CODE      |
|                | APPA_2         | (APPA)         |
+----------------+----------------+----------------+
| Unified        | UNIT_CODE      | Standard code: |
| Management     | (KIKA or SEQU) | SEKI           |
+----------------+----------------+----------------+

This is the first iteration of this newer process. Thus, this narrative,
and some of the code, is rough around the edges. Several chunks of code
that were used in exploring the lands dataset that might be useful if
that dataset changes substantially remain in this file, but are not
executed. Rerunning this against future Lands Division Releases wil
result in a more polished document and code.

# Approach

In order to run this script, you must first download the L[ands Division
zipped shapefile from
DataStore](https://irma.nps.gov/DataStore/Reference/Profile/2224545?lnv=True),
unzip it to a directory, then set the YAML parameter boundsshp to point
to the unzipped shapefile.

There are 5 steps to this process. The first simple step is to generate
SplitCodes for the Preserve half of Park/Monument + Preserve pairs, and
copy the park's UNIT_CODE value into GroupCode. The second simple step
is to create GroupCodes for KICA + SEQU = SEKI, and for CRMO + CRMP =
CRMO.

The next two steps are more complicated, so the records are spit into 3
groups: those with widely dispersed subunits they need to be split into,
those long diagonal that need to be cut into segments for multiple
bounding boxes, and those that need no modification. In theory the
natural subunits and segments units could be a single group, as now the
actions are the same. But, they require a bit more and different
checking, so they are treated as a separate group.

Finally, after processing, those groups of units need to be combined
back into a single sf polygons object and written in several spatial
data formats.

# Code with Commentary

The first chunk is boilerplate loading of needed packages, installing
those that are not already installed. This is old school, not dependent
on any packages aready being installed. pak::pkg_install() can do this
in a single call, but that requires installing pak if it isn't already
installed, and pak occasionally has issues with circular package
dependencies.

```{r setup}
#| echo: FALSE
#| message: FALSE
#| eval: TRUE 


# "e:/boundary/nps_boundary.shp"
# "d:/boundaries/Administrative Boundaries of National Park System Units.shp"


#-- Global settings
options(timeout = 600)
gc() # Call a garbage collection to free up memory


Sys.setenv("PROJ_LIB"="")

#-- Packages
# R packages used in this script
pkgLst <- c(
            "plyr",
            "tidyverse",
            "reshape2",
            "sf",          # spatial functions
            "sfheaders",   # helper functions
            "jsonlite",    # JSON parser
            "units",       # unit transformation from area to ha and acres
 #           "terra",       # SpatRaster
 #           "tidyterra",   # for geom_spat
            "tigris",       # for political boundaries
            "ggplot2",
            "ggforce",     # for facet_wrap_paginate
            "ggspatial",
            "knitr",
            "rmarkdown"   # for render()
            )

#-- Install
# Install packages if they aren't in your library
instPkgs <- pkgLst %in% rownames(installed.packages())
if (any(instPkgs == FALSE)) {
  install.packages(pkgLst[!instPkgs], 
                   lib =  .libPaths()[1], 
                   repos = "https://cloud.r-project.org",
 #                  type = 'source', 
 #                  quiet = TRUE,
                   dependencies = TRUE) 
}

# Load packages into work space
# Note: This script is written so the packages do not need to be loaded.
#     Comment out the next line if you want to supress loading packages.
invisible(lapply(pkgLst, library, character.only = TRUE))


#-- knitr
# Set global options for RMarkdown
knitr::opts_chunk$set(eval = TRUE, 
                      echo = FALSE, 
                      results = 'hide', 
                      comment = "",
                      message=FALSE, 
                      warning=FALSE,
                      fig.path = "Figures/",
                      tidy = TRUE, 
                      tidy.opts = list(width.cutoff = 60), 
                      cache = FALSE)
 
#-- ggplot default theme updates
if ("ggplot2" %in% pkgLst) {
  #-- ggplot2 theme
  ggplot2::theme_set(ggplot2::theme_bw())
  ggplot2::theme_update(plot.title = ggplot2::element_text(hjust = 0.5),
                        plot.subtitle = ggplot2::element_text(hjust = 0.5))
}

```

This next chunk loads some of Tom's spatial functions, in this script
used to calculate area for all features. Perhaps some day Tom will make
these a utility package, but until then the latest version of them
exists in TomsSpatial.R

```{r TomsSpatialFunctions}
#| echo: FALSE
#| message: TRUE
#| eval: TRUE 

#' function to determine UTM zone for spatial object
#' 
#' This function takes an sf spatial object and returns the relevant UTM Zone.
#' It does not safeguard against being passed spatial obects so large 
#' that they extend across multiple UTM zones.  Modest-sized spatial objects
#' that cross UTM zone boundaries should end up with the lower-numbered zone,
#' which is appropriate for the main use-case of this function: supporting
#' the buff_poly() function that needs to transform a spatial object into 
#' a projected CRS with units of meters to perform the buffering.
#'
#' @param sfo A an sf spatial object, typically polygons.
#'
#' @return A character string of the utm zone, such as 12N or 8S.
#' 
#' @examples get_zone(FOLSbounds)
get_zone <- function(sfo) {
   co <- sfo |> 
          sf::st_transform(crs = 4326) |>
          sf::st_centroid() |>   # this step require to work with MULTIPOLYGONS
		      sf::st_coordinates() 
   z <- floor(((co[1,1] + 180)/6) %% 60) + 1   
   if (co[1,2] > 0) return(paste0(z, " N")) else return(paste0(z, " S"))
}
 
#' function to generate an epsg for a given NAD83 UTM zone
#' 
#' This function takes a NAD83 UTM zone as a character string (e.g., "12N")
#'  and returns the epsg.
#'
#' @param z  a NAD83 UTM zone as a character string, e.g., "12N" or "8S".UnitCode #'
#' @return an epsg as an integer number, either 326## for N or 327## for S.
#' 
#' @examples utm2crs("12N")

utm2crs <- function(z) {
   zz <- strsplit(z, " ")[[1]]
   if (zz[2] == "N") return(32600 + as.numeric(zz[1])) else  return(32700 + as.numeric(zz[1]))
}
 
#' This function returns the area in m^2 of a spatial polygon object of 
#' class sf, sfc, or sfg.
#' Simply applying sf::st_area() to objects in any crs can generate invalid areas.
#' This function first reprojects the object into the appropriate UTM zone, 
#' then calculates the area using that projected crs.  This approach shoud be close enough for r
#' reasonable latitude values.
#'
#' @param sfobg  The sf polygon object: sf, sfc, or sfg
#'
#' @return numeric area in units of m^2.
#' 
#' @examples valid_area(FOLSbounds)

valid_area <- function(sfobj) {
   utmcrs <- sfobj |> 
             get_zone() |> 
             utm2crs()
   tmp <- sf::st_transform(sfobj, utmcrs)
   return(sf::st_area(tmp))
}  
 

#' function to return a buffered version of an sf polygon
#' This function takes an sf polygon object in any crs, and returns 
#' another sf object in that same crs, buffered by the specified distance.
#'
#' @param sfobg  The sf polygon object: sfg, sfc, or sf
#' @param buffdist The numeric buffer distance.  Negative values are allowed
#'                 for shrinking polygons.
#' @param bunits The units for the buffer distance, any character string recognized
#'               by the units package, most commonly "m" or "km".
#' @return The buffered polygon as a sfg, sfc, or sf object, the same as the
#'             original object.

buff_poly <- function(sfobj, buffdist = 3, bunits = "km") {
   inCRS <- st_crs(sfobj)
   utmcrs <- sfobj |> 
             get_zone() |> 
             utm2crs()
   tmp <- sfobj |>
          sf::st_transform(crs = utmcrs) |>
          sf::st_buffer(dist = units::set_units(buffdist, bunits, mode = "standard")) |>
          sf::st_transform(crs = inCRS)
   return(tmp)
}


```

## Units To Process

This next chunk specifies which units need to be split into disjunct
subunits, and which need to be spit into adjacent chunks to allow
staggered bounding boxes. These are the units where subdivision is
useful for generating estimates for the entire park. They are not
specifically driven by park managers wanting separate information at the
scale of subunits. If this work gets expanded to provide information at
the subunit scale, there are multiple parks like John Day Fossil Beds
National Monument (JODA) with multiple disjunct subunits at the scale of
tens of km that should be considered.

```{r setGroups}
#| echo: FALSE
#| message: TRUE
#| eval: TRUE 

toSpit <- c("CHIS", "GUIS", "KLGO", "MAPR", "MIIN", "NEPE", "NPSA", "SAGU", "THRO", "TILL")
tooLong <- c("APPA", "BLRI", "NATR", "NOCO")



```

This chunk loads the US Census TIGER state boundary polygons, so
subunits can have their SplitState assigned to them.\
Note that both the NPS boundaries and the Census TIGER state boundaries
claim EPSG = 4269, so the st_transform() in the next chunk should not be
necessary.

```{r loadStates}
states <- tigris::states()


```

## Load Lands Division Boundary File

As of July 2025, the Lands Division "Administrative Boundaries of
National Park System Units - National Geospatial Data Asset (NGDA) NPS
National Park Dataset" is available as a zipped shapefile or public or
NPS AGOL. In previous years the DataStore holding included an ESRI
geodatabase. Note that this script write the new version in multiple
formats, including geopackage.

```{r loadBoundsShp}
#| echo: FALSE
#| eval: TRUE

# allBounds <- sf::st_read("d:/boundaries/Administrative Boundaries of National Park System Units.shp")

allBounds <- sf::st_read(params$boundsshp) |>
             sf::st_make_valid() 
# table(allBounds$UNIT_CODE)

allBounds <- allBounds[order(allBounds$UNIT_CODE),]


if (! st_crs(allBounds) == st_crs(states)) {
  states <- sf::st_transform(states, sf::st_crs(allBounds))
}

```

The next chunk adds the ISO code for Northern Marianas (MP) as STATE for
AMME, following the practice for American Samoa (AS), GUam (GU), Puerto
Rico (PR), and US Virgin Isands (VI). In the June 2025 dataaset, STATE
was missing for AMME.

It then adds the 3 new variables and populates them with the default
values for units requiring no further processing.

```{r addVars}
#| echo: FALSE
#| eval: TRUE

# Add "state" to AMME to follow the AS, GU, PR, & VI pattern for other territories
allBounds$STATE[allBounds$UNIT_CODE == "AMME"] <- "MP"

allBounds$SplitCode <- allBounds$UNIT_CODE
allBounds$GroupCode <- allBounds$UNIT_CODE
allBounds$SplitState <- allBounds$STATE


```

For the already split park and preserve pairs, we want to change the
SplitCode value to be the first 3 letters of UNIT_CODE plus "P" for
preserve. This is how the CRMO + CRMP Craters of the Moon Monument and
Preserve already are in the Lands dataset.

This could be coded a bit cleaner using just the PandP values to create
xcodes, but that would depend on the current pattern of the only
instances of multiple features for the same UNIT_CODE value being
Park/Monument and Preserve pairs. And, there are several stand aone
National Preserves not paired with Parks or Monuments: Big Cypress, Big
Thicket, Mojave, Valles Caldera, etc. This two-step is a bit of
prevention to only do the \*\*\*P SplitCode for Preserves.

```{r dualPreserves}
#| echo: FALSE
#| eval: TRUE

tmp <-  allBounds |> sf::st_drop_geometry() |>
                     dplyr::count(UNIT_CODE)
PandP <- tmp$UNIT_CODE[tmp$n > 1]
PandP


xcodes <- allBounds$SplitCode[allBounds$UNIT_CODE %in% PandP & 
                              allBounds$UNIT_TYPE == "National Preserve"]

xcodes <- paste0(substring(xcodes, 1, 3), "P")

allBounds$SplitCode[allBounds$UNIT_CODE %in% PandP & 
                    allBounds$UNIT_TYPE == "National Preserve"] <- xcodes

rm(tmp, xcodes)

```

The next step is even simpler: creating GroupCode values for units split
into separate features with different verbatim UNIT_CODE values. The
initial set are KICA + SEQU to SEKI, and CRMO + CRMP to CRMO.

```{r addGroup}
#| echo: FALSE
#| eval: TRUE

allBounds$GroupCode[allBounds$UNIT_CODE %in% c("CRMO", "CRMP")] <- "CRMO"
allBounds$GroupCode[allBounds$UNIT_CODE %in% c("KICA", "SEQU")] <- "SEKI"


```

## Separate Out Units for Splitting into Disjunct v Diagonal Subsets

This next chunk divides the park features into the subset that need
splitting disjunct subunits, the long diagonal ones that need clipping
for staggered bounding boxes, and those that need no further processing.

```{r splitUp}
#| echo: FALSE
#| eval: TRUE


# subset to split and process
disjunct <- allBounds[allBounds$UNIT_CODE %in% toSpit,] |>
            sf::st_make_valid() |>
            sf::st_cast(to = "POLYGON")
disjunct$CREATED_BY <- paste(disjunct$CREATED_BY, " TP Split", format(Sys.time(), '%d %B, %Y'))
# Generate centroids of new polygons for semi-reproducible identification.
cent <- disjunct |>
        sf::st_centroid() |>
        sf::st_coordinates()
# add centroid coordinates as X and Y variables (in units of deg lon lat)
disjunct <- cbind(disjunct, cent)

##############################
# Subset to split into segments
diagonal <-  allBounds[allBounds$UNIT_CODE %in% tooLong,] |>
            sf::st_make_valid() |>
            sf::st_cast(to = "POLYGON")
diagonal$CREATED_BY <- "Lands TP Diagonal Split"
# Generate centroids of new polygons for semi-reproducible identification.
cent <- diagonal |>
        sf::st_centroid() |>
        sf::st_coordinates()
# add centroid coordinates as X and Y variables (in units of deg lon lat)
diagonal <- cbind(diagonal, cent)


###############################
# subset to use without further splitting
nochange <- allBounds[! allBounds$UNIT_CODE %in% c(toSpit, tooLong),] |>
            sf::st_make_valid() 
nochange$Area <- nochange |> 
             valid_area() |>
             set_units(value = "ha")

# Generate centroids of new polygons for semi-reproducible identification.
cent <- nochange |>
        sf::st_centroid() |>
        sf::st_coordinates()
# add centroid coordinates as X and Y variables (in units of deg lon lat)
nochange <- cbind(nochange, cent)


```

## Split Widespread Units into Natural Subunits

We don't want to simply split all multipolygon features into individual
features. Too many compact units have multiple polygons immediately
adjacent to each other.

We want to make separate records (features) for specific parks with
widely dispersed units:

-   MAPR: MAPR_Hanford, MAPR_LosAlamos, and MAPR_Oakridge

-   CHIS: CHIS_ADMIN CHIS_Anacapa CHIS_SanMiguel CHIS_SantaBarbara
    CHIS_SantaCruz CHIS_SantaRosa

-   GUIS: GUIS_Florida GUIS_Miss

-   KLGO: KLGO KLSE

-   MIIN: MIIN_Bainbridge MIIN_Idaho

-   NEPE: NEPE_BearPaw_MT NEPE_NezPerce_ID NEPE_NezPerce_OR

-   NPSA: NPSA_Ofu NPSA_Tau NPSA_Tutuila

-   SAGU: SAGU_Rincon SAGU_Tucson

-   THRO: THRO_Elkhorn THRO_North THRO_South

-   TILL: TILL_Chicago TILL_Miss

For the set of units getting their features split, we need a
reproducible way to put the right SplitCode on each separate new
feature. It is dangerous to assume that as updated versions of park
boundaries are produced by Lands, the split geometries will always come
back in the same order, or the centroids will remain constant when park
boundaries change. Therefore, SplitCode values are assigned by
inequlities for centroid X or Y values, which should remain valid except
for very large changes to park boundaries.

Also, note that for KLGO, MAPR, MIIN, and NEPE, we also want to change
some STATE values for the subunits as SplitState.

Identifying appropriate longitude or latitude cuts between subunits can
be done via the [Lands AGOL
application](https://nps.maps.arcgis.com/apps/webappviewer/index.html?id=43bc9db4736140e88e661c67460936e6),
but identifying which polygon in that map corresponds to which polygons
in the shapefile isn't easy. Therefore, this first chunk makes something
of a locator map for each disjunct unit.

```{r splitDispersed}
#| echo: FALSE
#| eval: TRUE



# for ease of use, sort into order and generate a unique identifier
disjunct <- disjunct[order(disjunct$UNIT_CODE, disjunct$X),]
disjunct$index <- as.factor(1:nrow(disjunct))

# plot(tmp["index"])

table(disjunct$STATE)

statemin <- states[,"STUSPS"]
   



tmpstate <-  sf::st_drop_geometry(sf::st_intersection(disjunct, statemin))
tmpstate <- tmpstate[!duplicated(tmpstate$index),c("index", "STUSPS")]
                     
# drop duplicates where 
  
disjunct$SplitState <-  tmpstate$STUSPS[match(disjunct$index, tmpstate$index)]

# CHIS order is Ventura, SantaBarbara, Anacapa, SantaCruz,  SantaRosa, SanMiguel
# NPSA is Tau, Tutuila, Ofu, 
disjunct[,c("UNIT_CODE", "X")]


# sf::st_write(disjunct, "toSpit.kml",  append = FALSE)


# Note that one can't use facet_wrap(..., scales = "free") with geom_sf to make
# a faceted set of local maps.  Therefore, the ggplot needs to loop over units
# and make a separate graph for each unit.

domap <- function(UnitCode) {
   map <- ggplot(data = disjunct[disjunct$UNIT_CODE == UnitCode,], aes(fill = index)) +
          geom_sf() +
          ggtitle(UnitCode)
   return(map)
}  

maps <- llply(toSpit, domap)
lapply(maps, FUN = "print")


# NEPE Special (minus BearPaw)

NEPEmap <- ggplot(data = disjunct[disjunct$UNIT_CODE == "NEPE" & disjunct$X < -110,], aes(fill = index)) +
          geom_sf() +
          ggtitle("NEPE")
NEPEmap

nepe <- disjunct[disjunct$UNIT_CODE == "NEPE",]


# sf::st_write(nepe, "nepe.kml", append = FALSE)


```

The only tricky one is CHIS. From W to E the islands are:
("CHIS_SanMiguel", "CHIS_SantaRosa", "CHIS_SantaCruz", "CHIS_Anacapa",
"CHIS_Ventura", "CHIS_SantaBarbara")

There's also a bit of an issue with GUIS. While the land for GUIS falls
in either MS or FL, the boundary for the easternmost part in Mississippi
(TRACT_ID 01-109) extends approximately 1 mile offshore. While the
eastern end of that offshore is cut sharp, not a curved buffer boundary,
a thin sliver of that offshore / underwater polygon falls into the Tiger
Polygon for Alabama. I don't care which is the official state boundary
offshore, and just glue that sliver back into GUIS_MS,

```{r labelSplits}
#| echo: FALSE
#| eval: TRUE


# CHIS
disjunct$SplitCode[disjunct$UNIT_CODE == "CHIS" & disjunct$X < -120.3] <- "CHIS_SanMiguel"
disjunct$SplitCode[disjunct$UNIT_CODE == "CHIS" & disjunct$X > -120.3 & disjunct$X < -120] <- "CHIS_SantaRosa"
disjunct$SplitCode[disjunct$UNIT_CODE == "CHIS" & disjunct$X > -120 & disjunct$X < -119.5] <- "CHIS_SantaCruz"
disjunct$SplitCode[disjunct$UNIT_CODE == "CHIS" & disjunct$X > -119.5 & disjunct$X < -119.3] <- "CHIS_Anacapa"
disjunct$SplitCode[disjunct$UNIT_CODE == "CHIS" & disjunct$X > -119.3 & disjunct$X < -119.1] <- "CHIS_Ventura"
disjunct$SplitCode[disjunct$UNIT_CODE == "CHIS" & disjunct$X > -119.1] <- "CHIS_SantaBarbara"

# GUIS
disjunct$SplitCode[disjunct$UNIT_CODE == "GUIS" & disjunct$X > -88] <- "GUIS_Florida"  
disjunct$SplitCode[disjunct$UNIT_CODE == "GUIS" & disjunct$X < -88] <- "GUIS_Miss"  
# fix sliver of GUIS_Miss in AL
disjunct$SplitState[disjunct$SplitCode == "GUIS_Miss"] <- "MS"

# KLGO
disjunct$SplitCode[disjunct$UNIT_CODE == "KLGO" & disjunct$X > -130] <- "KLSE"   # KLGO_Seattle


# MAPR
disjunct$SplitCode[disjunct$UNIT_CODE == "MAPR" & disjunct$X < -110] <- "MAPR_Hanford"  
disjunct$SplitCode[disjunct$UNIT_CODE == "MAPR" & disjunct$X > -110 & disjunct$X < -100] <- "MAPR_LosAlamos"  
disjunct$SplitCode[disjunct$UNIT_CODE == "MAPR" & disjunct$X > -100] <- "MAPR_Oakridge"  


# MIIN
disjunct$SplitCode[disjunct$UNIT_CODE == "MIIN" & disjunct$X > -115] <- "MIIN_Idaho"  
disjunct$SplitCode[disjunct$UNIT_CODE == "MIIN" & disjunct$X < -115] <- "MIIN_Bainbridge"  


# NPSA
disjunct$SplitCode[disjunct$UNIT_CODE == "NPSA" & disjunct$X > -169.5] <- "NPSA_Tau"
disjunct$SplitCode[disjunct$UNIT_CODE == "NPSA" & disjunct$X < -170] <- "NPSA_Tutuila"
disjunct$SplitCode[disjunct$UNIT_CODE == "NPSA" & disjunct$X > -170 & disjunct$X < -169.5] <- "NPSA_Ofu"


# SAGU
disjunct$SplitCode[disjunct$UNIT_CODE == "SAGU" & disjunct$X < -111] <- "SAGU_Tucson"  
disjunct$SplitCode[disjunct$UNIT_CODE == "SAGU" & disjunct$X > -111] <- "SAGU_Rincon"  

# THRO
disjunct$SplitCode[disjunct$UNIT_CODE == "THRO" & disjunct$Y > 47.4] <- "THRO_North"  
disjunct$SplitCode[disjunct$UNIT_CODE == "THRO" & disjunct$Y < 47.4 & disjunct$Y > 47] <- "THRO_Elkhorn"  
disjunct$SplitCode[disjunct$UNIT_CODE == "THRO" & disjunct$Y < 47] <- "THRO_South"  

# TILL
disjunct$SplitCode[disjunct$UNIT_CODE == "TILL" & disjunct$Y > 38] <- "TILL_Chicago"  
disjunct$SplitCode[disjunct$UNIT_CODE == "TILL" & disjunct$Y < 38] <- "TILL_Miss"  


# NEPE
disjunct$SplitCode[disjunct$UNIT_CODE == "NEPE" & disjunct$X > -110] <- "NEPE_BearPaw"
disjunct$SplitCode[disjunct$UNIT_CODE == "NEPE" & disjunct$SplitState == "ID"] <- "NEPE_Idaho"
disjunct$SplitCode[disjunct$UNIT_CODE == "NEPE" & disjunct$SplitState == "OR"] <- "NEPE_Oregon"
disjunct$SplitCode[disjunct$UNIT_CODE == "NEPE" & disjunct$SplitState == "WA"] <- "NEPE_Washington"
table(disjunct$SplitCode[disjunct$UNIT_CODE == "NEPE"])
table(disjunct$SplitState[disjunct$UNIT_CODE == "NEPE"])
# 4,5,6 are Spaulding


disjunct[disjunct$UNIT_CODE == "MAPR", c("SplitCode", "SplitState", "X", "Y")]

table(disjunct$SplitCode)

with(disjunct, table(SplitCode, SplitState))
# plot(tmp[tmp$SplitCode == "MAPR_LosAlamos", "index"])

# Now merge individual polygon features within a SplitCode into a
# single MULTIPOLYGON geometry

newgeom <- disjunct |>
  group_by(SplitCode) |> 
  dplyr::summarize()

newdata <- unique(sf::st_drop_geometry(disjunct)[,c(1:21)])

disjunct2 <- merge(newgeom, newdata, by = "SplitCode")


disjunct2$Area <- disjunct2 |> 
             valid_area() |>
             set_units(value = "ha")



tmpstate <- disjunct2 |>
            sf::st_intersection(y = statemin) |>
            sf::st_drop_geometry()
tmpstate <- tmpstate[!duplicated(tmpstate$SplitCode),
                     c("SplitCode", "STUSPS")]
                     
disjunct2$SplitState <- tmpstate$STUSPS[match(disjunct2$SplitCode, 
                                               tmpstate$SplitCode)]


#  junk <- st_drop_geometry(disjunct2[,c("UNIT_CODE", "SplitCode", "SplitState")])
# junk <- junk[order(junk$SplitCode),]

```

## Long Diagonal Units

Next, create segmented parks for the long diagonals: APPA, BLRI, and
NATR. These 2 parks and a trail are relatively mature, with legislated
boundaries.

There are at least 2 more National (Scenic) Trails to keep an eye on.
North Country National Scenic Trail extends from New York to North
Dakota. As of mid-2025, the Lands Division boundaries for it only
include 2 small compact parcels, but that might change. Ice Age National
Scenic Trail (IATR) is anther developing National Trail that currently
has a single small polygon in the Lands Division boundary data. It is
1200 miles long, but folds as an S curve within Wisconsin, so even if
additional parcels become NPS property, the bounding box might not be a
big problem.

While the original code to create stacks of compact bounding boxes left
each of these units as a single multipolygon feature, and split that via
a set of lines along latitudes, a different approach might be better. If
the constituent polygons for each of these diagonal parks are segments
of the park, the centroids of those segments could be used to clump
those segments into blocks with reasonable bounding boxes. The
alternative would be if one polygon was nearly the entire park, with the
other polygons smaller local additions.

The next code chunk breaks each of these multipolygon features into
separate features for each polygon: 193 for APPA, 8 for BLRI, and 10 for
NATR. If these chunks are not 1 or 2 long polygons plus a bunch of
smaller patches, splitting those individual polygons into groups by the
latitude of their centroids prevents any individual polygon from being
sliced in two.

```{r diagonals1}
#| echo: FALSE
#| eval: TRUE

# long linear diagonal units with huge bounding boxes
# these are being dealt with at the level of multi-part bounding boxes

table(diagonal$UNIT_CODE)


diagonal$Area <- diagonal |> 
             valid_area() |>
             set_units(value = "ha")
# bbox2 <- sf::st_bbox(diagonal)
diagonal$index <- paste0(diagonal$UNIT_CODE, "_", 1:nrow(diagonal))

fn0 <- function(sfx) {
  sfx |> sf::st_as_sf() |> 
         sf::st_cast(to = "MULTIPOLYGON") |>
         sf::st_bbox()
} 

bbox2 <- plyr::ddply(diagonal, .(index), fn0)
bbox2$yrange <- bbox2$ymax - bbox2$ymin
bbox2$xrange <- bbox2$xmax - bbox2$xmin
bbox2 <- merge(diagonal[,c("index", "Area", "X", "Y", "SplitState")], bbox2, by = "index")

# are there any really big/long polygons?

bbox2[bbox2$yrange > 1,]

# plot(diagonal[diagonal$UNIT_CODE == "NOCO", "index"])

ggplot(data= diagonal[diagonal$UNIT_CODE == "NOCO",], aes(color = index, fill = index)) +
  geom_sf()


```

So the answer is that APPA and BLRI each have at least 1 single polygon
\> 2.7 degrees latitude range. So plan B of splitting to individual
polygons then lumping into blocks will not reduce the size of bounding
boxes very much. It may work better to go ahead and split units by
latitude, even though that will split some individual polygons.

The coding trick here is that splitting polygons by a set of polylines
isn't a valid spatial operation. So, a set of polygon features are
generated as rectangles using consecutive latitude cutlines and extreme
longitudes, then st_intersection() can be used for the trail boundary
feature and those boxes.

```{r diagonals2}
#| echo: FALSE
#| eval: TRUE

rm(diagonal) # delete from above and start anew


# split APPA into 12 chunks at [34], 35, 36, ..., 44, [46.5]
# split BLRI into 3 chunks at [35], 36, 37, [38.5]
# split NATR into 5 chunks by latitude at [31], 32, 33, 34, 35, [36.5] deg lat
# all can be lon -60 to -170
# split NOCO (North Country National Scenic Trail) into separate chunks based on latitude but mostly longitude.

APPAcuts <- c(34:40, 40.2, 40.4, 40.6, 40.8, 41, 41.6, 43.5, 44, 45, 48) # 16 bboxes
BLRIcuts <- c(34, 35.25, 35.95, 36.4, 37, 39)
NATRcuts <- c(31.8, 33.4, 34.2, 34.7, 37)
NOCOcuts <- c(42, 43, 47)

allcuts <- data.frame(
            UnitCode = c(rep("APPA", length(APPAcuts)),
                         rep("BLRI", length(BLRIcuts)),
                         rep("NATR", length(NATRcuts)),
                         rep("NOCO", length(NOCOcuts))),
            lats = c(APPAcuts, BLRIcuts, NATRcuts, NOCOcuts)             
            )


forbox <- data.frame(
               UnitCode = allcuts$UnitCode,
               xmin = -90,
               ymax= allcuts$lats,
               xmax = -70
)
forbox$xmin[forbox$UnitCode == "NOCO"] <- -95
forbox$xmax[forbox$UnitCode == "NOCO"] <- -80

forbox$ymin <- dplyr::lag(forbox$ymax)

forbox <- forbox[duplicated(forbox$UnitCode),] # discard first for each unitcode

# This uses a syntatically odd form of st_as_sf, where coords as 4 numeric vars
# creates a box
boxes <- sf::st_as_sf(forbox, 
                      coords = c("xmin", "ymin", "xmax", "ymax"),
                      crs = st_crs(allBounds))
boxcount <- dplyr::count(st_drop_geometry(boxes), UnitCode)
boxes$boxID <- paste0(boxes$UnitCode, "_", c(1:boxcount[1,2], 1:boxcount[2,2], 1:boxcount[3, 2]))
  
boxes <- sf::st_make_valid(boxes)

db <- sf::st_make_valid(allBounds[allBounds$UNIT_CODE %in% tooLong,])


splitdiag <- rbind(sf::st_intersection(boxes[boxes$UnitCode == "APPA",],
                                       db[db$UNIT_CODE == "APPA",]),
                    sf::st_intersection(boxes[boxes$UnitCode == "BLRI",],
                                        db[db$UNIT_CODE == "BLRI",]),
                    sf::st_intersection(boxes[boxes$UnitCode == "NATR",],
                                        db[db$UNIT_CODE == "NATR",]),
                    sf::st_intersection(boxes[boxes$UnitCode == "NOCO",],
                                        db[db$UNIT_CODE == "NOCO",])
 
                   )


# nparts <- dplyr::count(st_drop_geometry(splitdiag), UnitCode)

splitdiag$SplitCode <- splitdiag$boxID

splitdiag$Area <- splitdiag |> 
             valid_area() |>
             set_units(value = "ha")
splitdiag$CREATED_BY <- paste(splitdiag$CREATED_BY, " TP Diagonal Split", format(Sys.time(), '%d %B, %Y'))

# statecode <- splitdiag |>
#         sf::st_centroid() |>
#        sf::st_intersection(statemin)



# splitdiag <- sf::st_intersection(splitdiag, statemin)
tmpstate <-  sf::st_drop_geometry(sf::st_intersection(splitdiag, statemin))
# drop duplicates
tmpstate <- tmpstate[!duplicated(tmpstate$boxID),c("boxID", "STUSPS")]
                     

splitdiag$SplitState <- tmpstate$STUSPS[match(splitdiag$boxID, tmpstate$boxID)]

# plot(splitdiag[splitdiag$UNIT_CODE == "NATR", "SplitCode"])


diagMaps <- function(SC) {
   zz <- splitdiag[splitdiag$UNIT_CODE == SC,]
   xx <- sf::st_intersection(boxes[boxes$UnitCode == SC,], st_as_sfc(st_bbox(zz)))
   
   map <- ggplot(data = zz, 
                 aes(color = SplitCode, fill = SplitCode)) +
#          facet_wrap(vars(SplitCode), ncol = 3, scales = "free") +
          geom_sf() +
     geom_sf(data = xx, aes(color = boxID), fill = NA) + 
          ggtitle(SC)
}
dm <- lapply(tooLong, FUN = diagMaps)


for (i in 1:length(dm))  plot(dm[[i]])



# alply(db, .margins = 1, .fun = terra::ext)   
# junk <- alply(splitdiag, .margins = 1, .fun = terra::ext)   

```

The following disabled chunk provided exploration and diagnostics on
these splits.

```{r diagjunk}
#| echo: FALSE
#| eval: FALSE
# st_write(splitdiag, "splitdiag.kml")

BLRImaps <- function(SC) {
           zz <- splitdiag[splitdiag$SplitCode == SC,]
           map <- ggplot(data = zz, 
                 aes(color = SplitCode, fill = SplitCode)) +
#          facet_wrap(vars(SplitCode), ncol = 3, scales = "free") +
          geom_sf() +
          ggtitle(SC)
}

blrim <- lapply(c("BLRI_1", "BLRI_2", "BLRI_3", "BLRI_4"), FUN = BLRImaps )

for (i in 1:length(blrim))  plot(blrim[[i]])

fb <- function(SC) {
  sf::st_bbox(splitdiag[splitdiag$SplitCode == SC,])
}
  
bbb <- lapply(c("BLRI_1", "BLRI_2", "BLRI_3", "BLRI_4"), FUN = fb )
bbb


```

## Recombine

Finally, rejoin the three parts.

```{r rejoin}
#| echo: FALSE
#| message: TRUE
#| eval: TRUE 

names(allBounds)
names(nochange)

names(disjunct2)
names(splitdiag)

nochange$UCaction <- "verbatim"
disjunct2$UCaction <- "split"
splitdiag$UCaction <- "segmented"


namelist <- c("OBJECTID", "UNIT_CODE", "GIS_Notes", "UNIT_NAME",
              "DATE_EDIT", "STATE", "REGION", "GNIS_ID", "UNIT_TYPE", "CREATED_BY",
              "METADATA", "PARKNAME", "CreationDa", "Creator", "EditDate", "Editor",
              "GlobalID", "Shape_Leng", "Shape_Area", 
              "SplitCode", "GroupCode", "SplitState", "Area", "UCaction", "geometry")

TPbounds <- rbind(nochange[, namelist], disjunct2[, namelist], splitdiag[, namelist])
TPbounds <-TPbounds[order(TPbounds$UNIT_CODE, TPbounds$SplitCode),]
```

And write out as geopackage, geojson, and kml.

```{r writeout}
#| echo: FALSE
#| message: TRUE
#| eval: TRUE 


sf::st_write(TPbounds, "BoundaryJuly2025.gpkg", layer = "boundary", append = FALSE)
sf::st_write(allBounds, "BoundaryJuly2025.geojson", layer = "boundary", append = FALSE)
sf::st_write(allBounds, "BoundaryJuly2025.kml", layer = "boundary", append = FALSE)


```

```{r spotCheck}
#| echo: FALSE
#| message: TRUE
#| eval: TRUE 
guis <- ggplot(data = TPbounds[TPbounds$UNIT_CODE == "GUIS",],
               aes(fill = SplitState, color = SplitState)) +
  geom_sf()

guis

guis2 <- ggplot(data = TPbounds[TPbounds$SplitCode == "GUIS_Miss",],
               aes(fill = SplitState, color = SplitState)) +
  geom_sf()

guis2

guis3 <- ggplot(data = TPbounds[TPbounds$UNIT_CODE == "GUIS" &
                                TPbounds$SplitState == "AL",],
               aes(fill = SplitState, color = SplitState)) +
  geom_sf()

guis3

sf::st_drop_geometry(TPbounds[TPbounds$UNIT_CODE == "GUIS",])[,c("SplitCode", "SplitState", "Area")]

```

# Appendices

## Appendix 1: UNIT_CODE,Spit Code, Group Code, UNIT_NAME

```{r ucTable}
kable(sf::st_drop_geometry(TPbounds)[,c("UNIT_CODE", "SplitCode", "GroupCode", "UNIT_NAME")], width = 1200)
```

## Appendix 2: R Code Listing

```{r Listing, ref.label=knitr::all_labels() ,echo=TRUE, eval=FALSE}
```

## Appendix 3: SessionInfo

```{r sessioninfo}

sessionInfo()

```