Week 4

Raster Data

Recap: sf and spatial concepts

Spatial phenomena can be thought of discrete objects with clear boundaries or as continuous phenomenon that can be observed everywhere, with no natural boundaries.
We have described these as objects and fields (Kuhn, 2012)
Objects are usually represented by vector data consisting of:
- a geometry (simple features (sfc, sfg))
- some attribute information (data.frame)
In R these are unified as a sf object
Field data is typically represented by raster data
For this, we will begin our discussions using the terra package

`terra`

Like sf, terra is an implementation of standard raster data model
The model is used by all GIS platforms
Represented continuous data either as continuous or categorical values as a regular set of cells in a grid (matrix)
cells have: (1) resolution, (2) infered cell coordinate (centroid) (3) the coordinate and value apply to the entire cell area

Recap: R data structures

Vector:
- A vector can have dimensions
  - A 1D vector in a collection of values
  - A 2D vector is a matrix
  - A 3D vector is an array
List: a collection of objects
data.frame: a list with requirement of equal length column (vectors)
data.frames and lists (sfc) defined our vector model
Arrays will define our raster model

Spatial Extent

One last topic with respect to vector data (that will carry us into raster) is the idea of an extent:

(ny <- AOI::aoi_get(state = "NY") |> 
   st_transform(5070) |> 
   dplyr::select(name))
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1319502 ymin: 2149150 xmax: 1997508 ymax: 2658543
#> Projected CRS: NAD83 / Conus Albers
#>       name                       geometry
#> 1 New York MULTIPOLYGON (((1661335 263...

In geometry, the minimum bounding box for a point set (stored as POINT, POLYLINE, POLYGON) in N dimensions is “…the box with the smallest measure within which all the points lie.”

We can extract bounding box coordinates with `st_bbox`

returns: an object of class bbox of length 4.

(bb = st_bbox(ny))
#>    xmin    ymin    xmax    ymax 
#> 1319502 2149150 1997508 2658543

class(bb)
#> [1] "bbox"

typeof(bb)
#> [1] "double"

There is a method for creating an sfc form a `bbox` object

(bb = st_as_sfc(bb))
#> Geometry set for 1 feature 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1319502 ymin: 2149150 xmax: 1997508 ymax: 2658543
#> Projected CRS: NAD83 / Conus Albers

class(bb)
#> [1] "sfc_POLYGON" "sfc"

typeof(bb)
#> [1] "list"

Result:

plot(bb, border = rgb(0,0,1))
plot(ny, add = TRUE, col = rgb(1,0,0, .5))

Extents can be discritized in a number of ways:

grid = st_make_grid(bb)
plot(ny$geometry)
plot(grid, add = TRUE)

grid1km = st_make_grid(bb, cellsize = 10000)
plot(ny$geometry)
plot(grid1km, add = TRUE)

What makes of a regular tesselation

length(grid1km) # how many grid tiles
#> [1] 3468

. . .

mapview::npts(grid1km) # how many points?
#> [1] 17340

. . .

mapview::npts(grid1km) * 2 # how many X and Y?
#> [1] 34680

. . .

mapview::npts(grid1km) / length(grid) # how many points per tile?
#> [1] 173.4

. . .

sqrt(st_area(grid1km)[1]) # length of each tile?
#> 10000 [m]

st_bbox(grid1km) # extent of grid
#>    xmin    ymin    xmax    ymax 
#> 1319502 2149150 1999502 2659150

Alternative representation

Regular grids can also be indexed by their centroids

cent <- st_centroid(grid1km)

plot(ny$geometry)
plot(cent, add = TRUE, pch = 16, cex = .25)

length(cent) # how many grid tiles
#> [1] 3468

mapview::npts(grid1km) # how many points?
#> [1] 17340

mapview::npts(grid1km) * 2 # how many X and Y?
#> [1] 34680

Equal area from centroid

We can use our voroni diagram the show that the area closest to a cell centroid is the cell itself.

vor = st_union(cent) |> 
  st_voronoi() |> 
  st_cast() |> 
  st_intersection(bb)

plot(ny$geometry); plot(vor, add = TRUE)

Raster Model

The raster model is one of the earliest and most widely used data models within geographic information systems (Tomlin, 1990; Goodchild, 1992, Maguire, 1992).
Typically used to record, analyze and visualize data with a continuous nature such as elevation, temperature (“GIS”), or reflected or emitted electromagnetic radiation (“Remote Sensing”)
Quotes are used because you’ll find from a data perspective these differences are artificial and a product of the ESRI/ENVI/ERDAS divide

. . .

The term raster originated from the German word for screen, implying a series of orthogonality oriented parallel lines.
Digital raster objects most often take the form of a regularly spaced, grid-like pattern of rows and columns
Each element referred to as a cell, pixel, or grid point.

Many terms mean the same thing …

The entire raster is sometimes referred to as an “image”, “array”, “surface”, “matrix”, or “lattice” (Wise, 2000).
The all mean the same thing…
Cells of the raster are most often square, but may be rectangular (with differing resolutions in x and y directions) or other shapes that can be tessellated such as triangles and hexagons (Figure below from Peuquet, 1984).

Photos and Computers …

Aerial Imagery (really just a photo 😄)

What is stored in these cells?

Categorical Values (integer/factor)

Continuous Values (numeric)

Spectral Values

Either Color, or sensor

Any digital image contains an RBG channel for color:

Red, green, and blue are the three additive colors (primary colors of light)
In R, colors can be defined using the RBG channels

(rgb(1,0,0)) # red
#> [1] "#FF0000"
(rgb(0,.54,.96)) # UCSB navy
#> [1] "#008AF5"
(rgb(254,188,17, maxColorValue = 255)) # UCSB navy
#> [1] "#FEBC11"

Pure RBG

par(mfrow = c(1,3), mar = c(0,0,0,0))
plot(ny$geometry,   col = rgb(1,0,0)) # red
plot(ny$geometry,   col = rgb(0,1,0)) # green
plot(ny$geometry,   col = rgb(0,0,1)) # blue

RGB and bytes/bits

The red, green and blue use 8 bits each (1 byte), which each have integer values from 0 to 255.
This makes 256^3 = 16,777,216 possible colors.
See more here
How does this relate to the selection of 256x256 pixel tiles in web maps?

The secondary colors in an RGB color wheel are cyan, magenta, and yellow because these are the three subtractive colors

Think of your printer and the CMYK ink cartridges! Where black, is the absence of color (0,0,0)

par(mfrow = c(1,4), mar = c(0,0,0,0))
plot(ny$geometry, col = rgb(0,1,1))  # cyan
plot(ny$geometry, col = rgb(1,0,1)) # Magenta
plot(ny$geometry, col = rgb(1,1,0)) # Yellow
plot(ny$geometry, col = rgb(0,0,0)) # Key (black)

Why do we care?

Pixels are the base unit of raster data and have a resolution
This is the X and the Y dimension of each cell in the units of the CRS

Resolution drives image clarity (granulairty)

Higher resolution (smaller cells) = more detail, but bigger data!

Raster images seek to discritize the real world into cell-based values

Again either integer (categorical), continuous, or signal

All rasters have an extent!

This is the same extent as a bounding box
Can be described as 4 values (xmin,ymin,xmax,ymax)

Implicity Coordinates

