From b385623dd24f354265d41457c170fc63ae7b235a Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 17 Nov 2024 10:46:12 -1000 Subject: [PATCH] correct spelling for WIMSI --- .../index/execute-results/html.json | 4 ++-- posts/2024-11-15-bcdata-ortho-historic/index.qmd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) 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 32e4e4c..f3a9e3d 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": "827cf180fd306465edd702eb4b5a8bda", + "hash": "1126034641864aacb56bbe09583fa03d", "result": { "engine": "knitr", - "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
\n\n \nThe table below can be used to filter photos from any time (provided it's location falls in the mapsheet areas we downloaded info for) and export the result to csv or excel file. This could be useful for narrowing down our search not only to timescales but also to specific mapsheets with options given for both NTS 1:50,000 as well as BCGS 1:20,000 grid mapsheets.\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## Export `csv` Files with Photo Information\nLet's burn out some csvs that can be used to find the imagery for:\n\n - \"before 1969-12-31\" \n - 1970-01-01 to 1975-12-31 timeline\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\n# path_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\nphotos1 |> \n readr::write_csv(\n lfile_name_photos(photos1), na =\"\"\n )\n\nphotos2 |> \n readr::write_csv(\n lfile_name_photos(photos2), na =\"\"\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## Pick a Horse\nWe are absolutely privaledged to potentially have the assistance of Mike Price to help us obtain this imagery from the UBC archives. Sooo... - unless he has a preference that 100% can override our recommendation - we would think that he would probably prefer clear direction from us. Soooo - to start - thinking the best plan is to go with the earliest imagery possible and we scan in the \"before 1969-12-31\" timeframe. The name of the csv with that information is `airphotos_1968-05-09_1968-07-31.csv` and the file can found along with the other export at the link above or individually [here](https://github.com/NewGraphEnvironment/new_graphiti/tree/main/posts/2024-11-15-bcdata-ortho-historic/exports/airphotos_1968-05-09_1968-07-31.csv).\n\n## Outstanding Issues\n\nWhy does the date range of the photos we can get IDs for through the [WIMSIN](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/digital-imagery/air-photos) online portal not overlap the dates listed in the BC Data Catalouge [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points) for the dates in the Historic Imagery Points that are listed as 1950 - 1963 vs our earliest records from the portal at 1968-05-09??!!\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 [WIMSI](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/digital-imagery/air-photos) 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
\n\n \nThe table below can be used to filter photos from any time (provided it's location falls in the mapsheet areas we downloaded info for) and export the result to csv or excel file. This could be useful for narrowing down our search not only to timescales but also to specific mapsheets with options given for both NTS 1:50,000 as well as BCGS 1:20,000 grid mapsheets.\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## Export `csv` Files with Photo Information\nLet's burn out some csvs that can be used to find the imagery for:\n\n - \"before 1969-12-31\" \n - 1970-01-01 to 1975-12-31 timeline\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\n# path_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\nphotos1 |> \n readr::write_csv(\n lfile_name_photos(photos1), na =\"\"\n )\n\nphotos2 |> \n readr::write_csv(\n lfile_name_photos(photos2), na =\"\"\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## Pick a Horse\nWe are absolutely privaledged to potentially have the assistance of Mike Price to help us obtain this imagery from the UBC archives. Sooo... - unless he has a preference that 100% can override our recommendation - we would think that he would probably prefer clear direction from us. Soooo - to start - thinking the best plan is to go with the earliest imagery possible and we scan in the \"before 1969-12-31\" timeframe. The name of the csv with that information is `airphotos_1968-05-09_1968-07-31.csv` and the file can found along with the other export at the link above or individually [here](https://github.com/NewGraphEnvironment/new_graphiti/tree/main/posts/2024-11-15-bcdata-ortho-historic/exports/airphotos_1968-05-09_1968-07-31.csv).\n\n## Outstanding Issues\n\nWhy does the date range of the photos we can get IDs for through the [WIMSI](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/digital-imagery/air-photos) online portal not overlap the dates listed in the BC Data Catalouge [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points) for the dates in the Historic Imagery Points that are listed as 1950 - 1963 vs our earliest records from the portal at 1968-05-09??!!\n\n", "supporting": [ "index_files" ], diff --git a/posts/2024-11-15-bcdata-ortho-historic/index.qmd b/posts/2024-11-15-bcdata-ortho-historic/index.qmd index d5394d2..602d5e2 100644 --- a/posts/2024-11-15-bcdata-ortho-historic/index.qmd +++ b/posts/2024-11-15-bcdata-ortho-historic/index.qmd @@ -297,7 +297,7 @@ It does however seem that if we go through the online interface, we can use a NT ## Download Photo Information -Next, 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: +Next, we go into the provincial [WIMSI](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/digital-imagery/air-photos) 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: - `Years from`: 1963 to 1980 - `Media`: Black and White - Film @@ -557,5 +557,5 @@ We are absolutely privaledged to potentially have the assistance of Mike Price t ## Outstanding Issues -Why does the date range of the photos we can get IDs for through the [WIMSIN](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/digital-imagery/air-photos) online portal not overlap the dates listed in the BC Data Catalouge [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points) for the dates in the Historic Imagery Points that are listed as 1950 - 1963 vs our earliest records from the portal at `r min(photos_aoi_refined$photo_date)`??!! +Why does the date range of the photos we can get IDs for through the [WIMSI](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/digital-imagery/air-photos) online portal not overlap the dates listed in the BC Data Catalouge [Historic Imagery Points](https://catalogue.data.gov.bc.ca/dataset/airborne-imagery-historical-index-map-points) for the dates in the Historic Imagery Points that are listed as 1950 - 1963 vs our earliest records from the portal at `r min(photos_aoi_refined$photo_date)`??!!