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 ac0ab14..32e4e4c 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": "3ae8cf90e00a8ff3848917feb8e050a6", + "hash": "827cf180fd306465edd702eb4b5a8bda", "result": { "engine": "knitr", - "markdown": "---\ntitle: \"Getting details of historic orthophoto imagery with R\"\nauthor: \"al\"\ndate: \"2024-11-15\"\ndate-modified: \"2024-11-15\"\ncategories: [fwapg, r, bcdata]\nimage: \"image.jpg\"\nparams:\n repo_owner: \"NewGraphEnvironment\"\n repo_name: \"new_graphiti\"\n post_name: \"2024-11-15-bcdata-ortho-historic\"\nformat: \n html:\n code-fold: true\n---\n\n\nWe need some historic ortho photo imagery so that we can have a look at historic watershed conditions compared to current. First thing we're gonna do is generate a area of interest for which we want all of the map sheets ideas for.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsuppressMessages(library(tidyverse))\nlibrary(ggplot2)\nlibrary(bcdata)\nlibrary(fwapgr)\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:::\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#get the bounding box of our aoi\n# aoi_bb <- sf::st_bbox(aoi)\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::: {.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_to_trim <- tibble::lst(\n !!!mget(ls)\n)\n\n\n# Function to validate and repair geometries\nvalidate_geometries <- function(layer) {\n layer <- sf::st_make_valid(layer)\n layer <- layer[sf::st_is_valid(layer), ]\n return(layer)\n}\n\n# Apply validation to the AOI and layers\naoi <- validate_geometries(aoi)\nlayers_to_trim <- purrr::map(layers_to_trim, validate_geometries)\n\n# clip them with purrr and sf\nlayers_trimmed <- purrr::map(\n layers_to_trim,\n ~ sf::st_intersection(.x, aoi)\n) \n```\n:::\n\nOK, so it looks like we cannot get historical air photo information from layers downloaded from BC data catalogue because the majority of the photos are not referenced. It does however seem that if we go through the online interface, we can use a NTS 1 to 50,000 map sheet to export a CSV file that does contain information about individual air photo locations. Our study area is huge so to begin with let's just have a look at the area from Macquarie Creek up to a boat the location of Talley. For this area we just have 21 to 50,000 map sheets. We can visualize which map she's those are by downloading that map sheet layer from BC data catalogue overlay it on our general study area.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap <- ggplot() +\n # tidyterra::geom_spatraster(\n # data = as.factor(land_cover_raster_resampled),\n # use_coltab = TRUE,\n # maxcell = Inf\n # ) +\n geom_sf(\n data = 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 = 1\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 # ggdark::dark_theme_void()\n\n\nmap\n```\n\n::: {.cell-output-display}\n![](index_files/figure-html/unnamed-chunk-4-1.png){width=672}\n:::\n:::\n\n\n\nSo looks like we have map sheet 093L09 as well as 0 93L08. Next, we're gonna go into the provincial [WIMSI9](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/digital-imagery/air-photos) service and manually select our NTS grids, state the date range we would like to search and export. We chose the following options:\n\n - Years from: 1963 to 1980\n\nIt looks like this.\n\n\n::: {.cell}\n\n```{.r .cell-code}\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/Screenshot1.png){width=743}\n:::\n\n```{.r .cell-code}\n# knitr::include_graphics(\"fig/Screenshot 2024-11-15 at 2.03.42 PM.png\")\n```\n:::\n\n\nLet's read in what we have, turn into spatial object, trim to overall study area and plot\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# list csvs\n\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\nphotos_aoi <- sf::st_intersection(photos_raw, aoi)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = photos_aoi,\n alpha = 0.25,\n ) \n```\n\n::: {.cell-output-display}\n![](index_files/figure-html/map2-1.png){width=672}\n:::\n\n```{.r .cell-code}\n # geom_sf_text(\n # data = photos_aoi,\n # aes(\n # label = \n # paste(\n # # roll_frame_series,\n # frame,\n # sep = \"-\"\n # )\n # ),\n # size = 2 # Adjust size of the text labels as needed\n # )\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nbuffer <- 3000\ndrm <- 263795\n```\n:::\n\n\n\nOk - so let's say we only want the ones that are within say 3000 of either side of the Bulkley River mainstem.\n\n\nLet's grab just the mainstem segments up to the top end of Bulkley Lake, merge, buffer the stream and clip our airphoto list by that.\n\nWe can use the `downstream_route_measure` of the stream layer to exclude areas upstream of Bulkley Lake (also known as Taman Creek). It looks like this in QGIS\n\n\n::: {.cell}\n\n```{.r .cell-code}\nknitr::include_graphics(fs::path(\n path_post,\n \"fig\",\n \"Screenshot2\",\n ext = \"png\"\n )\n)\n```\n\n::: {.cell-output-display}\n![](fig/Screenshot2.png){width=946}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\naoi_refined_raw <- layers_trimmed$l_streams |> \n dplyr::filter(gnis_name == \"Bulkley River\" & downstream_route_measure < drm) |> \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\n\naoi_refined_buffered <- sf::st_buffer(\n # we need to run st_sf or we get a sp object in a list...\n sf::st_sf(\n aoi_refined_raw\n ), \n buffer, endCapStyle = \"FLAT\"\n) \n\nphotos_aoi_refined <- sf::st_intersection(photos_raw, aoi_refined_buffered)\n```\n:::\n\n\nok lets plot again\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = photos_aoi_refined,\n alpha = 0.25,\n ) \n```\n\n::: {.cell-output-display}\n![](index_files/figure-html/map3-1.png){width=672}\n:::\n:::\n\nNow lets filter by date to see what we can get for the earliest photos possible. What is our range:\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] \"1968-05-09\" \"1980-09-07\"\n```\n\n\n:::\n:::\n\n\nWhat do we have for photos before 1970\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = photos_aoi_refined |> dplyr::filter(photo_date <= \"1969-12-31\"),\n alpha = 0.25,\n ) \n```\n\n::: {.cell-output-display}\n![](index_files/figure-html/map4-1.png){width=672}\n:::\n:::\n\nLooks like decent coverage so lets go with that.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nphotos <- photos_aoi_refined |> \n dplyr::filter(photo_date <= \"1969-12-31\")\n\nphotos |> \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\nThere are 104 photos. Let's burn out a csv that can be used to find them.\n\n\n::: {.cell}\n\n```{.r .cell-code}\npath_output <- fs::path(\n path_post,\n \"exports\",\n paste(\n \"airphotos\",\n paste(range(photos$photo_date), collapse = \"_\"),\n sep = \"_\"\n ),\n ext = \"csv\"\n )\n\nphotos |> \n readr::write_csv(\n path_output, na =\"\"\n )\n```\n:::\n\n\n\nWe can view and download exported csv file [here](https://github.com/NewGraphEnvironment/new_graphiti/tree/main/posts/2024-11-15-bcdata-ortho-historic/exports/)\n\n\n::: {.cell}\n\n```{.r .cell-code}\n leaflet::leaflet() |>\n leaflet::addTiles() |>\n leafem:::addCOG(\n url = url\n , group = \"COG\"\n , resolution = 512\n , autozoom = TRUE\n )\n```\n:::\n", + "markdown": "---\ntitle: \"Getting details of historic orthophoto imagery with R\"\nauthor: \"al\"\ndate: \"2024-11-15\"\ndate-modified: \"2024-11-16\"\ncategories: [fwapg, r, bcdata]\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)\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 5th 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 \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# burn all these to file so quick to repeat\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\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\n\nOK, seems we cannot get machine readable historical air photo information from layers downloaded from BC data catalogue - perhpas because the majority of the photos are not georeferenced? What we see in the map under the 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. From as far as I can tell - if we wanted to try to use these kmls or pdfs to select orthoimagery we would need to manually 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. 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\n\n::: {.cell}\n\n```{.r .cell-code}\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::: {.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,\n # alpha = 0,\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 # ggdark::dark_theme_void()\n\n\nmap\n```\n\n::: {.cell-output-display}\n![](index_files/figure-html/unnamed-chunk-1-1.png){width=672}\n:::\n:::\n\n\nIt does however seem that if we go through the online interface, we can use a NTS 1:50,000 map sheet to export a CSV file that does contain information about individual air photo locations (centre coordinate?). Our study area is huge so to begin with we can look at the area from somewhere around the confluence of McQuarie Creek and the Neexdzii Kwah up to say the upstream end of Bulkley Lake. For this area we just have two 1:50,000 map sheets (`093L09` and `093L08`). \n\n\n\n## Download Photo Information\n\nNext, we go into the provincial [WIMSIN](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/digital-imagery/air-photos) (getting that abbreviation wrong but it is close and the link works) service and for each NTS 1:50,000 grid - we state the date range we would like to search and export. Thinking we may be able to hit the API directly with a query in the future which would eliminate the painful point and click portion of this workflow. For now though - we are doing this piece manually and chose the following options:\n\n - `Years from`: 1963 to 1980\n - `Media`: Black and White - Film\n - `Only show georeferenced images`: - No\n\nRunning these options exports a csv for each query. We saved them in this repo with their original sane names `airphotos.{mapsheet}.csv`.\n\n\n::: {.cell}\n\n```{.r .cell-code}\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/Screenshot1.png){width=743}\n:::\n:::\n\n\nAt 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\n\n::: {.cell}\n\n```{.r .cell-code}\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 = photos_aoi,\n alpha = 0.25,\n ) \n```\n\n::: {.cell-output-display}\n![](index_files/figure-html/map2-1.png){width=672}\n:::\n\n```{.r .cell-code}\n # geom_sf_text(\n # data = photos_aoi,\n # aes(\n # label = \n # paste(\n # # roll_frame_series,\n # frame,\n # sep = \"-\"\n # )\n # ),\n # size = 2 # Adjust size of the text labels as needed\n # )\n```\n:::\n\n\nThat is a lot of photos! 849 photos to be exact!!!\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nq_buffer <- 1000\nq_drm_main <- 263795\nq_drm_other <- 5000\n```\n:::\n\n\n## Clip Photo Information with Streams Buffers\nHere are our query parameters to narrow down the area within our selected mapsheets in which we want to find photos for:\n\n - Buffer: 1000m - 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 + stream segments which begin before 5000 from the mainstem\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}\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 == \"Maxan Creek\") |> \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 photos_raw, \n aoi_refined_buffered\n )\n```\n:::\n\n\nLet's plot again and include our buffered areas around the first 5000m of streams (area in red) along with the location of the photo points that land within that area. Looks like this give us 234 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![](index_files/figure-html/map3-1.png){width=672}\n:::\n:::\n\n\n## Filter Photos by Date\nNow lets filter by date 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] \"1968-05-09\" \"1980-09-07\"\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nq_date_end1 <- \"1969-12-31\"\nq_date_start1 <- \"1970-01-01\"\nq_date_end2 <- \"1975-12-31\"\n```\n:::\n\n\nThere is lots going on here so we will map a comparison of a couple of options. First we will look at all photos available for before 1970 (red dots on map below). We will also see what photos we have for between 1970-01-01 and on or before 1975-12-31 (blue dots on map below)?\n\n\n::: {.cell}\n\n```{.r .cell-code}\nphotos1 <- photos_aoi_refined |> \n dplyr::filter(\n photo_date <= q_date_end1\n )\n\nphotos2 <- photos_aoi_refined |> \n dplyr::filter(\n photo_date >= q_date_start1 &\n photo_date <= q_date_end2\n )\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nmap +\n geom_sf(\n data = photos1,\n color = \"red\",\n alpha = 0.25,\n ) \n```\n\n::: {.cell-output-display}\n![](index_files/figure-html/map4-1.png){width=672}\n:::\n\n```{.r .cell-code}\nmap + \n geom_sf(\n data = photos2,\n color = \"blue\",\n alpha = 0.25,\n ) \n```\n\n::: {.cell-output-display}\n![](index_files/figure-html/map4-2.png){width=672}\n:::\n:::\n\n\nWell we have some tough decisions to make as we seem to have better coverage in the Maxan Creek drainage in the 70s (91 photos) but good coverage of the mainstem Neexdzii Kwah in the 60s (72 photos)...\n\n