Unlike vector data, the raster data model stores the coordinate of the grid cells indirectly
Coordinates are derived from the reference (Xmin,Ymin) the resolution, and the cell index (e.g. [100,150])
For example: If we want the coordinates of a value in the 3rd row and the 40th column of a raster matrix, we have to move from the origin (Xmin, Ymin) (3 x Xres) in x-direction and (40 x Yres) in y-direction

So, any image (.png, .tif, .gif) can be read as a raster…
The raster is defined by the extent and resolution of the cells
To be spatial, the extent (thus coordinates) must be grounded in an CRS

(img = terra::rast('images/17-raster-extent.png'))
#> class       : SpatRaster 
#> dimensions  : 788, 1067, 4  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 1067, 0, 788  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source      : 17-raster-extent.png 
#> names       : 17-rast~xtent_1, 17-rast~xtent_2, 17-rast~xtent_3, 17-rast~xtent_4

terra::plotRGB(img, r = 1, g = 2, b = 3)

Raster Data in R

A SpatRast represents single-layer (variable) raster data.

A SpatRast always stores the fundamental parameters that describe it. - The number of columns and rows, - The spatial extent - The Coordinate Reference System.

In addition, a SpatRast can store information about the file where raster values are stored (if there is such a file).

Here we construct an empty raster:

(r <- rast(ncol=20, nrow=20, xmax=-80, xmin=-120, ymin=20, ymax=60))
#> class       : SpatRaster 
#> dimensions  : 20, 20, 1  (nrow, ncol, nlyr)
#> resolution  : 2, 2  (x, y)
#> extent      : -120, -80, 20, 60  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)

Extract Diminisionality …

nrow(r)
#> [1] 20
ncol(r)
#> [1] 20
ncell(r)
#> [1] 400
nlyr(r)
#> [1] 1

Raster Values

values can be extracted with values()

head(values(r))
#>      lyr.1
#> [1,]   NaN
#> [2,]   NaN
#> [3,]   NaN
#> [4,]   NaN
#> [5,]   NaN
#> [6,]   NaN

Assigning a values as a vector

values(r) <- 1:ncell(r)

r
#> class       : SpatRaster 
#> dimensions  : 20, 20, 1  (nrow, ncol, nlyr)
#> resolution  : 2, 2  (x, y)
#> extent      : -120, -80, 20, 60  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source(s)   : memory
#> name        : lyr.1 
#> min value   :     1 
#> max value   :   400

head(values(r))
#>      lyr.1
#> [1,]     1
#> [2,]     2
#> [3,]     3
#> [4,]     4
#> [5,]     5
#> [6,]     6

plot(r)

Raster is S4

Compared to S3, the S4 object system is much stricter, and much closer to other OO systems.
What does this mean for us?
Data is structured with a “representation”
which is a list of slot (or attributes), giving their names and classes, accessed with @

str(r, max.level = 2)
#> S4 class 'SpatRaster' [package "terra"]
r@pntr
#> C++ object <0x11d191f10> of class 'SpatRaster' <0x10d012710>
r@pntr$ncol
#> Class method definition for method ncol()
#> function () 
#> {
#>     " unsigned long ncol()  \n   docstring : ncol"
#>     .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x1030ca1b0>, 
#>         dll = list(name = "Rcpp", path = "/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/libs/Rcpp.so", 
#>             dynamicLookup = TRUE, handle = <pointer: 0x7aab93b0>, 
#>             info = <pointer: 0x1030ba810>, forceSymbols = FALSE), 
#>         numParameters = -1L), <pointer: 0x10d012710>, <pointer: 0x10310f510>, 
#>         .pointer)
#> }
#> <environment: 0x10e2f3db0>

Access via Function

ext(r)
#> SpatExtent : -120, -80, 20, 60 (xmin, xmax, ymin, ymax)
crs(r)
#> [1] "GEOGCRS[\"WGS 84 (CRS84)\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"World\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"OGC\",\"CRS84\"]]"
nrow(r)
#> [1] 20
ncol(r)
#> [1] 20
head(values(r))
#>      lyr.1
#> [1,]     1
#> [2,]     2
#> [3,]     3
#> [4,]     4
#> [5,]     5
#> [6,]     6

Multi-layers

In many cases multi-variable raster data sets are used. Variables can related to time or measurements

A SpatRaster is a collection ofobjects with the same spatial extent and resolution.

In essence it is a list of objects, or a 3D array.

Remember our Array data?

(v = 1:27)
#>  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
#> [26] 26 27
(arr = array(v, dim = c(3,3,3)))
#> , , 1
#> 
#>      [,1] [,2] [,3]
#> [1,]    1    4    7
#> [2,]    2    5    8
#> [3,]    3    6    9
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3]
#> [1,]   10   13   16
#> [2,]   11   14   17
#> [3,]   12   15   18
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3]
#> [1,]   19   22   25
#> [2,]   20   23   26
#> [3,]   21   24   27

An array is a single file so it can be “rasterized”

b = rast(arr)
plot(b)

Core Raster Principles

1. Bounding Box / Extents

Geometries have extents that define the maximum and minimum coverage of the shape in a coordinate reference system

Image Source: National Ecological Observatory Network (NEON)

2. Extent

When dealing with objects, the extent (or bbox) is derived from the coordinate set
When dealing with raster data, the extent is a fondational component of the raster data structure
- That is, we need to know the area the raster is covering!

3. Discretization

Once we know the extent, we need to know how that space is split up

Two complimentary bit of information can tell us this:

Resolution (res)
Number of row and number of columns (nrow/ncol)

So,

A raster is made of an extent, and a resolution / row-column structure

A vector of values fill that structure (same way a vector in R can have diminisons)
- These values are often scaled to integers to reduce file size
Values are referenced in cartisian space, based on cell index
A CRS along with the extent, can provide spatial reference / coordinates

General Process

Almost all remote sensing / image analysis begins with the same basic steps:

Identifying an area of interest (AOI)
Identifying and downloading the relevant images or products
Analyzing the raster products

The definition of a AOI is critical because raster data in continuous, therefore we need to define the bounds of the study rather then the bounds of the objects

But, objects often (even typically) define our bounds

Find elevation data for Fort Collins:

Define the AOI

bb = read_csv("../labs/data/uscities.csv") |>
  st_as_sf(coords = c("lng", "lat"), crs = 4326) |> 
  filter(city == "Fort Collins") |> 
  st_transform(5070) |> 
  st_buffer(50000) |> 
  st_bbox() |> 
  st_as_sfc() |> 
  st_as_sf()

Read data from elevation map tiles, for a specific zoom, and crop to the AOI

elev = elevatr::get_elev_raster(bb, z = 11) |> crop(bb)
writeRaster(elev, filename = "data/foco-elev.tif", overwrite = TRUE)

The resulting raster …

(elev = rast("data/foco-elev.tif"))
#> class       : SpatRaster 
#> dimensions  : 3450, 3450, 1  (nrow, ncol, nlyr)
#> resolution  : 28.98965, 28.98965  (x, y)
#> extent      : -810156.6, -710142.3, 1934613, 2034627  (xmin, xmax, ymin, ymax)
#> coord. ref. : +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      : foco-elev.tif 
#> name        : foco-elev 
#> min value   :      1382 
#> max value   :      4346

Raster Values Continuity

v = values(elev)
class(v)
#> [1] "matrix" "array"
length(v)
#> [1] 11902500

