Ratified Raster Brick

This is a script that takes selected fields (selected_fields) from a shape file and turns each into a raster. If the field is character-based categories, it ratifies the raster. At the end, you have a list of rasters. You can then bind these into a brick or perform other list-wise operations on them.

# The input file geodatabase
fgdb <- "land.gdb"

# Read the feature class
fc <- st_read(dsn=fgdb,layer="all_land")
forest <- read.csv("IDs_to_keep.csv")
world.map <- fc[fc$gridcode %in% forest$gridcode,]

ext <- extent(world.map)
gridsize <- 30
r <- raster(ext, res=gridsize)

selected_fields <- c("Own", "Cover", "Secondary")

test <- lapply(selected_fields, function(x) {
   myraster <- fasterize(world.map, r, field=x)
   names(myraster) <- x
   if(is.factor(world.map[,x][[1]])) {
      field <- x
      myraster <- ratify(myraster)
      rat <- levels(myraster)[[1]]
      rat[[field]] <- levels(world.map[,x][[1]])[rat$ID]
   levels(myraster) <- rat}

Weekend GIS Fun: Part 3

I had been referencing stops from the original north shore geology field trip guide, until I found an updated version of the guide in a book (embedded below p.132-144)!

In this chapter, he lists the military grid reference system (MGRS) coordinates (Green 2011) so I entered them into a spreadsheet. I then downloaded it as a comma separated values (CSV) sheet to convert it to a spatial object. The main highlight of this script is how I dealt with different MGRS zones for points in the same data frame:

plot_locations_HARV <- read.csv("geology.csv")
co <- st_read(dsn="counties_in_minnesota", layer="mn_county_boundaries")
pts <- lapply(1:nrow(plot_locations_HARV),function(x){
   point <- plot_locations_HARV[x,]
   epsg <- ifelse(as.integer(gsub("[^[:digit:]]","",point$utm))==15,2027,2028)
   pt <- st_as_sf(point, coords = c("E", "N"))
   this_crs <- st_crs(paste0("+init=epsg:",epsg))
   st_crs(pt) <- this_crs
   pt_project <- st_transform(pt, st_crs(co))
df <- do.call(rbind, pts)
st_write(df, dsn="geology", layer="geology", driver= "ESRI Shapefile")

I assigned the numeric part of the MGRS zone as the Universal Transverse Mercator (UTM) zone. From there, I did my best assigning a coordinate reference system (CRS) to the points accordingly. I ultimately wanted it in the same projection as the MN county layer.

I’ll resume my efforts to “ground truth” the locations as soon as the conditions are right, and to add some additional points corresponding to details in the stop descriptions. If next weekend doesn’t allow the ice to thaw, that will probably mean in the spring!

Literature Cited

Green, John C., et al. “The North Shore Volcanic Group: Mesoproterozoic plateau volcanic rocks of the Midcontinent Rift system in northeastern Minnesota.” Archean to Anthropocene: Field Guides to the Geology of the Mid-Continent of North America: Geological Society of America Field Guide 24 (2011): 121-146.

Latest Error

I’m getting this error after a 137

/panfs/roc/msisoft/R/3.6.0/lib64/R/bin/BATCH: line 60: 30473 Killed ${R_HOME}/bin/R -f ${in} ${opts} ${R_BATCH_OPTIONS} > ${out} 2>&1

It ends on this line…

fr <- sapply(1:nrow(forest),function(y) fasterize(forest[y,], cr[[1]]))

Since smaller memory processes run without error, I can only assume that storing a vector of larger cropped rasters is the problem. So, onto a velox-based solution…

Parallel Raster Extraction

This is my job submission script for the raster extraction.

#!/bin/bash -l
#PBS -l walltime=48:00:00,nodes=1:ppn=24,mem=61gb
#PBS -m abe
#PBS -M myaddress@email.com
module load gdal/2.0.1
module load proj
module load R
module load ompi
module load udunits
module load geos

mpirun -v -hostfile $PBS_NODEFILE -np 1 R --slave CMD BATCH ./extract.R /dev/stdout

This is the R script.

all <- stack("all.tif")
forest <- st_read(dsn="./forest",layer="forest_simple")

raster_cells <- extract(all,forest)



Weekend GIS Fun: Part 2

See part 1 for the relevant shape file and description of what’s going on! In the last post, we identified the map codes of interest to subset the shape file. Here’s some code for working with it in R.

bedrock <- readOGR(dsn="~/Arrowhead_South_Bedrock_Geology","Arrowhead_South_Bedrock_Geology")
basalt = subset(bedrock,map_label=="Mnd" | map_label=="Mna")
writeOGR(basalt,dsn=".",layer="basalt",driver="ESRI Shapefile")

So now we have a shape file of the diabase-textured basalt and andesite! To take it in the field, I found this great iPhone app. Once I had the shape file loaded into a project, I could see my GPS location for navigation.

Here’s a non-GIS based hint of where to find pipe vesicles. In biological news, I saw a Franklin’s ground squirrel at Brighton Beach today!


I currently work on several projects…

  • a contracted project for The Nature Conservancy, relating LiDAR to ground-based measures
  • database & website management for the Coastal Wetland Monitoring (CWM) project
  • Lessard-Sams Outdoor Heritage Council
  • various database/web management needs for the bird & wildlife lab