Getting details of historic orthophoto imagery with R

We would like to obtain historic ortho photo imagery so that we can compare historic watershed conditions compared to current (ex. floodplain vegetation clearing, channel morphology, etc.). First we will generate an area of interest. In our first few code chunks we load our packages and load in some functions that will help us do this work.

Code

suppressMessages(library(tidyverse))
library(ggplot2)
library(bcdata)
library(fwapgr)
library(knitr)
suppressMessages(library(sf))
# library(leaflet)
# library(leafem)

Code

path_post <- fs::path(
  here::here(),
  "posts",
  params$post_name
)

Code

staticimports::import(
  dir = fs::path(
    path_post,
    "scripts"
  ),
  outfile = fs::path(
    path_post,
    "scripts",
    "staticimports",
    ext = "R"
  )
)

Code

source(
  fs::path(
    path_post,
    "scripts",
    "staticimports",
    ext = "R"
  )
)


lfile_name <- function(dat_name = NULL, ext = "geojson") {
  fs::path(
    path_post,
    "data",
    paste0(dat_name, ".", ext)
  )
}

lburn_sf <- function(dat = NULL, dat_name = NULL) {
  if (is.null(dat_name)) {
    cli::cli_abort("You must provide a name for the GeoJSON file using `dat_name`.")
  }
  
  dat |>
    sf::st_write(
      lfile_name(dat_name),
      delete_dsn = TRUE
      # append = FALSE
    )
}

# Function to validate and repair geometries
lngs_geom_validate <- function(layer) {
  layer <- sf::st_make_valid(layer)
  layer <- layer[sf::st_is_valid(layer), ]
  return(layer)
}

Download Spatial Data Layers

Here we download our area of interest which is the Neexdzii Kwah River (a.k.a Upper Bulkley River) which is located between Houston, BC (just south of Smithers) and Topley, BC which is east of Houston and north of Burns Lake, BC. We hit up our remote database managed by Simon Norris with a package built by Poisson Consulting specifically for the task. We use the downstream_route_measure of the Bulkley River (identified through a unique blue_line_key) to query the watershed area upstream of the point where the Neexdzii Kwah River enters the Wedzin Kwah River (a.k.a Morice River).

Code

# lets build a custom watersehed just for upstream of the confluence of Neexdzii Kwa and Wetzin Kwa
# blueline key
blk <- 360873822
# downstream route measure
drm <- 166030.4

aoi <- fwapgr::fwa_watershed_at_measure(blue_line_key = blk, 
                                        downstream_route_measure = drm) |> 
  sf::st_transform(4326)


#get the bounding box of our aoi
# aoi_bb <- sf::st_bbox(aoi)

#lets burn this so we don't need to download each time
aoi <- lngs_geom_validate(aoi)
lburn_sf(
  aoi,
  deparse(substitute(aoi)))

Next we grab a few key layers from the BC Data Catalougue API using convience functions in our rfp package (“Reproducable Field Products”) which wrap the provincially maintained bcdata package. We grab:

Railways
Streams in the Bulkley Watershed group that are 4th order or greater.
Historic Imagery Points
Historic Imagery Polygons
NTS 1:50,000 Grid (we will see why in a second)
Air Photo Centroids

Code

# grab all the railways
l_rail <- rfp::rfp_bcd_get_data(
    bcdata_record_id = "whse_basemapping.gba_railway_tracks_sp"
) |> 
  sf::st_transform(4326) |> 
  janitor::clean_names() 


# streams in the bulkley and then filter to just keep the big ones
l_streams <- rfp::rfp_bcd_get_data(
  bcdata_record_id = "whse_basemapping.fwa_stream_networks_sp",
  col_filter = "watershed_group_code",
  col_filter_value = "BULK",
  # grab a smaller object by including less columns
  col_extract = c("linear_feature_id", "stream_order", "gnis_name", "downstream_route_measure", "blue_line_key", "length_metre")
) |> 
  sf::st_transform(4326) |> 
  janitor::clean_names() |> 
  dplyr::filter(stream_order >= 4)

# historic orthophotos
# WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POLY
#https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points
l_imagery <- rfp::rfp_bcd_get_data(
  bcdata_record_id = "WHSE_IMAGERY_AND_BASE_MAPS.AIMG_ORTHOPHOTO_TILES_POLY") |> 
  sf::st_transform(4326) |> 
  janitor::clean_names()