elev
#> class       : SpatRaster 
#> dimensions  : 3450, 3450, 1  (nrow, ncol, nlyr)
#> resolution  : 28.98965, 28.98965  (x, y)
#> extent      : -810156.6, -710142.3, 1934613, 2034627  (xmin, xmax, ymin, ymax)
#> coord. ref. : +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      : foco-elev.tif 
#> name        : foco-elev 
#> min value   :      1382 
#> max value   :      4346

Raster Values

# The length of the vector is equal to the rows * columns
length(v) == nrow(elev) * ncol(elev)
#> [1] TRUE
# The span of the x extent divided by the resolution equals the raster rows
((xmax(elev) - xmin(elev)) / res(elev)[1]) == ncol(elev) 
#> [1] TRUE
# The span of the x extent divided by the number of rows equals the raster resolution
((xmax(elev) - xmin(elev)) / ncol(elev)) == res(elev)[1] 
#> [1] TRUE

All image files are the same!

download.file(url = "https://a.tile.openstreetmap.org/18/43803/104352.png",
              destfile = "data/104352.png")
img = png::readPNG("data/104352.png")
class(img)
#> [1] "array"
typeof(img)
#> [1] "double"
dim(img)
#> [1] 256 256   3

img[1,1,1:3]
#> [1] 1.0000000 1.0000000 0.9333333
rgb(1.0000000, 0.9607843, 0.8980392)
#> [1] "#FFF5E5"

Raster Algebra

So our raster data is stored as a large numeric array/vector
Many generic functions allow for simple algebra on Raster objects,
These include:
- normal algebraic operators such as +, -, *, /
- logical operators such as >, >=, <,==, !
- functions like abs, round, ceiling, floor, trunc, sqrt, log, log10, exp, cos, sin, atan, tan, max, min, range, prod, sum, any, all

Raster Algebra

In these functions you can mix SpatRast objects with numbers, as long as the first argument is a raster object.
That means you can add 100 to a raster object but not a raster object to 100

# GOOD
raster + 100

# BAD
100 + raster

For example:

elev + 100
#> class       : SpatRaster 
#> dimensions  : 3450, 3450, 1  (nrow, ncol, nlyr)
#> resolution  : 28.98965, 28.98965  (x, y)
#> extent      : -810156.6, -710142.3, 1934613, 2034627  (xmin, xmax, ymin, ymax)
#> coord. ref. : +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(s)   : memory
#> varname     : foco-elev 
#> name        : foco-elev 
#> min value   :      1482 
#> max value   :      4446
log10(elev)
#> class       : SpatRaster 
#> dimensions  : 3450, 3450, 1  (nrow, ncol, nlyr)
#> resolution  : 28.98965, 28.98965  (x, y)
#> extent      : -810156.6, -710142.3, 1934613, 2034627  (xmin, xmax, ymin, ymax)
#> coord. ref. : +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(s)   : memory
#> varname     : foco-elev 
#> name        : foco-elev 
#> min value   :  3.140508 
#> max value   :  3.638090

Replacement

Raster values can be replaced on a conditional statements
Doing this changes the underlying data!
If you want to retain the original data, you must make a copy of the base layer

plot(elev)

elev2 = elev #<<
elev2[elev2 <= 1500] = NA #<<
plot(elev2)

Modifying a raster

When we want to modify the extent of a raster we can clip it to a new bounds

crop: lets you reduce the extent of a raster to the extent of another, overlapping object:

#remotes::install_github("mikejohnson51/AOI")
fc = AOI::geocode("Fort Collins", bbox = TRUE) |> 
  st_transform(crs(elev))

plot(elev)
plot(fc, add = TRUE, col = NA)

fc_elev = crop(elev, fc) #<<
plot(fc_elev)

Modifying the underlying data:

mask: mask takes an input object (sf, sp, or raster) and set anything not undelying the input to a new value (default = NA)

library(osmdata)

osm = osmdata::opq(st_bbox(st_transform(fc,4326))) |> 
  add_osm_feature("water") |> 
  osmdata_sf()

(poly = osm$osm_polygons |> 
  st_transform(crs(elev)))
#> Simple feature collection with 181 features and 17 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -769181 ymin: 1977155 xmax: -753314.8 ymax: 2000444
#> Projected CRS: PROJCRS["unknown",
#>     BASEGEOGCRS["NAD83",
#>         DATUM["North American Datum 1983",
#>             ELLIPSOID["GRS 1980",6378137,298.257222101004,
#>                 LENGTHUNIT["metre",1]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         ID["EPSG",4269]],
#>     CONVERSION["Albers Equal Area",
#>         METHOD["Albers Equal Area",
#>             ID["EPSG",9822]],
#>         PARAMETER["Latitude of false origin",23,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8821]],
#>         PARAMETER["Longitude of false origin",-96,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8822]],
#>         PARAMETER["Latitude of 1st standard parallel",29.5,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8823]],
#>         PARAMETER["Latitude of 2nd standard parallel",45.5,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8824]],
#>         PARAMETER["Easting at false origin",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8826]],
#>         PARAMETER["Northing at false origin",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8827]]],
#>     CS[Cartesian,2],
#>         AXIS["easting",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]],
#>         AXIS["northing",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]]
#> First 10 features:
#>            osm_id              name         alt_name area basin boat
#> 6051656   6051656              <NA>             <NA> <NA>  <NA> <NA>
#> 23598705 23598705          Lee Lake             <NA> <NA>  <NA> <NA>
#> 23598868 23598868              <NA>             <NA> <NA>  <NA> <NA>
#> 23598945 23598945      College Lake             <NA> <NA>  <NA> <NA>
#> 23599224 23599224       Warren Lake             <NA> <NA>  <NA> <NA>
#> 23599227 23599227 Harmony Reservoir             <NA> <NA>  <NA> <NA>
#> 23599320 23599320              <NA>             <NA> <NA>  <NA> <NA>
#> 23599327 23599327              <NA>             <NA> <NA>  <NA> <NA>
#> 44074911 44074911     Lake Sherwood Nelson Reservoir <NA>  <NA> <NA>
#> 47471957 47471957              <NA>             <NA> <NA>  <NA> <NA>
#>          description  ele gnis:feature_id intermittent landuse layer natural
#> 6051656         <NA> <NA>            <NA>         <NA>    <NA>  <NA>    <NA>
#> 23598705        <NA> <NA>            <NA>         <NA>    <NA>  <NA>   water
#> 23598868        <NA> <NA>            <NA>         <NA>    <NA>  <NA>    <NA>
#> 23598945        <NA> <NA>            <NA>         <NA>    <NA>  <NA>   water
#> 23599224        <NA> <NA>            <NA>         <NA>    <NA>  <NA>   water
#> 23599227        <NA> <NA>            <NA>         <NA>    <NA>  <NA>   water
#> 23599320        <NA> <NA>            <NA>         <NA>    <NA>  <NA>   water
#> 23599327        <NA> <NA>            <NA>         <NA>    <NA>  <NA>   water
#> 44074911        <NA> 1512          201223         <NA>    <NA>  <NA>   water
#> 47471957        <NA> <NA>            <NA>         <NA>    <NA>  <NA>   water
#>          place salt source     water                       geometry
#> 6051656   <NA> <NA>   <NA>      <NA> POLYGON ((-762901.7 1988981...
#> 23598705  <NA> <NA>   <NA>      lake POLYGON ((-765209.9 1991246...
#> 23598868  <NA> <NA>   <NA>      <NA> POLYGON ((-767955.5 1985345...
#> 23598945  <NA> <NA>   <NA>      lake POLYGON ((-766871.6 1989267...
#> 23599224  <NA> <NA>   <NA> reservoir POLYGON ((-760465.5 1983486...
#> 23599227  <NA> <NA>   <NA> reservoir POLYGON ((-760046.6 1982348...
#> 23599320  <NA> <NA>   <NA>      pond POLYGON ((-762458.2 1991222...
#> 23599327  <NA> <NA>   <NA>      pond POLYGON ((-763462.8 1992438...
#> 44074911  <NA> <NA>   <NA> reservoir POLYGON ((-758628.3 1984729...
#> 47471957  <NA> <NA>   <NA> reservoir POLYGON ((-764440.4 1992962...

