Sample Workflow for Converting .gdb to .tif
Cameryn Brock
1/24/2022
Attach libraries
library(landscapemetrics)
library(raster)
library(fasterize)
library(sf)
library(rasterVis)
library(tidyverse)Read in data
Define template raster
# need template raster for fasterize and to define crs up front
# I use a FORE-SCE one here but whatever you decide is best
template <- raster("~/Conservation International/Data/FORE-SCE/CONUS_Landcover_B2/CONUS_B2_y2006.tif")
template## class : RasterLayer
## dimensions : 11747, 18563, 218059561 (nrow, ncol, ncell)
## resolution : 250, 250 (x, y)
## extent : -2357953, 2282797, 238542.6, 3175293 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
## source : CONUS_B2_y2006.tif
## names : CONUS_B2_y2006
## values : 0, 255 (min, max)WBD from USGS
GDB downloaded from USGS: https://prd-tnm.s3.amazonaws.com/index.html?prefix=StagedProducts/Hydrography/WBD/National/GDB/
# can look at layers in gdb
st_layers(dsn = "WBD_National_GDB/WBD_National_GDB.gdb")## Driver: OpenFileGDB
## Available layers:
## layer_name geometry_type features fields
## 1 WBDHU12 Multi Polygon 102968 20
## 2 NWISDrainageLine Multi Line String 127876 12
## 3 WBDHU8 Multi Polygon 2400 15
## 4 NonContributingDrainageLine Multi Line String 0 10
## 5 WBDHU4 Multi Polygon 244 15
## 6 WBDHU16 Multi Polygon 7202 20
## 7 NonContributingDrainageArea Multi Polygon 0 12
## 8 NWISDrainageArea Multi Polygon 2733 19
## 9 WBDHU6 Multi Polygon 404 15
## 10 WBDHU2 Multi Polygon 22 15
## 11 WBDHU14 Multi Polygon 8024 20
## 12 WBDLine Multi Line String 394614 8
## 13 WBDHU10 Multi Polygon 18853 17
## 14 FeatureToMetadata NA 424657 3
## 15 ExternalCrosswalk NA 0 7
## 16 HUMod NA 23 5
## 17 MetaProcessDetail NA 1559 16
## 18 MetaSourceDetail NA 1562 27
## 19 ProcessingParameters NA 2 2
## 20 UpdateStatus NA 0 2# read in gdb with sf
# assuming we care about HU12
watersheds <- sf::st_read(
dsn = "WBD_National_GDB/WBD_National_GDB.gdb",
layer = "WBDHU12"
) %>%
st_transform(crs = crs(template))## Reading layer `WBDHU12' from data source
## `C:\Users\cbrock\OneDrive - Conservation International Foundation\Documents\Conservation International\Arnhold SPRCLER-CC\arnhold_fellows\WBD_National_GDB\WBD_National_GDB.gdb'
## using driver `OpenFileGDB'
## Simple feature collection with 102968 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -179.2297 ymin: -14.61019 xmax: 179.8567 ymax: 71.43957
## Geodetic CRS: NAD83# look at variables
glimpse(watersheds)## Rows: 102,968
## Columns: 21
## $ tnmid <chr> "{AAF0D733-828B-4B8E-9E52-388A49AC0A23}", "{F~
## $ metasourceid <chr> "{33EA2180-A425-4AE6-951F-ABCBBF25B893}", "{3~
## $ sourcedatadesc <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
## $ sourceoriginator <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
## $ sourcefeatureid <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
## $ loaddate <dttm> 2017-05-25 02:05:31, 2017-05-25 02:05:42, 20~
## $ referencegnis_ids <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
## $ areaacres <dbl> 13211.47, 10514.34, 15561.22, 31033.21, 28839~
## $ areasqkm <dbl> 53.46, 42.55, 62.97, 125.59, 116.71, 198.76, ~
## $ states <chr> "IA,MN", "MN", "MN", "GA", "GA", "GA", "MN", ~
## $ huc12 <chr> "070200090402", "070200030503", "070200030602~
## $ name <chr> "Judicial Ditch Number Thirteen", "West Branc~
## $ hutype <chr> "S", "S", "S", "S", "S", "S", "S", "S", "S", ~
## $ humod <chr> "NM", "NM", "NM", "NM", "NM", "NM", "NM", "NM~
## $ tohuc <chr> "070200090403", "070200030705", "070200030603~
## $ noncontributingareaacres <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ~
## $ noncontributingareasqkm <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ~
## $ globalid <chr> "{991AACA7-E29C-11E2-8094-0021280458E6}", "{A~
## $ shape_Length <dbl> 0.6371651, 0.4568287, 0.7015241, 0.5549508, 0~
## $ shape_Area <dbl> 0.005950712, 0.004850696, 0.007178123, 0.0119~
## $ shape <MULTIPOLYGON [m]> MULTIPOLYGON (((159993 2287..., ~For example, California
Add unique identifier
# to test (efficiently), here I filter to just one state
# then I create a unique row identifier - this is optional
# otherwise can just use as.numeric(huc12) - but numbers are large and messy
ca <- watersheds %>%
filter(states == "CA") %>%
dplyr::select(huc12) %>%
mutate("id" = row_number())
# save id to huc12 'codebook' for safekeeping
codebook <- ca %>%
st_drop_geometry()
glimpse(codebook)## Rows: 4,222
## Columns: 2
## $ huc12 <chr> "180702010002", "180702020101", "180300122304", "180300122401", ~
## $ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1~write_csv(codebook,
"huc12_id_codebook.csv")
ggplot() +
geom_sf(data = ca,
aes(fill = id),
lwd = 0,
show.legend = FALSE) +
theme_void()
Convert to tif
ca_tif <- ca %>%
fasterize::fasterize(
raster = template,
field = "id"
)
# view & plot
ca_tif## class : RasterLayer
## dimensions : 11747, 18563, 218059561 (nrow, ncol, ncell)
## resolution : 250, 250 (x, y)
## extent : -2357953, 2282797, 238542.6, 3175293 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 4222 (min, max)rasterVis::levelplot(
crop(ca_tif, extent(ca)),
margin = FALSE
)
Save
# save raster
writeRaster(
ca_tif,
"test_ca_raster.tif",
overwrite = TRUE
)