Example data generation for earthpy

earthpy-data-gen.R

library(tidyverse)
library(raster)
library(sf)
library(rmapshaper)
library(elevatr)
library(USAboundaries)

dir.create('out')

# RGB satellite imagery over Rocky Mountain National Park -------------
download.file(url = 'https://eoimages.gsfc.nasa.gov/images/imagerecords/88000/88405/rmnp_oli_2014263_geo.tif', 
              destfile = 'rmnp_oli_2014263_geo.tif')
r <- stack('rmnp_oli_2014263_geo.tif')

# Coarsen and reproject to epsg 4326
pr <- projectRaster(r, 
                    crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs", 
                    res = 0.0015)
pr[] <- as.integer(pr[])
dataType(pr)

writeRaster(pr, 
            file.path('out', 'rmnp-rgb.tif'), 
            datatype = 'INT1U', 
            overwrite = TRUE)

# write separate R, G, and B bands
fnames <- c('red.tif', 'green.tif', 'blue.tif')
for (i in seq_along(fnames)) {
  writeRaster(pr[[i]], 
              file.path('out', fnames[i]), 
              datatype = 'INT1U', 
              overwrite = TRUE)
}

# RMNP park boundary ------------------------------------------------------
rmnp <- st_zm(st_read("https://opendata.arcgis.com/datasets/0b69abf463aa4a44b0ba45988fbba8af_0.geojson"))
rmnp <- ms_simplify(rmnp)
plot(rmnp, max.plot = 1)

# to avoid saving a multipolygon, split into separate polygon features, 
# and discard a little chunk of RMNP that will complicate examples
rmnp_poly <- st_cast(rmnp, "POLYGON")
bigger_poly <- which.max(st_area(rmnp_poly))
big_part <- rmnp_poly[bigger_poly, ]
names(big_part$geometry) <- '0'
st_write(big_part, file.path('out', 'rmnp.shp'))


# RMNP digital elevation model --------------------------------------------
rmnp_sp <- as(rmnp, 'Spatial')
dem <- get_elev_raster(rmnp_sp, z = 8, verbose = TRUE)
cdem <- crop(dem, rmnp_sp)
writeRaster(cdem, 
            filename = file.path('out', 'rmnp-dem.tif'), 
            overwrite = TRUE, 
            datatype = "INT2U")


# glacier point locations --------------------------------------------
download.file('http://www.glims.org/download/glims_db_20171027.zip', 
              'glims_db_20171027.zip')
unzip('glims_db_20171027.zip')
glaciers <- st_read(file.path('glims_download_07040', 'glims_points.shp'))
colorado <- st_read('https://opendata.arcgis.com/datasets/4402a8e032ed49eb8b37fd729e4e8f03_9.geojson')
co_glaciers <- st_crop(glaciers, colorado)
st_write(co_glaciers, file.path('out', 'colorado-glaciers.geojson'))


# Colorado counties ---------------------------------------------------
co_counties <- us_counties(resolution = "low", states = "CO")
split_counties <- st_cast(co_counties, "POLYGON")
st_write(split_counties, file.path('out', 'colorado-counties.geojson'))


# Colorado section of the continental divide trail ------------------------
download.file("https://www.dropbox.com/s/82vdy5jm392y2pp/CDT_Google_Earth_2017.kml?dl=1", 
              "CDT_Google_Earth_2017.kml")
cdt <- st_read("CDT_Google_Earth_2017.kml")
simple_cdt <- cdt %>%
  st_transform(crs = 4326) %>%
  ms_simplify(keep = 0.1) %>%
  dplyr::select(geometry) %>%
  mutate(name = 'Continental Divide Trail')
cdt_path <- file.path('out', 'continental-div-trail.geojson')
unlink(cdt_path)
st_write(simple_cdt, cdt_path)