plot(fc_elev)
plot(poly, col = "blue", add  = TRUE)

ma =  mask(fc_elev, poly)
ma2 = mask(fc_elev, poly, inverse = TRUE)
plot(c(fc_elev, ma, ma2))

What is `mask` doing?

NA * 7
#> [1] NA

. . .

mask_r = rasterize(poly, fc_elev, background = NA)
plot(mask_r)

base_mask = mask_r * fc_elev
plot(base_mask)

Crop or/and mask

Crop is more efficient then mask
Often you will want to mask and crop a raster
The correct way to do this is crop then mask

cm = crop(fc_elev, poly) |>  
  mask(poly)

plot(cm)

Aggregate and disaggregate

aggregate and disaggregate allow for changing the resolution of a Raster object.
This is similar to the zoom scaling on a web map except the scale factor is not set to 2
For aggregate, you need to specify a function determining what to do with the grouped cell values (default = mean).

plot(fc_elev)
plot(rpoly, add = T)

agg = aggregate(fc_elev, 10, fun = max) #<<
plot(agg)
plot(rpoly, add = T)

`app`

Just like a vector, we can apply functions over a raster with app

These types of formulas are very useful for thresholding analysis

Question: separate Fort Collins into the higher and lower elevations

FUN = function(x){ ifelse(x < mean(x), 1, 2) }

elev3 = app(elev, FUN) #<<
plot(elev3, col = c("red", "blue"))

Read in the saved raster file

(r = rast("data/foco-elev.tif"))
#> class       : SpatRaster 
#> dimensions  : 3450, 3450, 1  (nrow, ncol, nlyr)
#> resolution  : 28.98965, 28.98965  (x, y)
#> extent      : -810156.6, -710142.3, 1934613, 2034627  (xmin, xmax, ymin, ymax)
#> coord. ref. : +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      : foco-elev.tif 
#> name        : foco-elev 
#> min value   :      1382 
#> max value   :      4346

Create a conditional (threshold) mask

threshold = function(x) {ifelse(x <= 1520 , NA, 1)}

threshold(1600)
#> [1] 1
threshold(-100)
#> [1] NA

(m = app(r, threshold))
#> class       : SpatRaster 
#> dimensions  : 3450, 3450, 1  (nrow, ncol, nlyr)
#> resolution  : 28.98965, 28.98965  (x, y)
#> extent      : -810156.6, -710142.3, 1934613, 2034627  (xmin, xmax, ymin, ymax)
#> coord. ref. : +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(s)   : memory
#> name        : lyr.1 
#> min value   :     1 
#> max value   :     1

Results

plot(r)

plot(m)

Multiply cell-wise

algebraic, logical, and functional operations act on a raster cell-wise

elev_cut = m * r
plot(elev_cut, col = viridis::viridis(256))

Reclassify

Reclassify is a function that allows you to change the values of a raster based on a set of rules
The rules are defined in a data frame with three columns:
- min = the minimum value of the range
- max = the maximum value of the range
- lab = the new value to assign to that range

(rcl = data.frame(min = seq(1500,1590,10), max =  seq(1510,1600,10), lab = c(0:9)))
#>     min  max lab
#> 1  1500 1510   0
#> 2  1510 1520   1
#> 3  1520 1530   2
#> 4  1530 1540   3
#> 5  1540 1550   4
#> 6  1550 1560   5
#> 7  1560 1570   6
#> 8  1570 1580   7
#> 9  1580 1590   8
#> 10 1590 1600   9

(rc = classify(elev_cut, rcl, include.lowest = TRUE))
#> class       : SpatRaster 
#> dimensions  : 3450, 3450, 1  (nrow, ncol, nlyr)
#> resolution  : 28.98965, 28.98965  (x, y)
#> extent      : -810156.6, -710142.3, 1934613, 2034627  (xmin, xmax, ymin, ymax)
#> coord. ref. : +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(s)   : memory
#> name        : lyr.1 
#> min value   :     2 
#> max value   :  4346

(s = c(r, m, elev_cut, rc) |> 
  setNames(c("elevation", "elev-mask", "terrain", "topography")))
#> class       : SpatRaster 
#> dimensions  : 3450, 3450, 4  (nrow, ncol, nlyr)
#> resolution  : 28.98965, 28.98965  (x, y)
#> extent      : -810156.6, -710142.3, 1934613, 2034627  (xmin, xmax, ymin, ymax)
#> coord. ref. : +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 
#> sources     : foco-elev.tif  
#>               memory  
#>               memory  
#>               memory  
#> names       : elevation, elev-mask, terrain, topography 
#> min values  :      1382,         1,    1521,          2 
#> max values  :      4346,         1,    4346,       4346

plot(s, col = viridis::viridis(256))

Real Example: Classify Rainfall Regions of Colorado

#remotes::install_github("mikejohnson51/climateR")
library(climateR)

AOI = AOI::aoi_get(state = 'CO') 

system.time({ prcp = climateR::getTerraClim(AOI, "ppt", 
                                            startDate = "2000-01-01", endDate = '2005-12-31') })
#>    user  system elapsed 
#>   0.255   0.029   2.802

# More on global below ...
quarts = global(prcp$ppt, fivenum)

(quarts = colMeans(quarts))
#>         X1         X2         X3         X4         X5 
#>   6.271233  18.434247  27.995890  41.608219 105.642466

(rcl = data.frame(quarts[1:4], quarts[2:5], 1:4))
#>    quarts.1.4. quarts.2.5. X1.4
#> X1    6.271233    18.43425    1
#> X2   18.434247    27.99589    2
#> X3   27.995890    41.60822    3
#> X4   41.608219   105.64247    4

terra::classify(mean(prcp$ppt), rcl, include.lowest=TRUE) |> 
   plot(col = blues9)

Object/Field Interaction

For this example, we used OSM to extact the river data for the area of interest. We will talk more about OSM next week:

foco_rivers <- read_sf("data/foco-rivers-osm.gpkg")

Lets find the longest river segment IN our extent

river = foco_rivers |> 
  st_transform(crs(r)) |> 
  st_intersection(st_as_sfc(st_bbox(r))) %>% 
  mutate(length = st_length(.)) |> 
  slice_max(length, n = 1)

plot(r)
plot(river, add = TRUE, col = "blue", lwd = 2)

Value Extraction

Often, we want to know the profile and sinousity of a river
To do this, we need to know the inlet and outlet as well as the straight line connector

inlet  = head(st_cast(river, "POINT"), 1)
outlet = tail(st_cast(river, "POINT"), 1) 
pts    = bind_rows(inlet, outlet) 

line = st_cast(st_union(pts), "LINESTRING")