l_imagery_hist <- rfp::rfp_bcd_get_data(
  bcdata_record_id = "WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT") |> 
  sf::st_transform(4326) |> 
  janitor::clean_names()

l_imagery_grid <- rfp::rfp_bcd_get_data(
  bcdata_record_id = "WHSE_BASEMAPPING.NTS_50K_GRID") |> 
  sf::st_transform(4326) |> 
  janitor::clean_names()

Following download we run some clean up to ensure the geometry of our spatial files is “valid”, trim to our area of interest and burn locally so that every time we rerun iterations of this memo we don’t need to wait for the download process which takes a little longer than we want to wait.

Code

# get a list of the objects in our env that start with l_
ls <- ls()[stringr::str_starts(ls(), "l_")] 

layers_all <- tibble::lst(
  !!!mget(ls)
)

# Apply validation to the AOI and layers
layers_all <- purrr::map(
  layers_all, 
  lngs_geom_validate
  )

# clip them  with purrr and sf
layers_trimmed <- purrr::map(
  layers_all,
  ~ sf::st_intersection(.x, aoi)
) 

# Burn each `sf` object to GeoJSON
purrr::walk2(
  layers_trimmed,
  names(layers_trimmed),
  lburn_sf
)

Code

# lets use the nts mapsheet to query the photo centroids to avoid a massive file download
col_value <- layers_trimmed$l_imagery_grid |> 
  dplyr::pull(map_tile) 

  

l_photo_centroids <- rfp::rfp_bcd_get_data(
  bcdata_record_id = "WHSE_IMAGERY_AND_BASE_MAPS.AIMG_PHOTO_CENTROIDS_SP",
  col_filter = "nts_tile",
  col_filter_value = col_value) |> 
  sf::st_transform(4326) |> 
  janitor::clean_names()

# Apply validation to the AOI and layers
l_photo_centroids <-lngs_geom_validate(l_photo_centroids)

# clip to aoi - can use  layers_trimmed$aoi 
l_photo_centroids <- sf::st_intersection(l_photo_centroids, aoi)


lburn_sf(l_photo_centroids, "l_photo_centroids")

Next - we read the layers back in. The download step is skipped now unless we turn it on again by changing the update_gis param in our memo yaml header to TRUE.

Code

# now we read in all the sf layers that are local so it is really quick
layers_to_load <- fs::dir_ls(
  fs::path(
    path_post,
    "data"),
  glob = "*.geojson"
)

layers_trimmed <- layers_to_load |>
  purrr::map(
    ~ sf::st_read(
      .x, quiet = TRUE)
  ) |> 
  purrr::set_names(
    nm =tools::file_path_sans_ext(
      basename(
        names(
          layers_to_load
        )
      )
    )
  )

OK, seems we cannot get machine readable historical air photo information from the downloaded from the BC data catalogue Historic Imagery Points layer perhaps because the majority of the photos are not georeferenced? What we see in the map and table below (red dot on map) is one point which contains 8 records including links to pdfs and kmls which are basically a georeferenced drawing of where the imagery overlaps (Figure 1). From as far as I can tell - if we wanted to try to use the kmls or pdfs linked in the attribute tables of the “Historic Imagery Points” layer to select orthoimagery we would need to eyeball where the photo polygons overlap where we want to see imagery for and manually write down identifiers for photo by hand. Maybe I am missing something but it sure seems that way.

Code

map <- ggplot() +
  geom_sf(
      data = layers_trimmed$aoi,
      fill = "transparent",
      color = "black",
      linewidth = .5
  ) +
  geom_sf(
    data = layers_trimmed$l_streams,
    color = "blue",
    size = 1
  ) +
  geom_sf(
    data = layers_trimmed$l_rail,
    color = "black",
    size = 1
  ) +
  geom_sf(
    data = layers_trimmed$l_imagery_hist,
    color = "red",
    size = 2
  ) +
  geom_sf(
    data = layers_trimmed$l_imagery_grid,
    alpha = 0.25,
  ) +
  geom_sf_text(
    data = layers_trimmed$l_imagery_grid,
    aes(label = map_tile),
    size = 3  # Adjust size of the text labels as needed
  )

map +
  geom_sf_text(
    data = layers_trimmed$l_streams |> dplyr::distinct(gnis_name, .keep_all = TRUE),
    aes(
      label = gnis_name
    ),
    size = 2  # Adjust size of the text labels as needed
  )

