Giter Club home page Giter Club logo

Comments (11)

bevingtona avatar bevingtona commented on September 7, 2024 4

Some ideas arround grabbing the .dem resources and stitching them together using an sf object as a filter...

library(sf)
library(dplyr)
library(rvest)
library(bcdata)
library(bcmaps)

# DEFINE YOUR AREA OF INTEREST AOI
my_aoi <- bcmaps::bc_cities() %>% filter(NAME == "Atlin") %>% st_buffer(50000)

# DOWNLOAD AND MOSAIC FUNCTION
my_download_and_mosaic <- function(my_aoi, 
                                   tempfolder,
                                   destination_tif) {

  dir.create(tempfolder, showWarnings = F)
  
  # 1:250k TILES THAT INTERSECT MY_AOI
  my_tiles <- bcdc_query_geodata("c50803db-7645-4da8-bce5-1bc5b501a7b3") %>% 
    dplyr::filter(INTERSECTS(my_aoi)) %>% 
    collect()
  print(paste("There area a total of", nrow(my_tiles), "1:250k tiles"))
  
  # TILE LIST
  my_tile_list <- my_tiles$MAP_TILE_DISPLAY_NAME
  
  # GET LINKS PER MAP SHEET
  lapply(my_tile_list, function(tilename){
    # tilename=my_tile_list[1]
    print(paste(tilename, "start"))
    
    # CREATE OUTPUT FOLDER
    dir.create(paste0(tempfolder, tilename), showWarnings = F)
    tile_files <- list.files(paste0(tempfolder, tilename))
    
    # MAP SHEET FOLDER
    URL <- paste0("https://pub.data.gov.bc.ca/datasets/175624/",tilename,"/")
    # READ PAGE
    PAGE <- html_attr(html_nodes(read_html(URL), "a"), "href")
    # LIST OF ZIP
    LIST <- PAGE[grep(pattern = "zip$", x = PAGE)]
    # URL OF ZIPS
    ZIPS <- paste0(URL,LIST)
    
    # DONWLOAD LIST
    lapply(1:length(LIST), function(i){
      if(sub(pattern = ".zip", replacement = "", LIST[i]) %in% tile_files){}else{
        print(paste(tilename, "start", i, "of", length(LIST)))
        download.file(ZIPS[i], destfile = paste0(tempfolder, tilename, "/", LIST[i]), quiet = T)
        unzip(paste0(tempfolder,tilename, "/", LIST[i]), exdir = paste0(tempfolder,tilename))
      }})
    
    # MAKE VRT OF 1:250k SUBTILES
    gdalUtils::gdalbuildvrt(gdalfile = list.files(path = paste0(tempfolder,tilename), pattern = "*.dem$", full.names = T), 
                            output.vrt = paste0(tempfolder,tilename,".vrt"), separate = 0)
    
  })
    
  
  # VRT of ALL VRTS  
  gdalUtils::gdalbuildvrt(gdalfile = list.files(path = paste0(tempfolder), pattern = "*.vrt$", full.names = T, recursive = F), 
                          output.vrt = paste0(tempfolder,"output.vrt"), separate = 0)
  # VRT TO TIF
  gdalUtils::gdal_translate(src_dataset = paste0(tempfolder,"output.vrt"), 
                            dst_dataset = paste0(tempfolder,destination_tif), ot = "Int16")
}  

# RUN
my_download_and_mosaic(my_aoi = my_aoi, 
                       tempfolder = "C:/Users/bevington/Downloads/trim/",
                       destination_tif = "Atlin.tif")    

stars::read_stars("C:/Users/bevington/Downloads/trim/Atlin.tif") %>% plot()

from bcmaps.

bevingtona avatar bevingtona commented on September 7, 2024

Could add a clip to aoi parameter, and a delete intermediate data parameter

from bcmaps.

bevingtona avatar bevingtona commented on September 7, 2024