plot(r)
plot(river, add = TRUE, col = "blue", lwd = 2)
plot(line, add = TRUE, col = "black", lwd = 2)
plot(outlet$geom, add = TRUE, pch = 16, col = "red")
plot(inlet$geom,  add = TRUE, pch = 16, col = "green")

Sinuosity

Channel sinuosity is calculated by dividing the length of the stream channel by the straight line distance between the end points of the selected channel reach.

(sin = st_length(river) / st_length(line))
#> 1.400515 [1]

River Slope:

The change in elevation between the inlet/outlet divided by the length (rise/run) give us the slope of the river:

To calculate this, we must extract elevation values at the inlet and outlet:

(elev = extract(r, pts))
#>   ID foco-elev
#> 1  1      1813
#> 2  2      1675

100 * (elev$`foco-elev`[1] - elev$`foco-elev`[2]) / units::drop_units(st_length(river))
#> [1] 0.600446

River profile

What does the elevation profile of the river look like?

profile = extract(r, river)$`foco-elev`

plot(profile, type = "l")
lines(zoo::rollmean(profile,k = 10), 
      col = "darkred", lwd = 3)

Map Algebra

Dana Tomlin (Tomlin 1990) defined a framework for the analyizing field data stored as grided values.
He called this framework map algebra.
Map algebra operations and functions are broken down into four types:
- local
- focal
- zonal
- global

Local

Local operations and functions are applied to each individual cell and only involve those cells sharing the same location.
More than one raster can be involved in a local operation.
For example, rasters can be summed ( each overlapping pixels is added)
Local operations also include reclassification of values.

s = c(mean(prcp$ppt), app(prcp$ppt, sd), min(prcp$ppt), max(prcp$ppt)) |> 
  setNames(c("Mean", "StDev", "Min", "Max"))

rasterVis::levelplot(s)

Focal

Also referred to as “neighborhood” operations.
Assigns summary values to the output cells based on the neighboring cells in the input raster.
For example, a cell output value can be the average of 9 neighboring input cells (including the center cell) - this acts as a smoothing function.

Focal

Focal operations require a window (also known as a kernel) to work over
Additionally a kernel also defines the weight each neighboring cell contributes to the summary statistic.
For example, all cells in a 3x3 neighbor could each contribute 1/9th of their value to the summarized value (i.e. equal weight).
The weight can take on a more complex form defined by a function; such weights are defined by a kernel function.
One popular function is a Gaussian weighted function which assigns greater weight to nearby cells than those further away (Toblers first law)

Example: Focal

Lets apply a smoothing kernel to our Fort Collins elevation data over an 25x25 window, using the mean operator

foco = AOI::geocode("Fort Collins", bbox = TRUE) |> st_transform(crs(r))

foco_elev = crop(r, foco)
f1 <- focal(foco_elev, w= matrix(1,nrow=25,ncol=25), fun=mean)

Results

plot(foco_elev)

plot(f1)

What did we do?

matrix(1,nrow=25,ncol=25)
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
#>  [1,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>  [2,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>  [3,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>  [4,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>  [5,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>  [6,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>  [7,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>  [8,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>  [9,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [10,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [11,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [12,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [13,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [14,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [15,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [16,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [17,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [18,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [19,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [20,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [21,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [22,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [23,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [24,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#> [25,]    1    1    1    1    1    1    1    1    1     1     1     1     1
#>       [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
#>  [1,]     1     1     1     1     1     1     1     1     1     1     1     1
#>  [2,]     1     1     1     1     1     1     1     1     1     1     1     1
#>  [3,]     1     1     1     1     1     1     1     1     1     1     1     1
#>  [4,]     1     1     1     1     1     1     1     1     1     1     1     1
#>  [5,]     1     1     1     1     1     1     1     1     1     1     1     1
#>  [6,]     1     1     1     1     1     1     1     1     1     1     1     1
#>  [7,]     1     1     1     1     1     1     1     1     1     1     1     1
#>  [8,]     1     1     1     1     1     1     1     1     1     1     1     1
#>  [9,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [10,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [11,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [12,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [13,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [14,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [15,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [16,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [17,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [18,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [19,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [20,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [21,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [22,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [23,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [24,]     1     1     1     1     1     1     1     1     1     1     1     1
#> [25,]     1     1     1     1     1     1     1     1     1     1     1     1

mean(foco_elev[1:25, 1:25][,1])
#> [1] 1610.022

na.omit(values(f1))[1]
#> [1] 1610.022

Zonal

Zonal operations compute a summary values (such as the mean) from cells aggregated to some zonal unit.
Like focal operations, a zone and a mediating function must be defined
The most basis example of a zonal function is aggregation!

aggregate(foco_elev, 10) |> plot()

Zonal Statisics (More advanced)

For more complicated object zones, exactextractr is a fast and effiecient R utility that binds the C++ exactextract tool.
What is the county level mean January rainfall in California?

AOI = AOI::aoi_get(state = "CO", county = "all")
AOI$janPTT = exactextractr::exact_extract(prcp$ppt$`ppt_2000-01-01_total`, AOI, "mean", progress = FALSE)
plot(AOI['janPTT'])

What about the US?

counties <-  AOI::aoi_get(state = "conus", county = "all")
  
jan <- climateR::getTerraClim(counties, "ppt", startDate = "2000-01-01") 
  
counties$janPTT <-  exactextractr::exact_extract(jan$ppt, counties, "mean", progress = FALSE)

plot(counties['janPTT'], border = NA, key.pos = 4)

Global

Global operations make use of some or all input cells when computing an output cell value.
They are a special case of zonal operations with the entire raster represents a single zone.
Examples include generating descriptive statistics for the entire raster dataset

Summary Values

global: computes statistics for the values of each layer in a Raster* object.

elev <- rast('data/foco-elev.tif')
global(elev, mean)
#>               mean
#> foco-elev 1887.874
mean(values(elev), na.rm = TRUE)
#> [1] 1887.874

Why not just `mean()`

In the terra package, functions like max, min, and mean, return a new SpatRast* object (with a value computed for each cell).

In contrast, global returns a single value, computed from the all the values of a layer.

s = c(elev, elev^2, elev*.5)
mean(s) |> plot()

global(s, mean)
#>                     mean
#> foco-elev      1887.8738
#> foco-elev.1 3802035.3412
#> foco-elev.2     943.9369

Mean Monthly Rainfall for Colorado

global

plot(global(prcp$ppt, max)$max, type = "l", 
     ylab = "rainfall", xlab = "month since 2000-01")
lines(global(prcp$ppt, min)$min, type = "l", col = "blue")
lines(global(prcp$ppt, mean)$mean, type = "l", col = "darkred", lwd = 2)

mean()

plot(mean(prcp$ppt), col = blues9)
plot(AOI, add =TRUE, col = NA, lwd = 2)

Kmeans over to Raster Data!

each layer of a SpatRaster is a layer
each layer is a vector of values

library(climateR)
params <-  c("ppt", "tmax", "tmin", "srad", "q")
AOI   <-  AOI::aoi_get(state = "CO")

co <-  climateR::getTerraClim(AOI, params, startDate = "2018-10-01")  %>% 
  unlist() %>% 
  rast() |> 
  setNames(params)

Colorado October 2018 climate

plot(co)

Raster layers are vectors!