Code

#This what the information in the [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points) layer looks like.

layers_trimmed$l_imagery_hist |> 
  sf::st_drop_geometry() |> 
  knitr::kable()

id	historical_index_map_id	scale_category	geoextent_mapsheet	map_tag	start_year	end_year	kml_url	pdf_url	objectid	se_anno_cad_data	area_ha	localcode_ltree	refine_method	wscode_ltree
WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT.557	557	large	093l_e	093l_e_1	1950	1963	https://openmaps.gov.bc.ca/flight_indices/kml/large_scale/093/93l_e_1_1950_1963.kml	https://openmaps.gov.bc.ca/flight_indices/pdf/large_scale_block_10000_to_25000/093_nts/93l_e/93l_e_sheet_1_20ch_1950_1963.pdf	1860	NA	231944.7	400.431358.585806	CUT	400.431358
WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT.558	558	large	093l_e	093l_e_2	1971	1974	https://openmaps.gov.bc.ca/flight_indices/kml/large_scale/093/93l_e_2_1971_1974.kml	https://openmaps.gov.bc.ca/flight_indices/pdf/large_scale_block_10000_to_25000/093_nts/93l_e/93l_e_sheet_2_20ch_1971_1974.pdf	1861	NA	231944.7	400.431358.585806	CUT	400.431358
WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT.559	559	large	093l_e	093l_e_3	1975	1975	https://openmaps.gov.bc.ca/flight_indices/kml/large_scale/093/93l_e_3_1975.kml	https://openmaps.gov.bc.ca/flight_indices/pdf/large_scale_block_10000_to_25000/093_nts/93l_e/93l_e_sheet_3_20k_1975.pdf	1862	NA	231944.7	400.431358.585806	CUT	400.431358
WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT.560	560	large	093l_e	093l_e_4	1980	1980	https://openmaps.gov.bc.ca/flight_indices/kml/large_scale/093/93l_e_4_1980.kml	https://openmaps.gov.bc.ca/flight_indices/pdf/large_scale_block_10000_to_25000/093_nts/93l_e/93l_e_sheet_4_20k_1980.pdf	1863	NA	231944.7	400.431358.585806	CUT	400.431358
WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT.561	561	large	093l_e	093l_e_5	1980	1980	https://openmaps.gov.bc.ca/flight_indices/kml/large_scale/093/93l_e_5_1980.kml	https://openmaps.gov.bc.ca/flight_indices/pdf/large_scale_block_10000_to_25000/093_nts/93l_e/93l_e_sheet_5_10k_1980.pdf	1864	NA	231944.7	400.431358.585806	CUT	400.431358
WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT.562	562	large	093l_e	093l_e_6	1981	1983	https://openmaps.gov.bc.ca/flight_indices/kml/large_scale/093/93l_e_6_1981_1983.kml	https://openmaps.gov.bc.ca/flight_indices/pdf/large_scale_block_10000_to_25000/093_nts/93l_e/93l_e_sheet_6_20k_1981_1983.pdf	1865	NA	231944.7	400.431358.585806	CUT	400.431358
WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT.563	563	large	093l_e	093l_e_7	1989	1989	https://openmaps.gov.bc.ca/flight_indices/kml/large_scale/093/93l_e_7_1989.kml	https://openmaps.gov.bc.ca/flight_indices/pdf/large_scale_block_10000_to_25000/093_nts/93l_e/93l_e_sheet_7_15k_1989.pdf	1866	NA	231944.7	400.431358.585806	CUT	400.431358
WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT.564	564	large	093l_e	093l_e_8	1990	1990	https://openmaps.gov.bc.ca/flight_indices/kml/large_scale/093/93l_e_8_1990.kml	https://openmaps.gov.bc.ca/flight_indices/pdf/large_scale_block_10000_to_25000/093_nts/93l_e/93l_e_sheet_8_10k_15k_1990_bw_1990_colour.pdf	1867	NA	231944.7	400.431358.585806	CUT	400.431358

Code

my_caption <- "Screenshot of kml downloaded from link provided in Historic Imagery Points."
knitr::include_graphics(fs::path(
  path_post,
  "fig",
  "Screenshot1",
  ext = "png"
  )
)

Figure 1: Screenshot of kml downloaded from link provided in Historic Imagery Points.

For the Historic Imagery Polygons layer the range of year_operational is 1996, 2013. This is not as far back as we would prefer to be looking.

