diff --git a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/execute-results/html.json b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/execute-results/html.json
index b5f7a85..2489851 100644
--- a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/execute-results/html.json
+++ b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/execute-results/html.json
@@ -1,8 +1,8 @@
{
- "hash": "a1b58ca05772f760e1a2a30d7ff2c945",
+ "hash": "42885393bce698c9ec287fa2ea31ce0f",
"result": {
"engine": "knitr",
- "markdown": "---\ntitle: \"Getting details of historic orthophoto imagery with R\"\nauthor: \"al\"\ndate: \"2024-11-15\"\ndate-modified: \"2024-11-17\"\ncategories: [fwapg, r, bcdata, imagery, api]\nimage: \"image.jpg\"\nparams:\n repo_owner: \"NewGraphEnvironment\"\n repo_name: \"new_graphiti\"\n post_name: \"2024-11-15-bcdata-ortho-historic\"\n update_gis: FALSE\nformat: \n html:\n code-fold: true\n---\n\n\nWe 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.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsuppressMessages(library(tidyverse))\nlibrary(ggplot2)\nlibrary(bcdata)\nlibrary(fwapgr)\nlibrary(knitr)\nsuppressMessages(library(sf))\n# library(leaflet)\n# library(leafem)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\npath_post <- fs::path(\n here::here(),\n \"posts\",\n params$post_name\n)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nstaticimports::import(\n dir = fs::path(\n path_post,\n \"scripts\"\n ),\n outfile = fs::path(\n path_post,\n \"scripts\",\n \"staticimports\",\n ext = \"R\"\n )\n)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nsource(\n fs::path(\n path_post,\n \"scripts\",\n \"staticimports\",\n ext = \"R\"\n )\n)\n\n\nlfile_name <- function(dat_name = NULL, ext = \"geojson\") {\n fs::path(\n path_post,\n \"data\",\n paste0(dat_name, \".\", ext)\n )\n}\n\nlburn_sf <- function(dat = NULL, dat_name = NULL) {\n if (is.null(dat_name)) {\n cli::cli_abort(\"You must provide a name for the GeoJSON file using `dat_name`.\")\n }\n \n dat |>\n sf::st_write(\n lfile_name(dat_name),\n delete_dsn = TRUE\n # append = FALSE\n )\n}\n\n# Function to validate and repair geometries\nlngs_geom_validate <- function(layer) {\n layer <- sf::st_make_valid(layer)\n layer <- layer[sf::st_is_valid(layer), ]\n return(layer)\n}\n```\n:::\n\n\n\n## Download Spatial Data Layers\nHere 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).\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# lets build a custom watersehed just for upstream of the confluence of Neexdzii Kwa and Wetzin Kwa\n# blueline key\nblk <- 360873822\n# downstream route measure\ndrm <- 166030.4\n\naoi <- fwapgr::fwa_watershed_at_measure(blue_line_key = blk, \n downstream_route_measure = drm) |> \n sf::st_transform(4326)\n\n\n#get the bounding box of our aoi\n# aoi_bb <- sf::st_bbox(aoi)\n\n#lets burn this so we don't need to download each time\naoi <- lngs_geom_validate(aoi)\nlburn_sf(\n aoi,\n deparse(substitute(aoi)))\n```\n:::\n\n\nNext 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:\n\n - Railways\n - Streams in the Bulkley Watershed group that are 4th order or greater.\n - [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points)\n - [Historic Imagery Polygons](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-polygons)\n - [NTS 1:50,000 Grid](https://catalogue.data.gov.bc.ca/) (we will see why in a second)\n - [Air Photo Centroids](https://catalogue.data.gov.bc.ca/dataset/airphoto-centroids)\n \n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# grab all the railways\nl_rail <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"whse_basemapping.gba_railway_tracks_sp\"\n) |> \n sf::st_transform(4326) |> \n janitor::clean_names() \n\n\n# streams in the bulkley and then filter to just keep the big ones\nl_streams <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"whse_basemapping.fwa_stream_networks_sp\",\n col_filter = \"watershed_group_code\",\n col_filter_value = \"BULK\",\n # grab a smaller object by including less columns\n col_extract = c(\"linear_feature_id\", \"stream_order\", \"gnis_name\", \"downstream_route_measure\", \"blue_line_key\", \"length_metre\")\n) |> \n sf::st_transform(4326) |> \n janitor::clean_names() |> \n dplyr::filter(stream_order >= 4)\n\n# historic orthophotos\n# WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POLY\n#https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points\nl_imagery <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"WHSE_IMAGERY_AND_BASE_MAPS.AIMG_ORTHOPHOTO_TILES_POLY\") |> \n sf::st_transform(4326) |> \n janitor::clean_names()\n\nl_imagery_hist <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT\") |> \n sf::st_transform(4326) |> \n janitor::clean_names()\n\nl_imagery_grid <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"WHSE_BASEMAPPING.NTS_50K_GRID\") |> \n sf::st_transform(4326) |> \n janitor::clean_names()\n```\n:::\n\n\nFollowing 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.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# get a list of the objects in our env that start with l_\nls <- ls()[stringr::str_starts(ls(), \"l_\")] \n\nlayers_all <- tibble::lst(\n !!!mget(ls)\n)\n\n# Apply validation to the AOI and layers\nlayers_all <- purrr::map(\n layers_all, \n lngs_geom_validate\n )\n\n# clip them with purrr and sf\nlayers_trimmed <- purrr::map(\n layers_all,\n ~ sf::st_intersection(.x, aoi)\n) \n\n# Burn each `sf` object to GeoJSON\npurrr::walk2(\n layers_trimmed,\n names(layers_trimmed),\n lburn_sf\n)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n# lets use the nts mapsheet to query the photo centroids to avoid a massive file download\ncol_value <- layers_trimmed$l_imagery_grid |> \n dplyr::pull(map_tile) \n\n \n\nl_photo_centroids <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"WHSE_IMAGERY_AND_BASE_MAPS.AIMG_PHOTO_CENTROIDS_SP\",\n col_filter = \"nts_tile\",\n col_filter_value = col_value) |> \n sf::st_transform(4326) |> \n janitor::clean_names()\n\n# Apply validation to the AOI and layers\nl_photo_centroids <-lngs_geom_validate(l_photo_centroids)\n\n# clip to aoi - can use layers_trimmed$aoi \nl_photo_centroids <- sf::st_intersection(l_photo_centroids, aoi)\n\n\nlburn_sf(l_photo_centroids, \"l_photo_centroids\")\n```\n:::\n\n\nNext - 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`.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# now we read in all the sf layers that are local so it is really quick\nlayers_to_load <- fs::dir_ls(\n fs::path(\n path_post,\n \"data\"),\n glob = \"*.geojson\"\n)\n\nlayers_trimmed <- layers_to_load |>\n purrr::map(\n ~ sf::st_read(\n .x, quiet = TRUE)\n ) |> \n purrr::set_names(\n nm =tools::file_path_sans_ext(\n basename(\n names(\n layers_to_load\n )\n )\n )\n )\n```\n:::\n\n\nOK, seems we cannot get machine readable historical air photo information from the downloaded from the BC data catalogue [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-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 (@fig-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. \n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap <- ggplot() +\n geom_sf(\n data = layers_trimmed$aoi,\n fill = \"transparent\",\n color = \"black\",\n linewidth = .5\n ) +\n geom_sf(\n data = layers_trimmed$l_streams,\n color = \"blue\",\n size = 1\n ) +\n geom_sf_text(\n data = layers_trimmed$l_streams |> dplyr::distinct(gnis_name, .keep_all = TRUE),\n aes(\n label = gnis_name\n ),\n size = 2 # Adjust size of the text labels as needed\n ) +\n geom_sf(\n data = layers_trimmed$l_rail,\n color = \"black\",\n size = 1\n ) +\n geom_sf(\n data = layers_trimmed$l_imagery_hist,\n color = \"red\",\n size = 2\n ) +\n geom_sf(\n data = layers_trimmed$l_imagery_grid,\n alpha = 0.25,\n ) +\n geom_sf_text(\n data = layers_trimmed$l_imagery_grid,\n aes(label = map_tile),\n size = 3 # Adjust size of the text labels as needed\n )\n\nmap\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n#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.\n\nlayers_trimmed$l_imagery_hist |> \n sf::st_drop_geometry() |> \n knitr::kable()\n```\n\n::: {.cell-output-display}\n\n\n|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 |\n|:---------------------------------------------------------|-----------------------:|:--------------|:------------------|:--------|----------:|--------:|:-----------------------------------------------------------------------------------|:------------------------------------------------------------------------------------------------------------------------------------------|--------:|:----------------|--------:|:-----------------|:-------------|:------------|\n|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 |\n|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 |\n|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 |\n|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 |\n|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 |\n|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 |\n|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 |\n|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 |\n\n\n:::\n:::\n\n \n
\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmy_caption <- \"Screenshot of kml downloaded from link provided in Historic Imagery Points.\"\nknitr::include_graphics(fs::path(\n path_post,\n \"fig\",\n \"Screenshot1\",\n ext = \"png\"\n )\n)\n```\n\n::: {.cell-output-display}\n{#fig-1 width=1492}\n:::\n:::\n\n\n\n\n
\n \nFor the [Historic Imagery Polygons](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-polygons) layer the range of `year_operational` is 1996, 2013. This is not as far back as we would prefer to be looking.\n \n\n
\n\nIt does however seem that each of the [Air Photo Centroids](https://catalogue.data.gov.bc.ca/dataset/airphoto-centroids) are georeferenced with a date range of:\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrange(layers_trimmed$l_photo_centroids$photo_date)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] \"1963-08-07\" \"2019-09-18\"\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n# 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.\n# list csvs\nls <- fs::dir_ls(\n fs::path(\n path_post,\n \"data\"),\n glob = \"*.csv\"\n)\n\nphotos_raw <- ls |> \n purrr::map_df(\n readr::read_csv\n ) |> \n sf::st_as_sf(\n coords = c(\"Longitude\", \"Latitude\"), crs = 4326\n ) |> \n janitor::clean_names() |> \n dplyr::mutate(photo_date = lubridate::mdy(photo_date)) \n\n\nphotos_aoi <- sf::st_intersection(\n photos_raw, \n layers_trimmed$aoi |> st_make_valid()\n )\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = layers_trimmed$l_photo_centroids,\n alpha = 0.25\n ) \n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\nThat is a lot of photos! 7910 photos to be exact!!!\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# amount to buffer all stream segments\nq_buffer <- 500\n# q_drm_main <- 263795\n\n# length of streams other than selected explicity to buffer\nq_drm_other <- 3000\n```\n:::\n\n\n## Clip Photo Information with Streams Buffers\nHere are our query parameters to narrow down the area within our study are watershed in which we want to find photos for:\n\n - Buffer: 500m - size of buffer used on either side of stream lines selected\n - Stream segments: \n + Bulkley River (`gnis_name` in the stream layer)\n + Maxan Creek\n + Buck Creek\n + for each remaining stream - segments of that stream which begin before 3000m from the downstream system (i.e. the first 3km) of stream.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# 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.\nknitr::include_graphics(fs::path(\n path_post,\n \"fig\",\n \"Screenshot2\",\n ext = \"png\"\n )\n)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nr_streams <- c(\"Maxan Creek\", \"Buck Creek\")\n\naoi_refined_raw <- layers_trimmed$l_streams |> \n # removed & downstream_route_measure < q_drm_main for bulkley as doestn't cahnge 1960s query and increases beyond just by 5 photos\n dplyr::filter(gnis_name == \"Bulkley River\"|\n gnis_name != \"Bulkley River\" & downstream_route_measure < q_drm_other |\n gnis_name %in% r_streams) |> \n # dplyr::arrange(downstream_route_measure) |>\n # calculate when we get to length_m by adding up the length_metre field and filtering out everything up to length_m\n # dplyr::filter(cumsum(length_metre) <= length_m) |>\n sf::st_union() |> \n # we need to run st_sf or we get a sp object in a list...\n sf::st_sf()\n \naoi_refined_buffered <- sf::st_buffer(\n aoi_refined_raw,\n q_buffer, endCapStyle = \"FLAT\"\n) \n\nphotos_aoi_refined <- sf::st_intersection(\n layers_trimmed$l_photo_centroids, \n aoi_refined_buffered\n )\n```\n:::\n\n\nLet'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.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = aoi_refined_buffered,\n color = \"red\",\n alpha= 0\n ) +\n geom_sf(\n data = photos_aoi_refined,\n alpha = 0.25,\n ) \n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\nThat is not as many photos - but still quite a few (1361). The table below can be used to filter these photos from any time and/or mapsheet and export the result to csv or excel file. \n\n\n::: {.cell}\n\n```{.r .cell-code}\nphotos_aoi_refined |> \n my_dt_table(cols_freeze_left = 2)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n\n\n```\n\n:::\n:::\n\n\n## Filter Photos by Date\nNow lets map by year to see what our options are including the earliest photos possible. Here is our range to choose from:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrange(photos_aoi_refined$photo_date)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] \"1963-08-07\" \"2019-09-18\"\n```\n\n\n:::\n:::\n\n`\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\ngeom_sf(\n data = photos_aoi_refined |> dplyr::filter(photo_year <= \"1975\")\n ) +\n facet_wrap(~ photo_year)\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\nWell - 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...\n\n
\n\nThinking the ideal thing to do is to grab the photos from:\n\n - 1968 all\n - 1971 for the Buck Creek and Maxan Creek areas only\n - 1975 Maxan Creek only\n\n
\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# spatially represent just Buck and Maxan, buffer and clip the 1971 photos\n# \"r_\" is for \"refine\"\nr_year1 <- \"1968\"\nr_year2 <- \"1971\"\nr_year3 <- \"1975\"\n\nr_streams2 <- c(\"Maxan Creek\")\n\nl_streams_refined1 <- layers_trimmed$l_streams |> \n # we defined r_streams in chunk way above \n dplyr::filter(gnis_name %in% r_streams) |> \n sf::st_union() |> \n # we need to run st_sf or we get a sp object in a list...\n sf::st_sf()\n \naoi_refined_buffered2 <- sf::st_buffer(\n l_streams_refined1,\n q_buffer, endCapStyle = \"FLAT\"\n) \n\nl_streams_refined2 <- layers_trimmed$l_streams |> \n # we defined r_streams in chunk way above \n dplyr::filter(gnis_name %in% r_streams2) |> \n sf::st_union() |> \n # we need to run st_sf or we get a sp object in a list...\n sf::st_sf()\n \naoi_refined_buffered3 <- sf::st_buffer(\n l_streams_refined2,\n q_buffer, endCapStyle = \"FLAT\"\n) \n\n# filter first year\nphotos1 <- photos_aoi_refined |> \n dplyr::filter(\n photo_year == r_year1\n )\n\n# filter second year using just the streams we want to include\nphotos2 <- sf::st_intersection(\n layers_trimmed$l_photo_centroids |> dplyr::filter(photo_year == r_year2), \n aoi_refined_buffered2\n )\n\n# filter second year using just the streams we want to include\nphotos3 <- sf::st_intersection(\n layers_trimmed$l_photo_centroids |> dplyr::filter(photo_year == r_year3), \n aoi_refined_buffered3\n )\n\nphotos_all <- dplyr::bind_rows(photos1, photos2, photos3)\n```\n:::\n\n\n\nNow let's have a look at what we have\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = photos_all\n )\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n## Export `csv` with Photo Information\nLet's burn out csv that can be used to find the imagery for the 136 photos above.\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlfile_name_photos <- function(dat = NULL){\n fs::path(\n path_post,\n \"exports\",\n paste(\n \"airphotos\",\n paste(range(dat$photo_date), collapse = \"_\"),\n sep = \"_\"\n ),\n ext = \"csv\"\n )\n}\n\nphotos_all |> \n readr::write_csv(\n lfile_name_photos(photos_all), na =\"\"\n )\n\n\nlpath_link <- function(dat = NULL){\n paste0(\n \"https://github.com/NewGraphEnvironment/new_graphiti/tree/main/posts/2024-11-15-bcdata-ortho-historic/exports/\",\n basename(\n lfile_name_photos(dat)\n )\n )\n}\n```\n:::\n\n\nWe can view and download exported csv files [here](https://github.com/NewGraphEnvironment/new_graphiti/tree/main/posts/2024-11-15-bcdata-ortho-historic/exports/). \n\n\n\n\n",
+ "markdown": "---\ntitle: \"Getting details of historic orthophoto imagery with R\"\nauthor: \"al\"\ndate: \"2024-11-15\"\ndate-modified: \"2024-11-17\"\ncategories: [fwapg, r, bcdata, imagery, api]\nimage: \"image.jpg\"\nparams:\n repo_owner: \"NewGraphEnvironment\"\n repo_name: \"new_graphiti\"\n post_name: \"2024-11-15-bcdata-ortho-historic\"\n update_gis: FALSE\nformat: \n html:\n code-fold: true\n---\n\n\nWe 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.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsuppressMessages(library(tidyverse))\nlibrary(ggplot2)\nlibrary(bcdata)\nlibrary(fwapgr)\nlibrary(knitr)\nsuppressMessages(library(sf))\n# library(leaflet)\n# library(leafem)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\npath_post <- fs::path(\n here::here(),\n \"posts\",\n params$post_name\n)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nstaticimports::import(\n dir = fs::path(\n path_post,\n \"scripts\"\n ),\n outfile = fs::path(\n path_post,\n \"scripts\",\n \"staticimports\",\n ext = \"R\"\n )\n)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nsource(\n fs::path(\n path_post,\n \"scripts\",\n \"staticimports\",\n ext = \"R\"\n )\n)\n\n\nlfile_name <- function(dat_name = NULL, ext = \"geojson\") {\n fs::path(\n path_post,\n \"data\",\n paste0(dat_name, \".\", ext)\n )\n}\n\nlburn_sf <- function(dat = NULL, dat_name = NULL) {\n if (is.null(dat_name)) {\n cli::cli_abort(\"You must provide a name for the GeoJSON file using `dat_name`.\")\n }\n \n dat |>\n sf::st_write(\n lfile_name(dat_name),\n delete_dsn = TRUE\n # append = FALSE\n )\n}\n\n# Function to validate and repair geometries\nlngs_geom_validate <- function(layer) {\n layer <- sf::st_make_valid(layer)\n layer <- layer[sf::st_is_valid(layer), ]\n return(layer)\n}\n```\n:::\n\n\n\n## Download Spatial Data Layers\nHere 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).\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# lets build a custom watersehed just for upstream of the confluence of Neexdzii Kwa and Wetzin Kwa\n# blueline key\nblk <- 360873822\n# downstream route measure\ndrm <- 166030.4\n\naoi <- fwapgr::fwa_watershed_at_measure(blue_line_key = blk, \n downstream_route_measure = drm) |> \n sf::st_transform(4326)\n\n\n#get the bounding box of our aoi\n# aoi_bb <- sf::st_bbox(aoi)\n\n#lets burn this so we don't need to download each time\naoi <- lngs_geom_validate(aoi)\nlburn_sf(\n aoi,\n deparse(substitute(aoi)))\n```\n:::\n\n\nNext 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:\n\n - Railways\n - Streams in the Bulkley Watershed group that are 4th order or greater.\n - [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points)\n - [Historic Imagery Polygons](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-polygons)\n - [NTS 1:50,000 Grid](https://catalogue.data.gov.bc.ca/) (we will see why in a second)\n - [Air Photo Centroids](https://catalogue.data.gov.bc.ca/dataset/airphoto-centroids)\n \n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# grab all the railways\nl_rail <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"whse_basemapping.gba_railway_tracks_sp\"\n) |> \n sf::st_transform(4326) |> \n janitor::clean_names() \n\n\n# streams in the bulkley and then filter to just keep the big ones\nl_streams <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"whse_basemapping.fwa_stream_networks_sp\",\n col_filter = \"watershed_group_code\",\n col_filter_value = \"BULK\",\n # grab a smaller object by including less columns\n col_extract = c(\"linear_feature_id\", \"stream_order\", \"gnis_name\", \"downstream_route_measure\", \"blue_line_key\", \"length_metre\")\n) |> \n sf::st_transform(4326) |> \n janitor::clean_names() |> \n dplyr::filter(stream_order >= 4)\n\n# historic orthophotos\n# WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POLY\n#https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points\nl_imagery <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"WHSE_IMAGERY_AND_BASE_MAPS.AIMG_ORTHOPHOTO_TILES_POLY\") |> \n sf::st_transform(4326) |> \n janitor::clean_names()\n\nl_imagery_hist <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"WHSE_IMAGERY_AND_BASE_MAPS.AIMG_HIST_INDEX_MAPS_POINT\") |> \n sf::st_transform(4326) |> \n janitor::clean_names()\n\nl_imagery_grid <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"WHSE_BASEMAPPING.NTS_50K_GRID\") |> \n sf::st_transform(4326) |> \n janitor::clean_names()\n```\n:::\n\n\nFollowing 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.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# get a list of the objects in our env that start with l_\nls <- ls()[stringr::str_starts(ls(), \"l_\")] \n\nlayers_all <- tibble::lst(\n !!!mget(ls)\n)\n\n# Apply validation to the AOI and layers\nlayers_all <- purrr::map(\n layers_all, \n lngs_geom_validate\n )\n\n# clip them with purrr and sf\nlayers_trimmed <- purrr::map(\n layers_all,\n ~ sf::st_intersection(.x, aoi)\n) \n\n# Burn each `sf` object to GeoJSON\npurrr::walk2(\n layers_trimmed,\n names(layers_trimmed),\n lburn_sf\n)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n# lets use the nts mapsheet to query the photo centroids to avoid a massive file download\ncol_value <- layers_trimmed$l_imagery_grid |> \n dplyr::pull(map_tile) \n\n \n\nl_photo_centroids <- rfp::rfp_bcd_get_data(\n bcdata_record_id = \"WHSE_IMAGERY_AND_BASE_MAPS.AIMG_PHOTO_CENTROIDS_SP\",\n col_filter = \"nts_tile\",\n col_filter_value = col_value) |> \n sf::st_transform(4326) |> \n janitor::clean_names()\n\n# Apply validation to the AOI and layers\nl_photo_centroids <-lngs_geom_validate(l_photo_centroids)\n\n# clip to aoi - can use layers_trimmed$aoi \nl_photo_centroids <- sf::st_intersection(l_photo_centroids, aoi)\n\n\nlburn_sf(l_photo_centroids, \"l_photo_centroids\")\n```\n:::\n\n\nNext - 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`.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# now we read in all the sf layers that are local so it is really quick\nlayers_to_load <- fs::dir_ls(\n fs::path(\n path_post,\n \"data\"),\n glob = \"*.geojson\"\n)\n\nlayers_trimmed <- layers_to_load |>\n purrr::map(\n ~ sf::st_read(\n .x, quiet = TRUE)\n ) |> \n purrr::set_names(\n nm =tools::file_path_sans_ext(\n basename(\n names(\n layers_to_load\n )\n )\n )\n )\n```\n:::\n\n\nOK, seems we cannot get machine readable historical air photo information from the downloaded from the BC data catalogue [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-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 (@fig-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. \n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap <- ggplot() +\n geom_sf(\n data = layers_trimmed$aoi,\n fill = \"transparent\",\n color = \"black\",\n linewidth = .5\n ) +\n geom_sf(\n data = layers_trimmed$l_streams,\n color = \"blue\",\n size = 1\n ) +\n geom_sf(\n data = layers_trimmed$l_rail,\n color = \"black\",\n size = 1\n ) +\n geom_sf(\n data = layers_trimmed$l_imagery_hist,\n color = \"red\",\n size = 2\n ) +\n geom_sf(\n data = layers_trimmed$l_imagery_grid,\n alpha = 0.25,\n ) +\n geom_sf_text(\n data = layers_trimmed$l_imagery_grid,\n aes(label = map_tile),\n size = 3 # Adjust size of the text labels as needed\n )\n\nmap +\n geom_sf_text(\n data = layers_trimmed$l_streams |> dplyr::distinct(gnis_name, .keep_all = TRUE),\n aes(\n label = gnis_name\n ),\n size = 2 # Adjust size of the text labels as needed\n ) \n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n#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.\n\nlayers_trimmed$l_imagery_hist |> \n sf::st_drop_geometry() |> \n knitr::kable()\n```\n\n::: {.cell-output-display}\n\n\n|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 |\n|:---------------------------------------------------------|-----------------------:|:--------------|:------------------|:--------|----------:|--------:|:-----------------------------------------------------------------------------------|:------------------------------------------------------------------------------------------------------------------------------------------|--------:|:----------------|--------:|:-----------------|:-------------|:------------|\n|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 |\n|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 |\n|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 |\n|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 |\n|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 |\n|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 |\n|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 |\n|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 |\n\n\n:::\n:::\n\n \n
\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmy_caption <- \"Screenshot of kml downloaded from link provided in Historic Imagery Points.\"\nknitr::include_graphics(fs::path(\n path_post,\n \"fig\",\n \"Screenshot1\",\n ext = \"png\"\n )\n)\n```\n\n::: {.cell-output-display}\n{#fig-1 width=1492}\n:::\n:::\n\n\n\n\n
\n \nFor the [Historic Imagery Polygons](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-polygons) layer the range of `year_operational` is 1996, 2013. This is not as far back as we would prefer to be looking.\n \n\n
\n\nIt does however seem that each of the [Air Photo Centroids](https://catalogue.data.gov.bc.ca/dataset/airphoto-centroids) are georeferenced with a date range of:\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrange(layers_trimmed$l_photo_centroids$photo_date)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] \"1963-08-07\" \"2019-09-18\"\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n# 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.\n# list csvs\nls <- fs::dir_ls(\n fs::path(\n path_post,\n \"data\"),\n glob = \"*.csv\"\n)\n\nphotos_raw <- ls |> \n purrr::map_df(\n readr::read_csv\n ) |> \n sf::st_as_sf(\n coords = c(\"Longitude\", \"Latitude\"), crs = 4326\n ) |> \n janitor::clean_names() |> \n dplyr::mutate(photo_date = lubridate::mdy(photo_date)) \n\n\nphotos_aoi <- sf::st_intersection(\n photos_raw, \n layers_trimmed$aoi |> st_make_valid()\n )\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = layers_trimmed$l_photo_centroids,\n alpha = 0.25\n ) \n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\nThat is a lot of photos! 7910 photos to be exact!!!\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# amount to buffer all stream segments\nq_buffer <- 500\n# q_drm_main <- 263795\n\n# length of streams other than selected explicity to buffer\nq_drm_other <- 3000\n```\n:::\n\n\n## Clip Photo Information with Streams Buffers\nHere are our query parameters to narrow down the area within our study are watershed in which we want to find photos for:\n\n - Buffer: 500m - size of buffer used on either side of stream lines selected\n - Stream segments: \n + Bulkley River (`gnis_name` in the stream layer)\n + Maxan Creek\n + Buck Creek\n + for each remaining stream - segments of that stream which begin before 3000m from the downstream system (i.e. the first 3km) of stream.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# 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.\nknitr::include_graphics(fs::path(\n path_post,\n \"fig\",\n \"Screenshot2\",\n ext = \"png\"\n )\n)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nr_streams <- c(\"Maxan Creek\", \"Buck Creek\")\n\naoi_refined_raw <- layers_trimmed$l_streams |> \n # removed & downstream_route_measure < q_drm_main for bulkley as doestn't cahnge 1960s query and increases beyond just by 5 photos\n dplyr::filter(gnis_name == \"Bulkley River\"|\n gnis_name != \"Bulkley River\" & downstream_route_measure < q_drm_other |\n gnis_name %in% r_streams) |> \n # dplyr::arrange(downstream_route_measure) |>\n # calculate when we get to length_m by adding up the length_metre field and filtering out everything up to length_m\n # dplyr::filter(cumsum(length_metre) <= length_m) |>\n sf::st_union() |> \n # we need to run st_sf or we get a sp object in a list...\n sf::st_sf()\n \naoi_refined_buffered <- sf::st_buffer(\n aoi_refined_raw,\n q_buffer, endCapStyle = \"FLAT\"\n) \n\nphotos_aoi_refined <- sf::st_intersection(\n layers_trimmed$l_photo_centroids, \n aoi_refined_buffered\n )\n```\n:::\n\n\nLet'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.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = aoi_refined_buffered,\n color = \"red\",\n alpha= 0\n ) +\n geom_sf(\n data = photos_aoi_refined,\n alpha = 0.25,\n ) +\n geom_sf_text(\n data = layers_trimmed$l_streams |> dplyr::distinct(gnis_name, .keep_all = TRUE),\n aes(\n label = gnis_name\n ),\n size = 2 # Adjust size of the text labels as needed\n ) \n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\nThat is not as many photos - but still quite a few (1361). @fig-dt1 below can be used to filter these photos from any time and/or mapsheet and export the result to csv or excel file. \n\n\n::: {.cell tbl-cap='All photo centroids located with watershed study area.'}\n\n```{.r .cell-code}\nphotos_aoi_refined |> \n dplyr::select(-id) |> \n my_dt_table(cols_freeze_left = 0)\n```\n\n::: {#fig-dt1 .cell-output-display}\n\n```{=html}\n\n\n```\n\n:::\n:::\n\n\n## Filter Photos by Date\nNow lets map by year to see what our options are including the earliest photos possible. Here is our range to choose from:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrange(photos_aoi_refined$photo_date)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] \"1963-08-07\" \"2019-09-18\"\n```\n\n\n:::\n:::\n\n`\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\ngeom_sf(\n data = photos_aoi_refined |> dplyr::filter(photo_year <= \"1975\")\n ) +\n facet_wrap(~ photo_year)\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\nWell - 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...\n\n
\n\nThinking the ideal thing to do is to grab the photos from:\n\n - 1968 all\n - 1971 for the Buck Creek and Maxan Creek areas only\n - 1975 Maxan Creek only\n\n
\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# spatially represent just Buck and Maxan, buffer and clip the 1971 photos\n# \"r_\" is for \"refine\"\nr_year1 <- \"1968\"\nr_year2 <- \"1971\"\nr_year3 <- \"1975\"\n\nr_streams2 <- c(\"Maxan Creek\")\n\nl_streams_refined1 <- layers_trimmed$l_streams |> \n # we defined r_streams in chunk way above \n dplyr::filter(gnis_name %in% r_streams) |> \n sf::st_union() |> \n # we need to run st_sf or we get a sp object in a list...\n sf::st_sf()\n \naoi_refined_buffered2 <- sf::st_buffer(\n l_streams_refined1,\n q_buffer, endCapStyle = \"FLAT\"\n) \n\nl_streams_refined2 <- layers_trimmed$l_streams |> \n # we defined r_streams in chunk way above \n dplyr::filter(gnis_name %in% r_streams2) |> \n sf::st_union() |> \n # we need to run st_sf or we get a sp object in a list...\n sf::st_sf()\n \naoi_refined_buffered3 <- sf::st_buffer(\n l_streams_refined2,\n q_buffer, endCapStyle = \"FLAT\"\n) \n\n# filter first year\nphotos1 <- photos_aoi_refined |> \n dplyr::filter(\n photo_year == r_year1\n )\n\n# filter second year using just the streams we want to include\nphotos2 <- sf::st_intersection(\n layers_trimmed$l_photo_centroids |> dplyr::filter(photo_year == r_year2), \n aoi_refined_buffered2\n )\n\n# filter second year using just the streams we want to include\nphotos3 <- sf::st_intersection(\n layers_trimmed$l_photo_centroids |> dplyr::filter(photo_year == r_year3), \n aoi_refined_buffered3\n )\n\nphotos_all <- dplyr::bind_rows(photos1, photos2, photos3)\n```\n:::\n\n\n\nNow let's have a look at the individual year components (@fig-map5) as well as the whole dataset (@fig-map6). We are privileged to potentially have the assistance of [Mike Price](https://www.researchgate.net/profile/Michael-Price-14) 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.... \n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmy_caption <- \"Amalgamated photo points presented by year.\"\nmap +\n geom_sf(\n data = photos_all\n ) +\n facet_wrap(~ photo_year)\n```\n\n::: {.cell-output-display}\n{#fig-map5 width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nmy_caption <- \"Amalgamated photo points\"\nmap +\n geom_sf(\n data = photos_all\n ) +\n geom_sf_text(\n data = layers_trimmed$l_streams |> dplyr::distinct(gnis_name, .keep_all = TRUE),\n aes(\n label = gnis_name\n ),\n size = 2 # Adjust size of the text labels as needed\n ) \n```\n\n::: {.cell-output-display}\n{#fig-map6 width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nphotos_all |> \n dplyr::select(-id) |> \n my_dt_table(cols_freeze_left = 0)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n\n\n```\n\n:::\n:::\n\n\n\n## Export `csv` with Photo Information\nLet's burn out csv that can be used to find the imagery for the 136 photos above.\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlfile_name_photos <- function(dat = NULL){\n fs::path(\n path_post,\n \"exports\",\n paste(\n \"airphotos\",\n paste(range(dat$photo_date), collapse = \"_\"),\n sep = \"_\"\n ),\n ext = \"csv\"\n )\n}\n\nphotos_all |> \n readr::write_csv(\n lfile_name_photos(photos_all), na =\"\"\n )\n\n\nlpath_link <- function(dat = NULL){\n paste0(\n \"https://github.com/NewGraphEnvironment/new_graphiti/tree/main/posts/2024-11-15-bcdata-ortho-historic/exports/\",\n basename(\n lfile_name_photos(dat)\n )\n )\n}\n```\n:::\n\n\nWe can view and download exported csv files [here](https://github.com/NewGraphEnvironment/new_graphiti/tree/main/posts/2024-11-15-bcdata-ortho-historic/exports/). \n\n\n\n\n",
"supporting": [
"index_files"
],
diff --git a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/fig-map5-1.png b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/fig-map5-1.png
new file mode 100644
index 0000000..c1f81f8
Binary files /dev/null and b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/fig-map5-1.png differ
diff --git a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/fig-map6-1.png b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/fig-map6-1.png
new file mode 100644
index 0000000..a1acb8a
Binary files /dev/null and b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/fig-map6-1.png differ
diff --git a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map1-1.png b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map1-1.png
index 9163181..8421bef 100644
Binary files a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map1-1.png and b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map1-1.png differ
diff --git a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map2-1.png b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map2-1.png
index 88ea1c8..03bc72d 100644
Binary files a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map2-1.png and b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map2-1.png differ
diff --git a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map3-1.png b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map3-1.png
index 7ab0e87..d71156d 100644
Binary files a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map3-1.png and b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map3-1.png differ
diff --git a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map4-1.png b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map4-1.png
index 0b3b5c3..adbd01b 100644
Binary files a/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map4-1.png and b/_freeze/posts/2024-11-15-bcdata-ortho-historic/index/figure-html/map4-1.png differ
diff --git a/posts/2024-11-15-bcdata-ortho-historic/index.qmd b/posts/2024-11-15-bcdata-ortho-historic/index.qmd
index d155194..9608dab 100644
--- a/posts/2024-11-15-bcdata-ortho-historic/index.qmd
+++ b/posts/2024-11-15-bcdata-ortho-historic/index.qmd
@@ -279,13 +279,6 @@ map <- ggplot() +
color = "blue",
size = 1
) +
- 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
- ) +
geom_sf(
data = layers_trimmed$l_rail,
color = "black",
@@ -306,7 +299,14 @@ map <- ggplot() +
size = 3 # Adjust size of the text labels as needed
)
-map
+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
+ )
```
```{r tab-imagery-points, eval=TRUE}
@@ -467,14 +467,24 @@ map +
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 (`r nrow(photos_aoi_refined)`). The table below can be used to filter these photos from any time and/or mapsheet and export the result to csv or excel file.
+That is not as many photos - but still quite a few (`r nrow(photos_aoi_refined)`). @fig-dt1 below can be used to filter these photos from any time and/or mapsheet and export the result to csv or excel file.
-```{r tab1}
+```{r }
+#| label: fig-dt1
+#| tbl-cap: "All photo centroids located with watershed study area."
photos_aoi_refined |>
- my_dt_table(cols_freeze_left = 2)
+ dplyr::select(-id) |>
+ my_dt_table(cols_freeze_left = 0)
```
## Filter Photos by Date
@@ -562,15 +572,45 @@ photos_all <- dplyr::bind_rows(photos1, photos2, photos3)
```
-Now let's have a look at what we have
-```{r map5}
+Now let's have a look at the individual year components (@fig-map5) as well as the whole dataset (@fig-map6). We are privileged to potentially have the assistance of [Mike Price](https://www.researchgate.net/profile/Michael-Price-14) 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....
+
+
+```{r fig-map5, fig.cap=my_caption}
+
+my_caption <- "Amalgamated photo points presented by year."
+map +
+ geom_sf(
+ data = photos_all
+ ) +
+ facet_wrap(~ photo_year)
+
+```
+
+```{r fig-map6, fig.cap=my_caption}
+
+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
+ )
```
+```{r dt2}
+
+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 `r nrow(photos_all)` photos above.