values = values(co)
head(values)
#>       ppt tmax  tmin srad    q
#> [1,] 37.8  1.9 125.1 12.1 -2.4
#> [2,] 39.2  2.0 124.2 10.8 -2.8
#> [3,] 39.2  2.0 124.3 10.6 -2.9
#> [4,] 37.7  1.9 124.6 11.6 -2.4
#> [5,] 36.3  1.8 125.0 12.2 -2.2
#> [6,] 34.8  1.7 125.0 13.0 -2.1

Data Prep

Identify NA indices for latter reference
Remove NA values
Scale

idx <- which(!apply(is.na(values), 1, any))
v   <- na.omit(values)
vs  <- scale(v)

(E <- kmeans(vs, 5, iter.max = 100))
#> K-means clustering with 5 clusters of sizes 1471, 3619, 542, 6521, 4677
#> 
#> Cluster means:
#>          ppt       tmax       tmin       srad          q
#> 1  1.7381744  0.2854842 -0.1759393 -0.3587139 -0.0132010
#> 2 -0.2805363 -0.1935889  1.4262188  0.5651088  0.5636098
#> 3  3.2106644  4.8126020 -0.5710265 -1.7692979 -1.2335931
#> 4 -0.5799378 -0.2419152 -0.1384006  0.6815967  0.6486097
#> 5  0.1069063 -0.1604128 -0.7891111 -1.0697441 -1.1933422
#> 
#> Clustering vector:
#>     [1] 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>    [37] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>    [73] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [109] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [145] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 4 4 4 4
#>   [181] 4 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5
#>   [217] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>   [253] 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [289] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [325] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [361] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>   [397] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4
#>   [433] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [469] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [505] 4 4 4 4 4 4 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [541] 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>   [577] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [613] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [649] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [685] 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5
#>   [721] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>   [757] 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [793] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [829] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 4 4 4 4 4 4
#>   [865] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>   [901] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>   [937] 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>   [973] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1009] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1045] 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [1081] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4
#>  [1117] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1153] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1189] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5
#>  [1225] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [1261] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1297] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1333] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 4 4 4 4 4 4 4
#>  [1369] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [1405] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [1441] 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1477] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1513] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1549] 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [1585] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4
#>  [1621] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1657] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1693] 4 4 4 4 4 4 4 4 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5
#>  [1729] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [1765] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4
#>  [1801] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1837] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1873] 4 4 4 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5
#>  [1909] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [1945] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [1981] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2017] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2053] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2089] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2125] 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2161] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2197] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2233] 4 5 5 5 5 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2269] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4
#>  [2305] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2341] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2377] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5
#>  [2413] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2449] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2485] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2521] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 4 4
#>  [2557] 4 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2593] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2629] 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2665] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2701] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4
#>  [2737] 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2773] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2809] 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2845] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [2881] 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5
#>  [2917] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [2953] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4
#>  [2989] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3025] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3061] 4 4 4 4 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3097] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3133] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3169] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3205] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3241] 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3277] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3313] 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3349] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3385] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3421] 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3457] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4
#>  [3493] 4 4 2 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3529] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3565] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5
#>  [3601] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3637] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 2 2 2 2 2 2 2 2 4
#>  [3673] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3709] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3745] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3781] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3817] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4
#>  [3853] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3889] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [3925] 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3961] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [3997] 5 5 5 4 4 4 4 4 4 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4033] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4069] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5
#>  [4105] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4141] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 2 2
#>  [4177] 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4213] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4249] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4285] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4321] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 2 2 2 2 2 2 2 2 2 4 4 4 4
#>  [4357] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4393] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4429] 4 4 4 5 5 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4465] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4501] 5 5 5 5 5 5 5 5 5 5 4 4 4 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4537] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4573] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 4 4 4 4 4
#>  [4609] 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4645] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4681] 4 4 4 4 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4717] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4753] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5
#>  [4789] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4825] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 2 2 2 2 2 2
#>  [4861] 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4897] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [4933] 4 4 4 4 4 4 4 4 4 4 5 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [4969] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5005] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 2 2 2 2 2 2 2 2 2 2 4 4 2 4 4 4 4
#>  [5041] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5077] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5113] 5 5 4 4 5 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5149] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5185] 5 5 5 5 5 5 5 4 4 2 2 2 2 2 2 2 2 2 2 2 4 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5221] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5257] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 4 5 5 5 5 5
#>  [5293] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5329] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4
#>  [5365] 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5401] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5437] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5473] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5509] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 2 2 2 2 2 2 2 2 2 2
#>  [5545] 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5581] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5617] 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5653] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5689] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4
#>  [5725] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5761] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5
#>  [5797] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5833] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5869] 5 5 5 5 4 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5905] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [5941] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [5977] 5 5 5 5 4 4 4 5 5 5 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [6013] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 2 2 2 2
#>  [6049] 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6085] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6121] 4 1 4 4 4 4 4 4 4 4 5 5 5 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4
#>  [6157] 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [6193] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [6229] 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6265] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6301] 4 5 4 4 4 4 4 5 4 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5
#>  [6337] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [6373] 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4
#>  [6409] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6445] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 4 4 4 4 4 4 5
#>  [6481] 4 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [6517] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [6553] 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4
#>  [6589] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6625] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 4 4 4 4 4 4 5 4 4 4 4 5 5 5 5
#>  [6661] 5 5 5 1 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [6697] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2
#>  [6733] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6769] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6805] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [6841] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [6877] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [6913] 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6949] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [6985] 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7021] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7057] 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 2 2 4 4
#>  [7093] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7129] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7165] 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7201] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 5
#>  [7237] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 2 4 4 4 4 4 4 4 4 4 4 4
#>  [7273] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7309] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 5 5 5 5 5
#>  [7345] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7381] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 5 5 2 2 2 2 2 2 2 2 2
#>  [7417] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7453] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7489] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7525] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7561] 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [7597] 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7633] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7669] 4 4 4 4 4 4 4 4 1 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7705] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7741] 2 2 2 2 2 2 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4
#>  [7777] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7813] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 5 5
#>  [7849] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [7885] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 5 5 2 2 5 5 5 2
#>  [7921] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 2 4 4 4 4 4 4 4 4 4
#>  [7957] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [7993] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8029] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8065] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 5 5 5 5 2 2 2 2 2 2 2 2 2 2
#>  [8101] 2 2 2 2 2 2 2 2 2 2 2 2 2 4 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [8137] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [8173] 4 4 4 4 4 4 4 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8209] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8245] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [8281] 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [8317] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5
#>  [8353] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8389] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8425] 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4
#>  [8461] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [8497] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 5 5 5 5 5 5 5 5 5 5
#>  [8533] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8569] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2
#>  [8605] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [8641] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [8677] 1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 1 5 5 5 1 4 4 1 1 1 5 5 5 5 5 5 1 5 5 5 5
#>  [8713] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8749] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [8785] 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [8821] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 4 4 4 4
#>  [8857] 4 4 4 4 4 1 1 1 1 4 4 4 4 4 4 4 5 5 1 1 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8893] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [8929] 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [8965] 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [9001] 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [9037] 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [9073] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [9109] 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4
#>  [9145] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [9181] 4 4 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [9217] 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [9253] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2
#>  [9289] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [9325] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1
#>  [9361] 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5
#>  [9397] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [9433] 5 5 5 5 5 5 5 5 5 2 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [9469] 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [9505] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4
#>  [9541] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [9577] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2
#>  [9613] 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [9649] 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [9685] 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4 4
#>  [9721] 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#>  [9757] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 5 5 5 5 5 5 5 2 2
#>  [9793] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 2 4 4 4 4 4 4 4
#>  [9829] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 4 4 4
#>  [9865] 4 4 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 5 5 5 5 5 5 5 5
#>  [9901] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 5 5 5
#>  [9937] 5 5 5 5 5 5 5 5 5 5 2 5 5 2 2 2 2 2 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2
#>  [9973] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [10009] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 4 4 4 4 1 1 1 1 1 1 1
#> [10045] 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4 4 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [10081] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [10117] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [10153] 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [10189] 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4
#> [10225] 4 4 4 4 4 4 4 4 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [10261] 5 5 5 5 5 5 5 5 5 5 5 5 4 2 2 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2
#> [10297] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [10333] 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [10369] 4 4 1 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 1 5
#> [10405] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [10441] 5 5 4 2 2 2 5 2 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [10477] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4
#> [10513] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 4 1 4 4 4 4
#> [10549] 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 4 1 5 5 5 5 5 5 5 5 5 5
#> [10585] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 5
#> [10621] 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [10657] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [10693] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 4 4 4 4 4 4 1 1 1 1 1 1 1 1
#> [10729] 1 1 1 1 1 1 1 4 4 4 4 4 4 4 1 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [10765] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 5 2 2 2 2 2 2 2 2 2
#> [10801] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [10837] 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [10873] 4 4 4 4 4 4 4 4 1 1 4 4 1 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4
#> [10909] 4 4 4 4 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [10945] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [10981] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [11017] 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1
#> [11053] 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 4 4 1 1 5 5 5 5
#> [11089] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [11125] 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [11161] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4
#> [11197] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 4 4 4 4 1 4 4 4 4
#> [11233] 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [11269] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 5 5 5 5 5 5 5 5 5 2 2 2
#> [11305] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [11341] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [11377] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 4 4 1 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1
#> [11413] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 1 5 5 5 5
#> [11449] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2
#> [11485] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [11521] 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [11557] 4 4 4 4 1 4 4 4 4 1 1 1 1 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [11593] 1 1 1 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
#> [11629] 5 5 5 5 5 5 5 2 2 5 5 3 5 5 5 2 2 2 2 2 2 2 2 2 5 5 2 2 2 2 2 2 2 2 2 2
#> [11665] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4
#> [11701] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 4 4 4 1
#> [11737] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 5 1 5 5
#> [11773] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2
#> [11809] 5 5 3 5 5 5 2 2 2 2 2 2 2 2 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [11845] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4
#> [11881] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 4 1 1 1 1 1 1 1 1 1 1 1
#> [11917] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 1 5 5 3 3 3 5 5 5 5 5
#> [11953] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 5 1 5 5 5 2 2 2 2
#> [11989] 2 2 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [12025] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [12061] 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [12097] 1 1 1 1 1 1 1 1 1 1 3 3 1 1 3 3 5 5 3 3 3 3 5 5 1 1 3 3 5 5 5 5 5 5 5 5
#> [12133] 5 5 5 5 5 5 4 4 4 4 2 2 2 2 2 2 2 5 1 1 5 5 5 2 2 2 2 2 5 5 5 5 5 5 2 2
#> [12169] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [12205] 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [12241] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [12277] 3 3 1 3 3 3 5 5 3 3 3 5 5 5 3 3 3 3 3 1 1 1 1 5 5 5 5 5 5 5 5 5 5 2 2 2
#> [12313] 2 2 2 2 2 2 2 2 1 3 5 5 5 5 2 2 2 2 2 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2
#> [12349] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4
#> [12385] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1
#> [12421] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 1 1 3 3 1 3 3 3 1 5 3 3 3
#> [12457] 3 3 3 3 3 3 1 3 3 3 5 5 5 1 5 5 5 5 5 1 5 5 5 4 2 2 2 2 2 2 2 2 2 2 2 5
#> [12493] 3 5 5 5 2 2 2 2 2 2 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [12529] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [12565] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [12601] 1 1 1 1 1 1 1 1 1 3 3 3 3 3 1 3 1 1 3 3 3 3 1 1 3 3 3 3 3 3 3 1 1 3 3 3
#> [12637] 3 1 1 3 3 1 1 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5 5 5 1 5 2 2 2 2 2
#> [12673] 2 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [12709] 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [12745] 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [12781] 1 3 3 3 3 3 1 1 3 3 3 3 5 3 3 3 3 1 1 3 3 3 3 1 3 3 3 3 3 3 3 3 1 5 5 1
#> [12817] 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5 3 1 1 5 5 2 2 2 2 2 2 5 1 1 5 2 2 2 2
#> [12853] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [12889] 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1
#> [12925] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 1 3 3 3 1 3 1 1 3
#> [12961] 3 3 3 3 3 3 3 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 1 5 5 5 5 5 2 2 2 2 2 2 2 2
#> [12997] 2 2 2 2 2 2 2 5 3 3 3 5 2 2 2 2 2 2 2 1 5 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [13033] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4
#> [13069] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [13105] 1 1 1 1 1 1 3 1 1 3 3 3 3 3 1 3 3 3 3 3 3 3 3 3 1 1 1 3 3 3 3 3 1 3 3 3
#> [13141] 5 5 1 3 3 3 3 3 3 3 3 3 3 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5 3
#> [13177] 3 5 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [13213] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [13249] 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1
#> [13285] 1 3 3 3 3 3 3 3 3 3 3 3 3 1 1 1 3 3 3 3 1 1 5 1 3 3 5 5 5 5 1 3 3 3 3 3
#> [13321] 1 1 1 1 1 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5 5 5 2 2 2 2 2 2 2 2
#> [13357] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [13393] 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [13429] 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3
#> [13465] 3 3 1 1 3 3 3 1 3 1 1 1 5 5 1 1 1 3 3 1 1 1 3 3 3 3 1 5 1 1 5 2 2 2 2 2
#> [13501] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [13537] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4
#> [13573] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1
#> [13609] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 3 3 3 3 3 3 3 3 1 3 1 1 3 3 3 3 1 1
#> [13645] 1 1 3 1 5 5 1 3 3 1 1 3 5 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [13681] 2 2 2 2 2 2 1 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [13717] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4
#> [13753] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [13789] 1 1 1 1 1 1 1 1 3 1 3 3 1 1 3 3 1 3 1 3 3 3 3 3 3 3 3 3 3 3 1 5 1 3 3 1
#> [13825] 3 3 1 1 1 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5 3 5 5 2
#> [13861] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [13897] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [13933] 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [13969] 1 1 1 1 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 1 3 1 3 3 3 3 3 3 3 1 1 1 1 1 1 5
#> [14005] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 3 1 1 1 5 1 5 2 2 2 2 2 2 2
#> [14041] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [14077] 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 1
#> [14113] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3
#> [14149] 3 3 3 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 3 3 3 3 3 1 1 1 2 2 2 2 2 2 2 2
#> [14185] 2 2 2 2 2 2 2 2 2 2 2 3 1 1 3 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [14221] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4
#> [14257] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 1 1 1 1 1 1 1 1 1
#> [14293] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 1 3 3 1 3 1 1 1 1 3 3
#> [14329] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [14365] 2 3 3 3 1 1 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [14401] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [14437] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [14473] 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 1 3 1 1 3 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3
#> [14509] 3 3 3 3 3 3 3 3 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2 2 2
#> [14545] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [14581] 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [14617] 4 4 4 4 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [14653] 1 1 1 3 1 1 1 1 3 1 1 1 1 1 1 1 3 3 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1
#> [14689] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2
#> [14725] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [14761] 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 2 2
#> [14797] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [14833] 1 1 1 1 1 1 3 3 1 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 2 2 2 2 2 2 2 2 2
#> [14869] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [14905] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4
#> [14941] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 2 2 2 1 1 1 2 2 1 1 1 1
#> [14977] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [15013] 1 1 1 3 3 3 3 3 3 3 3 3 3 3 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [15049] 2 2 2 1 1 1 1 2 2 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [15085] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [15121] 4 4 4 4 4 4 4 4 4 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1
#> [15157] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3
#> [15193] 3 3 3 3 3 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 3 1 2 1 1
#> [15229] 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [15265] 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [15301] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [15337] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 2 2 2 2
#> [15373] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 3 3 1 1 1 2 2 2 2 2 2 2 2 2 2
#> [15409] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4
#> [15445] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 2 2 1 1 2 2
#> [15481] 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [15517] 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [15553] 2 2 2 2 2 2 2 2 1 3 3 3 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [15589] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [15625] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 2 2 1 1 2 2 2 1 1 1 1 1 1 1 1 1
#> [15661] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [15697] 3 3 3 3 3 3 3 3 3 3 3 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1
#> [15733] 3 3 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [15769] 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [15805] 4 4 4 4 4 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1
#> [15841] 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3
#> [15877] 3 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 3 3 1 2 2 2 2 2 2
#> [15913] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [15949] 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2
#> [15985] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 2 2 2 2 2 2 2 2 2 1 1 1
#> [16021] 2 2 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 1 3 3 3 3 3 3 3 3 3 1 1 1 2 2 2 2 2 2
#> [16057] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 3 3 3 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [16093] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4
#> [16129] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [16165] 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 1 1 1 2 1 1 1 1 1 1 1 2 2
#> [16201] 1 1 1 1 1 3 3 1 3 3 3 3 3 3 3 3 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [16237] 2 2 2 2 2 1 3 1 3 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [16273] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [16309] 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [16345] 2 2 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 1 2 1 2 2 2 2 2 1 1 1 1 1 1 1 1 1 3 3
#> [16381] 3 3 3 3 3 3 3 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 3 3 3 1 1
#> [16417] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [16453] 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> [16489] 4 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [16525] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 3 3 1 3 3 3 3 3 3 3 1 2
#> [16561] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 3 3 3 1 1 1 2 2 2 2 2 2 2 2 2 2
#> [16597] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4
#> [16633] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 2 2 2 2
#> [16669] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2
#> [16705] 2 2 1 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 1 1 1 2 2 2 2 2 2 2 2 2 2
#> [16741] 2 2 2 2 2 2 2 2 1 3 3 3 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#> [16777] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 4
#> [16813] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
#> 
#> Within cluster sum of squares by cluster:
#> [1] 3516.041 4545.222 3785.104 5524.327 7478.029
#>  (between_SS / total_SS =  70.5 %)
#> 
#> Available components:
#> 
#> [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
#> [6] "betweenss"    "size"         "iter"         "ifault"