It does however seem that each of the Air Photo Centroids are georeferenced with a date range of:

Code

range(layers_trimmed$l_photo_centroids$photo_date)

[1] "1963-08-07" "2019-09-18"

Code

# At this point we have downloaded two csvs (one for each NTS 1:50,000 mapsheet of course) with information about the airphotos including UTM coordinates that we will assume for now are the photo centres. In our next steps we read in what we have, turn into spatial object, trim to overall study area and plot.
# list csvs
ls <- fs::dir_ls(
  fs::path(
    path_post,
    "data"),
  glob = "*.csv"
)

photos_raw <- ls |> 
  purrr::map_df(
    readr::read_csv
  ) |> 
  sf::st_as_sf(
    coords = c("Longitude", "Latitude"), crs = 4326
  ) |> 
  janitor::clean_names() |> 
  dplyr::mutate(photo_date = lubridate::mdy(photo_date)) 


photos_aoi <- sf::st_intersection(
  photos_raw, 
  layers_trimmed$aoi |> st_make_valid()
  )

Code

map +
  geom_sf(
    data = layers_trimmed$l_photo_centroids,
    alpha = 0.25
  )

That is a lot of photos! 7910 photos to be exact!!!

Code

# amount to buffer all stream segments
q_buffer <- 500
# q_drm_main <- 263795

# length of streams other than selected explicity to buffer
q_drm_other <- 3000

Clip Photo Information with Streams Buffers

Here are our query parameters to narrow down the area within our study are watershed in which we want to find photos for:

Buffer: 500m - size of buffer used on either side of stream lines selected
Stream segments:
- Bulkley River (gnis_name in the stream layer)
- Maxan Creek
- Buck Creek
- for each remaining stream - segments of that stream which begin before 3000m from the downstream system (i.e. the first 3km) of stream.

Code

# We use the `downstream_route_measure` of the stream layer to exclude areas upstream of Bulkley Lake (also known as Taman Creek).  We find it in QGIS by highlighting the stream layer and clicking on our segment of interest while we have the information tool selected - the resulting pop-up looks like this in QGIS.
knitr::include_graphics(fs::path(
  path_post,
  "fig",
  "Screenshot2",
  ext = "png"
  )
)

Code

r_streams <- c("Maxan Creek", "Buck Creek")

aoi_refined_raw <- layers_trimmed$l_streams |> 
  # removed  & downstream_route_measure < q_drm_main for bulkley as doestn't cahnge 1960s query and increases beyond just by 5 photos
  dplyr::filter(gnis_name == "Bulkley River"|
                  gnis_name != "Bulkley River" & downstream_route_measure < q_drm_other |
                  gnis_name %in% r_streams) |> 
  # dplyr::arrange(downstream_route_measure) |>
  # calculate when we get to length_m by adding up the length_metre field and filtering out everything up to length_m
  # dplyr::filter(cumsum(length_metre) <= length_m) |>
  sf::st_union() |> 
  # we need to run st_sf or we get a sp object in a list...
  sf::st_sf()
  
aoi_refined_buffered <- sf::st_buffer(
  aoi_refined_raw,
  q_buffer, endCapStyle = "FLAT"
) 

photos_aoi_refined <- sf::st_intersection(
  layers_trimmed$l_photo_centroids, 
  aoi_refined_buffered
  )

Let’s plot again and include our buffered areas around the first 3000m of streams (area in red) along with the location of the photo points that land within that area. Looks like this give us 1361 photos.

Code

map +
  geom_sf(
    data = aoi_refined_buffered,
    color = "red",
    alpha= 0
  ) +
  geom_sf(
    data = photos_aoi_refined,
    alpha = 0.25,
  ) +
  geom_sf_text(
    data = layers_trimmed$l_streams |> dplyr::distinct(gnis_name, .keep_all = TRUE),
    aes(
      label = gnis_name
    ),
    size = 2  # Adjust size of the text labels as needed
  )

That is not as many photos - but still quite a few (1361). Figure 2 below can be used to filter these photos from any time and/or mapsheet and export the result to csv or excel file.

Code

photos_aoi_refined |> 
  dplyr::select(-id) |> 
  my_dt_table(cols_freeze_left = 0)

Figure 2

Filter Photos by Date

Now lets map by year to see what our options are including the earliest photos possible. Here is our range to choose from:

Code

range(photos_aoi_refined$photo_date)