Also worth noting: getSpatialData [https://github.com/16EAGLE/getSpatialData] has nice functionality to grab SRTM.

from bcmaps.

boshek avatar boshek commented on September 7, 2024

High level approach (function provisionally called cded):

from bcmaps.

ateucher avatar ateucher commented on September 7, 2024

@boshek you are going to create a new branch and start the PR with the mapsheets_250 function, correct? Then @bevingtona, @dfilatow and @HunterGleason can start on the cded_*functions?

Also, it was discussed that it would be nice for the core cded function to be user-facing so that it can return the path to the *.vrt file that a user could then use in other applications. It could have a path arg that defaults to tempfile(filext = ".vrt")...

from bcmaps.

bevingtona avatar bevingtona commented on September 7, 2024

We might want to keep some options for where the vrt gets stored .. For example, I would likely also like to use the vrt generated by the cded function in QGIS.

One thought is that the vrt would get created by cded, then loaded into R using cded_raster or cded_stars and then the vrt would be deleted... but I like vrts!

Perhaps an option to export the save the vrt file in a specific location could work... like:
my_raster <- cded_stars(aoi, clip, export, out_location, out_name, out_format)
where:
aoi = sf or sp polygon
clip = logical
export = logical, if false this function will load the mosaic (vrt) generated by the cded function as a stars object
out_location = path
out_name = character
out_format = "tif", "vrt", etc.

from bcmaps.

ateucher avatar ateucher commented on September 7, 2024

In terms of consistency and ease of maintenance, I would suggest adding those extra arguments to the core cded function. Then cded_stars and cded_raster could have a ... argument where the user could pass those through to cded.

I'm not sure if we want the function to do the clipping, or allow the user to do that on their own? What happens if export = TRUE?

I think most of what you mention could be accomplished by something like this?

cded <- function(aoi, path = tempfile(fileext = ".vrt") {
  # do stuff
  # write to "path" in format based on file extension 
  return(path)
}

cded_stars <- function(aoi, ...) {
  out_file <- cded(aoi, ...)
  read_stars(out_file)
}

cded_raster <- function(aoi, ...) {
  out_file <- cded(aoi, ...)
  read_raster(out_file)
}

from bcmaps.

boshek avatar boshek commented on September 7, 2024

I wonder if offloading any writing to individuals packages is a bit easier to maintain. So for example you can imagine:

cded_stars(aoi = aoi_object) %>%
   write_stars("filename.tif")

Then cded only returns the location of the .vrt?

from bcmaps.

HunterGleason avatar HunterGleason commented on September 7, 2024

Should the function download the tiles and store them in a permanent location as they are prompted for. Then there will be a single vrt in that directory that represents of all the tiles in the directory, i.e., a mosaic. This would be updated every time a new tile gets downloaded. The cded function would read in the vrt as a STARS or raster object, clip to the 'tiles' sf object, and return the clipped version of the vrt to the R environment, as either STARS or raster. Then no need to write a version of the clipped aoi to disk? If someone then wanted the vrt it would always be in the same location, and steadily grow as the user downloads more tiles? I guess the whole vrt is loaded to memory tho?

from bcmaps.

bevingtona avatar bevingtona commented on September 7, 2024

Added some comments, and ran the function for all of BC...! ~10 minutes without overviews

library(sf)
library(dplyr)
library(rvest)
library(bcdata)
library(bcmaps)

#### SET: CACHE FOLDER ####
dir = "C:/Users/bevin/BCMAPSCACHE_DEM/"
dir.create(dir)

#### SET: PROJECT FOLDER ####
dir_proj = "C:/Users/bevin/Documents"
filename = "BC_DEM_MOSAIC"

#### DEFINE YOUR AREA OF INTEREST AOI ####
my_aoi <- bcmaps::bc_cities() %>% filter(NAME == "Vancouver") %>% st_buffer(5000)
my_aoi <- bcmaps::bc_bound()

##### DEFINE DOWNLOAD AND MOSAIC FUNCTION ####
my_download_and_mosaic <- function(my_aoi, 
                                   tempfolder,
                                   destination_folder,
                                   destination_name) {
  
  # CREATE CACHE DIRECTORY
  dir.create(tempfolder, showWarnings = F)
  
  # TILES THAT INTERSECT MY_AOI
  my_tiles <- bcdc_query_geodata("c50803db-7645-4da8-bce5-1bc5b501a7b3") %>% 
    dplyr::filter(INTERSECTS(my_aoi)) %>% 
    collect()
    print(paste("There area a total of", nrow(my_tiles), "1:250k tiles"))
  
  # GET LIST OF TILE NAMES
  my_tile_list <- my_tiles$MAP_TILE_DISPLAY_NAME
  
  # GET LINKS PER MAP SHEET
  lapply(1:length(my_tile_list), function(j){
    # PRINT PROGRESS
    print(paste(round(100*j/nrow(my_tiles),0),"%"))    
    # GET TILENAME OF INTEREST
    tilename=my_tile_list[j]
    # CREATE OUTPUT FOLDER FOR TILE
    dir.create(paste0(tempfolder, tilename), showWarnings = F)
    # MAKE LIST OF EXISTING TILES IN 
    tile_files <- list.files(paste0(tempfolder, tilename))
    
    #### HTTP STUFF ####
    
    # MAP SHEET FOLDER
    URL <- paste0("https://pub.data.gov.bc.ca/datasets/175624/",tilename,"/")
    # READ PAGE
    PAGE <- html_attr(html_nodes(read_html(URL), "a"), "href")
    # LIST OF ZIP
    LIST <- PAGE[grep(pattern = "zip$", x = PAGE)]
    # LIST OF ZIP URLS
    ZIPS <- paste0(URL,LIST)
    
    #### DOWNLOAD MAGIC ####
    
    # DONWLOAD LIST
    lapply(1:length(LIST), function(i){
      if(sub(pattern = ".zip", replacement = "", LIST[i]) %in% tile_files){
        # print(paste(tilename, "[", j, "of", nrow(my_tiles), "250k tiles ]", "start", i, "of", length(LIST), "files", "ALREADY EXISTS"))
        }else{
          # print(paste(tilename, "[", j, "of", nrow(my_tiles), "250k tiles ]", "start", i, "of", length(LIST), "files"))
          download.file(ZIPS[i], destfile = paste0(tempfolder, tilename, "/", LIST[i]), quiet = T)
          unzip(paste0(tempfolder,tilename, "/", LIST[i]), exdir = paste0(tempfolder,tilename))
        }})
    
    # MAKE VRT PER 250k TILES
    gdalUtils::gdalbuildvrt(gdalfile = list.files(path = paste0(tempfolder,tilename), pattern = "*.dem$", full.names = T),
                            output.vrt = paste0(tempfolder,tilename,".vrt"), separate = 0)
    
  })
  
  # VRT of ALL VRTS  
  gdalUtils::gdalbuildvrt(gdalfile = list.files(path = paste0(tempfolder), pattern = "*.vrt$", full.names = T, recursive = F), 
                          output.vrt = paste0(destination_folder,"/",destination_name, ".vrt"), separate = 0, overwrite = T)
  
}  

#### OPTIONAL: BUILD OVERVIEWS FOR SPEED IN QGIS / ARCGIS #### 
# Takes time, but very helpful for pan/zoom speed 
#gdalUtils::gdaladdo(paste0(dir_proj,"/",filename, ".vrt"), levels = c(2,8,16))

#### RUN FUNCTION ####
system.time(my_download_and_mosaic(my_aoi = my_aoi,
                                   tempfolder = dir,
                                   destination_folder = dir_proj,
                                   destination_name = filename))
#> [1] "There area a total of 89 1:250k tiles"
#> [1] "1 %"
#> [1] "2 %"
#> [1] "3 %"
#> [1] "4 %"
#> [1] "6 %"
#> [1] "7 %"
#> [1] "8 %"
#> [1] "9 %"
#> [1] "10 %"
#> [1] "11 %"
#> [1] "12 %"
#> [1] "13 %"
#> [1] "15 %"
#> [1] "16 %"
#> [1] "17 %"
#> [1] "18 %"
#> [1] "19 %"
#> [1] "20 %"
#> [1] "21 %"
#> [1] "22 %"
#> [1] "24 %"
#> [1] "25 %"
#> [1] "26 %"
#> [1] "27 %"
#> [1] "28 %"
#> [1] "29 %"
#> [1] "30 %"
#> [1] "31 %"
#> [1] "33 %"
#> [1] "34 %"
#> [1] "35 %"
#> [1] "36 %"
#> [1] "37 %"
#> [1] "38 %"
#> [1] "39 %"
#> [1] "40 %"
#> [1] "42 %"
#> [1] "43 %"
#> [1] "44 %"
#> [1] "45 %"
#> [1] "46 %"
#> [1] "47 %"
#> [1] "48 %"
#> [1] "49 %"
#> [1] "51 %"
#> [1] "52 %"
#> [1] "53 %"
#> [1] "54 %"
#> [1] "55 %"
#> [1] "56 %"
#> [1] "57 %"
#> [1] "58 %"
#> [1] "60 %"
#> [1] "61 %"
#> [1] "62 %"
#> [1] "63 %"
#> [1] "64 %"
#> [1] "65 %"
#> [1] "66 %"
#> [1] "67 %"
#> [1] "69 %"
#> [1] "70 %"
#> [1] "71 %"
#> [1] "72 %"
#> [1] "73 %"
#> [1] "74 %"
#> [1] "75 %"
#> [1] "76 %"
#> [1] "78 %"
#> [1] "79 %"
#> [1] "80 %"
#> [1] "81 %"
#> [1] "82 %"
#> [1] "83 %"
#> [1] "84 %"
#> [1] "85 %"
#> [1] "87 %"
#> [1] "88 %"
#> [1] "89 %"
#> [1] "90 %"
#> [1] "91 %"
#> [1] "92 %"
#> [1] "93 %"
#> [1] "94 %"
#> [1] "96 %"
#> [1] "97 %"
#> [1] "98 %"
#> [1] "99 %"
#> [1] "100 %"
#>    user  system elapsed 
#>    4.09    7.69   45.00

#### READ OUTPUT RASTER ####
system.time(dem <- stars::read_stars(paste0(dir_proj,"/",filename,".vrt"), proxy = T))
#>    user  system elapsed 
#>    0.03    0.06    0.10

#### PLOT OUTPUT RASTER ####
system.time(plot(dem))

#>    user  system elapsed 
#>  289.20   13.33  302.53

Created on 2020-10-28 by the reprex package (v0.3.0)

from bcmaps.

bevingtona avatar bevingtona commented on September 7, 2024

New changes:

  1. gdalUtils:: changed to sf::gdal_utils
  2. The script deletes the .zip and converts the .dem to .tif to save ~45% space
  3. Name changes to clarify what is "cache" and what is "project"
library(sf)
library(dplyr)
library(rvest)
library(bcdata)
library(bcmaps)

#### SET: CACHE FOLDER ####
cachedir = "C:/Users/bevin/BCMAPSCACHE_DEM2/"
dir.create(cachedir)

#### SET: PROJECT FOLDER ####
projdir = "C:/Users/bevin/MyProjectFolder/"
dir.create(projdir)
filename = "PG_DEM_MOSAIC"

#### DEFINE YOUR AREA OF INTEREST AOI ####
my_aoi <- bcmaps::bc_cities() %>% filter(NAME == "Prince George") %>% st_buffer(5000)

##### DEFINE DOWNLOAD AND MOSAIC FUNCTION ####
cded <- function(my_aoi, 
                 cachefolder,
                 projectfolder,
                 exportname) {
  
  # CREATE CACHE DIRECTORY
  dir.create(cachefolder, showWarnings = F)
  # TILES THAT INTERSECT MY_AOI
  my_tiles <- bcdc_query_geodata("c50803db-7645-4da8-bce5-1bc5b501a7b3") %>% 
    dplyr::filter(INTERSECTS(my_aoi)) %>%
    collect()
  write_sf(my_tiles, "C:/Users/bevin/250.sqlite")
  print(paste("There area a total of", nrow(my_tiles), "250k tiles"))
  
  # GET LIST OF TILE NAMES
  my_tile_list <- my_tiles$MAP_TILE_DISPLAY_NAME
  
  # GET LINKS PER MAP SHEET
  vrtlist <- lapply(1:length(my_tile_list), function(j){
    # GET TILENAME OF INTEREST
    tilename=my_tile_list[j]
    # CREATE OUTPUT FOLDER FOR TILE
    dir.create(paste0(cachefolder, tilename), showWarnings = F)
    # MAKE LIST OF EXISTING TILES IN 
    tile_files <- list.files(paste0(cachefolder, tilename))
    
    #### HTTP STUFF ####
    
    # MAP SHEET FOLDER
    URL <- paste0("https://pub.data.gov.bc.ca/datasets/175624/",tilename,"/")
    # READ PAGE
    PAGE <- html_attr(html_nodes(read_html(URL), "a"), "href")
    # LIST OF ZIP
    LIST <- PAGE[grep(pattern = "dem.zip$", x = PAGE)]
    # LIST OF ZIP URLS
    ZIPS <- paste0(URL,LIST)
    
    #### DOWNLOAD MAGIC ####
    
    # DONWLOAD LIST
    subvrtlist <- lapply(1:length(LIST), function(i){
      if(sub(pattern = ".dem.zip", replacement = ".tif", LIST[i]) %in% tile_files){
        print(paste(tilename, "[", j, "of", nrow(my_tiles), "250k tiles ]", "start", i, "of", length(LIST), "files", "ALREADY EXISTS"))
        tile <- sub(".dem", ".tif", sub(".zip", "", paste0(cachefolder,tilename, "/", LIST[i])))
        tile
        }else{
        print(paste(tilename, "[", j, "of", nrow(my_tiles), "250k tiles ]", "start", i, "of", length(LIST), "files", "DOWNLOADING"))
        # Download
        download.file(ZIPS[i], destfile = paste0(cachefolder, tilename, "/", LIST[i]), quiet = T)
        # UNZIP
        unzip(paste0(cachefolder,tilename, "/", LIST[i]), exdir = paste0(cachefolder,tilename))
        # REMOVE ZIP
        file.remove(paste0(cachefolder,tilename, "/", LIST[i]))
        # TRANSLATE .DEM to .TIF
        sf::gdal_utils(util = "translate", 
                       source = sub(".zip", "", paste0(cachefolder,tilename, "/", LIST[i])), 
                       destination = sub(".dem", ".tif", sub(".zip", "", paste0(cachefolder,tilename, "/", LIST[i]))), 
                       options = c("-ot","Int16", 
                                   "-of", "GTiff"))
        # REMOVE .DEM
        file.remove(sub(".zip", "", paste0(cachefolder,tilename, "/", LIST[i])))
        
        tile <- sub(".dem", ".tif", sub(".zip", "", paste0(cachefolder,tilename, "/", LIST[i])))
        tile
      }})
    
    print(paste(round(100*j/nrow(my_tiles),0),"%"))    
    
    do.call(rbind, subvrtlist)
    })
  
  vrtlist <- do.call(rbind, vrtlist)

  # VRT of ALL VRTS  
  sf::gdal_utils(util = "buildvrt", 
                 source = vrtlist[,1], #list.files(path = paste0(cachefolder), pattern = "*.vrt$", full.names = T, recursive = F), 
                 destination = paste0(projectfolder,exportname, ".vrt"))
  
}  

#### OPTIONAL: BUILD OVERVIEWS FOR SPEED IN QGIS / ARCGIS #### 
# Takes time, but very helpful for pan/zoom speed 
#gdalUtils::gdaladdo(paste0(dir_proj,"/",filename, ".vrt"), levels = c(2,8,16))

#### RUN FUNCTION ####
system.time(cded(my_aoi = my_aoi,
                 cachefolder = cachedir,
                 projectfolder = projdir,
                 exportname = filename))


#### READ OUTPUT RASTER ####
system.time(dem <- stars::read_stars(paste0(projdir,"/",filename,".vrt"), proxy = T))

#### PLOT OUTPUT RASTER ####
system.time(plot(dem))

from bcmaps.

Related Issues (20)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.