Copying a raster structure

clus_raster <- co$tmax
values(clus_raster) <- NA

Assign values

clus_raster[idx] <- E$cluster
plot(clus_raster, col = RColorBrewer::brewer.pal(5, "Spectral"))

Merging across data sources

# Get elevations data
elev = elevatr::get_elev_raster(AOI, z = 5) %>% 
  crop(AOI) |> 
  rast()

# Align Raster extents and resolutions
elev = project(elev, co$ppt)

# Extract Values
values = c(co$ppt, elev) %>% values()

# Prep data
idx = which(!apply(is.na(values), 1, any))
v = na.omit(values)
vs = scale(v)

# Cluster
E = kmeans(vs, 5, iter.max = 100)

clus_raster = elev
values(clus_raster) = NA
clus_raster[idx] <- E$cluster

par(mfrow = c(2,1))
plot(elev)
plot(co$ppt)

plot(clus_raster, col = RColorBrewer::brewer.pal(5, "Spectral"))

Clustering at larger scales

counties =  AOI::aoi_get(state = "conus", county = "all")
params = c("tmax", "tmin", "ppt", "srad")

dat = climateR::getTerraClim(counties, params, startDate = "2018-06-01") %>% 
  unlist() |> 
  rast() %>% 
  setNames(params) |> 
  exactextractr::exact_extract(counties, "mean", progress = FALSE)