[1] "1963-08-07" "2019-09-18"

Code

map +
geom_sf(
  data = photos_aoi_refined |> dplyr::filter(photo_year <= "1975")
  ) +
  facet_wrap(~ photo_year)

Well - looks like we get really good coverage of the Bulkley River mainstem in 1968 then much better coverage of the Buck Creek drainage and Maxan Creek in 1971. For 1975 - the coverage of the Bulkley mainstem and Maxan Creek is pretty good…

Thinking the ideal thing to do is to grab the photos from:

1968 all
1971 for the Buck Creek and Maxan Creek areas only
1975 Maxan Creek only

Code

# spatially represent just Buck and Maxan, buffer and clip the 1971 photos
# "r_" is for "refine"
r_year1 <- "1968"
r_year2 <- "1971"
r_year3 <- "1975"

r_streams2 <- c("Maxan Creek")

l_streams_refined1 <- layers_trimmed$l_streams |> 
  # we defined r_streams in chunk way above 
  dplyr::filter(gnis_name %in% r_streams) |> 
  sf::st_union() |> 
  # we need to run st_sf or we get a sp object in a list...
  sf::st_sf()
  
aoi_refined_buffered2 <- sf::st_buffer(
  l_streams_refined1,
  q_buffer, endCapStyle = "FLAT"
) 

l_streams_refined2 <- layers_trimmed$l_streams |> 
  # we defined r_streams in chunk way above 
  dplyr::filter(gnis_name %in% r_streams2) |> 
  sf::st_union() |> 
  # we need to run st_sf or we get a sp object in a list...
  sf::st_sf()
  
aoi_refined_buffered3 <- sf::st_buffer(
  l_streams_refined2,
  q_buffer, endCapStyle = "FLAT"
) 

# filter first year
photos1 <- photos_aoi_refined |> 
  dplyr::filter(
      photo_year == r_year1
  )

# filter second year using just the streams we want to include
photos2 <- sf::st_intersection(
  layers_trimmed$l_photo_centroids |> dplyr::filter(photo_year == r_year2), 
  aoi_refined_buffered2
  )

# filter second year using just the streams we want to include
photos3 <- sf::st_intersection(
  layers_trimmed$l_photo_centroids |> dplyr::filter(photo_year == r_year3), 
  aoi_refined_buffered3
  )

photos_all <- dplyr::bind_rows(photos1, photos2, photos3)

Now let’s have a look at the individual year components (Figure 3) as well as the whole dataset (Figure 4). We are privileged to potentially have the assistance of Mike Price to help us obtain this imagery from the UBC archives. If there are too many photos to grab as is - the table below can be filtered by photo_year to reduce the number of photos. The resulting filtered dataset can then be downloaded by pressing the CSV or Excel buttons at the bottom of the table….

Code

my_caption <- "Amalgamated photo points presented by year."
map +
  geom_sf(
  data = photos_all
  ) +
  facet_wrap(~ photo_year)

Figure 3: Amalgamated photo points presented by year.

Code

my_caption <- "Amalgamated photo points"
map +
  geom_sf(
  data = photos_all
  ) +
  geom_sf_text(
    data = layers_trimmed$l_streams |> dplyr::distinct(gnis_name, .keep_all = TRUE),
    aes(
      label = gnis_name
    ),
    size = 2  # Adjust size of the text labels as needed
  )

Code

photos_all |> 
  dplyr::select(-id) |> 
  my_dt_table(cols_freeze_left = 0)

Export `csv` with Photo Information

Let’s burn out csv that can be used to find the imagery for the 136 photos above.

Code

lfile_name_photos <- function(dat = NULL){
  fs::path(
      path_post,
      "exports",
      paste(
        "airphotos",
        paste(range(dat$photo_date), collapse = "_"),
        sep = "_"
      ),
      ext = "csv"
    )
}

photos_all |> 
  readr::write_csv(
    lfile_name_photos(photos_all), na =""
  )


lpath_link <- function(dat = NULL){
  paste0(
    "https://github.com/NewGraphEnvironment/new_graphiti/tree/main/posts/2024-11-15-bcdata-ortho-historic/exports/",
    basename(
      lfile_name_photos(dat)
    )
  )
}

We can view and download exported csv files here.

Download Spatial Data Layers

Clip Photo Information with Streams Buffers

Filter Photos by Date

Export csv with Photo Information

Export `csv` with Photo Information