dat = scale(dat)

counties$clust8 = kmeans(dat, 8)$cluster

plot(st_transform(counties["clust8"], 5070), border = NA)

Extra! Animating Rasters…

Install and load the gifski package

save_gif: combines many individual plots
A for loop build the plots
The plot is what we have been doing all along (if you want a ggplot you must print the object!)
gif_file: the path to save the image
width/height: the image dimensions
delay: the pause between frames
loop: should the gif play over and over?

library(gifski)

save_gif(
  {for(i in 1:nlyr(prcp$ppt)) {
      plot(prcp$ppt[[i]], col = blues9, 
           legend = FALSE, 
           main = names(prcp$ppt)[i])}
  }, 
  gif_file = "images/ppt.gif", 
  width = 800, height = 600, delay = .33, loop = TRUE)

Recap: sf and spatial concepts

terra

Recap: R data structures

Spatial Extent

We can extract bounding box coordinates with st_bbox

There is a method for creating an sfc form a bbox object

Result:

Extents can be discritized in a number of ways:

What makes of a regular tesselation

Alternative representation

Equal area from centroid

Raster Model

Many terms mean the same thing …

Photos and Computers …

Aerial Imagery (really just a photo 😄)

What is stored in these cells?

Categorical Values (integer/factor)

Continuous Values (numeric)

Spectral Values

Any digital image contains an RBG channel for color:

Pure RBG

RGB and bytes/bits

Why do we care?

Resolution drives image clarity (granulairty)

Raster images seek to discritize the real world into cell-based values

All rasters have an extent!

Implicity Coordinates

Raster Data in R

Extract Diminisionality …

Raster Values

Assigning a values as a vector

Raster is S4

Access via Function

Multi-layers

Remember our Array data?

An array is a single file so it can be “rasterized”

Core Raster Principles

1. Bounding Box / Extents

2. Extent

3. Discretization

So,

General Process

Find elevation data for Fort Collins:

Raster Values Continuity

Raster Values

All image files are the same!

Raster Algebra

Raster Algebra

For example:

Replacement

Modifying a raster

Modifying the underlying data:

What is mask doing?

Crop or/and mask

Aggregate and disaggregate

app

Read in the saved raster file

Create a conditional (threshold) mask

Results

Multiply cell-wise

Reclassify

Real Example: Classify Rainfall Regions of Colorado

Object/Field Interaction

Lets find the longest river segment IN our extent

Value Extraction

Sinuosity

River Slope:

River profile

Map Algebra

Local

Focal

Focal

Example: Focal

Results

What did we do?

Zonal

Zonal Statisics (More advanced)

What about the US?

Global

Summary Values

`terra`

We can extract bounding box coordinates with `st_bbox`

There is a method for creating an sfc form a `bbox` object

What is `mask` doing?

`app`

Why not just `mean()`