ECON6027 Cheatsheet
Gertrude Cheung
2023-11-09
Packages & libraries
Spatial
install.packages(“sf”, dependencies = T) install.packages(“sp”,
dependencies = T) install.packages(“spatstat”, dependencies = T)
install.packages(“gstat”, dependencies = T) install.packages(“terra”,
repos=‘https://rspatial.r-universe.dev’)
install.packages(“raster”, repos=‘https://rspatial.r-universe.dev’)
install.packages(“dismo”, dependencies = T)
install.packages(“tmap”, dependencies = T) install.packages(“automap”,
dependencies = T)
Non-spatial
install.packages(“tidyverse”, dependencies = T)
Lecture 8a
1 Nearest neighbour interpolation
Fulmaris glacialis data
Load the data:
## year x y depth coast fulmar
## 1 1998 614192.3 5875490 6 3.445307 0
## 4 1998 613150.8 5872947 6 3.445307 0
## 5 1998 619644.7 5888800 9 3.042531 0
## 6 1998 608858.6 5839623 16 3.845921 0
## 7 1998 556583.0 5738519 1 2.701462 0
## 8 1998 553423.7 5735916 3 1.807041 0Create a Voronoi tessallation for the fulmar dataframe to provide nearest neighbours interpolation
- Convert dataset to sf object:
- Extract the 1999 data and plot the Fulmar Sighting Density:
fulmar.sf.1999 = fulmar.sf[fulmar.sf$year == 1999, ]
tm_shape(fulmar.sf.1999) +
tm_dots(col = "fulmar",
palette = "Set1") +
tm_layout(main.title = "Fulmar Sighting Density (1999)",
legend.outside = T)
- Convert to sp object as ‘dismo’ only works on sp objects. Plot the Voronoi tessellation:
fulmar.sp.1999 = as(fulmar.sf.1999,
"Spatial")
fulmar.1999.voronoi = voronoi(fulmar.sp.1999,
ext = st_bbox(fulmar.sf.1999))
class(fulmar.1999.voronoi)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"## Simple feature collection with 729 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 476209.6 ymin: 5695049 xmax: 735071.7 ymax: 6143369
## Projected CRS: +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs
## First 10 features:
## year depth coast fulmar geometry
## 1 1999 39 74.32973 0.00000 POLYGON ((656993.9 5988304,...
## 2 1999 39 66.93047 0.00000 POLYGON ((621770.4 5972015,...
## 3 1999 40 87.58816 6.83660 POLYGON ((587962.6 5975436,...
## 4 1999 33 66.16789 0.00000 POLYGON ((519371.2 5800328,...
## 5 1999 26 51.32900 0.00000 POLYGON ((593794.3 5935591,...
## 6 1999 21 33.62153 0.00000 POLYGON ((512644.7 5738338,...
## 7 1999 44 85.68657 28.42134 POLYGON ((585502.3 5998381,...
## 8 1999 36 140.23980 0.00000 POLYGON ((527315.8 5992129,...
## 9 1999 31 63.64503 0.00000 POLYGON ((486951.4 5739733,...
## 10 1999 25 41.07539 0.00000 POLYGON ((596477.3 5923979,...tm_shape(fulmar.1999.voronoi.sf) +
tm_polygons(col = "fulmar",
palette = "Set3") +
tm_layout(main.title = "Voronoi tessellation for Fulmar Density",
legend.outside = T)
Observations of the Voronoi map:
- Areas of the polygons vary with the density of points (for internal points)
- Edge points have polygons of infinite area (trimmed to an enclosing rectangle)
- Edge polygons somewhat dominate the interpolation visually
- Increase in density of sightings towards the North-East direction
- Discontinuous nature of the interpolated surfaces which creates a distorting effect of the variety of shapes and sizes, thus harder to identify a subtler pattern
- Interpolated surfaces are discontinuous which is not usual for geospatial data in general
California Ozone Layer
Load the data:
## Reading layer `ca_ozone_pts' from data source
## `/cloud/project/econ6027 cheatsheet/ca_ozone/ca_ozone_pts.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 193 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -2301588 ymin: -355613.7 xmax: -1795559 ymax: 780085.2
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic## Simple feature collection with 193 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -2301588 ymin: -355613.7 xmax: -1795559 ymax: 780085.2
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
## First 10 features:
## LATITUDE LONGITUDE ELEVATION OZONE geometry
## 1 39.1447 -123.2065 194 0.04650 POINT (-2298092 515557.4)
## 2 39.4030 -123.3491 420 0.04969 POINT (-2301588 546772.7)
## 3 37.7661 -122.3978 5 0.05000 POINT (-2273948 347691.4)
## 4 37.9508 -122.3569 23 0.05799 POINT (-2264847 366623.2)
## 5 36.6986 -121.6354 36 0.05860 POINT (-2241776 214412.1)
## 6 35.3668 -120.8450 18 0.05889 POINT (-2212812 51979.78)
## 7 37.8014 -122.2672 7 0.05990 POINT (-2261907 348380.1)
## 8 37.0120 -122.1929 0 0.06060 POINT (-2279746 261490.2)
## 9 38.4436 -122.7092 49 0.06350 POINT (-2279005 428141)
## 10 41.7293 -122.6354 802 0.06369 POINT (-2171057 780085.2)##
## TRUE
## 193## Reading layer `ca_outline' from data source
## `/cloud/project/econ6027 cheatsheet/ca_ozone/ca_outline.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2354901 ymin: -364877.4 xmax: -1646015 ymax: 845703.7
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic## Simple feature collection with 1 feature and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2354901 ymin: -364877.4 xmax: -1646015 ymax: 845703.7
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
## AREA STATE_NAME STATE_FIPS SUB_REGION STATE_ABBR
## 1 408634904360 California 06 Pacific CA
## geometry
## 1 MULTIPOLYGON (((-2200259 37...##
## TRUE
## 1- Check and set the CRS if necessary:
## Coordinate Reference System:
## User input: USA_Contiguous_Albers_Equal_Area_Conic
## wkt:
## PROJCRS["USA_Contiguous_Albers_Equal_Area_Conic",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["USA_Contiguous_Albers_Equal_Area_Conic",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",37.5,
## 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["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["USA - CONUS - onshore"],
## BBOX[24.41,-124.79,49.38,-66.91]],
## ID["ESRI",102003]]## [1] "+proj=aea +lat_0=37.5 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"## Coordinate Reference System:
## User input: USA_Contiguous_Albers_Equal_Area_Conic
## wkt:
## PROJCRS["USA_Contiguous_Albers_Equal_Area_Conic",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["USA_Contiguous_Albers_Equal_Area_Conic",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",37.5,
## 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["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["USA - CONUS - onshore"],
## BBOX[24.41,-124.79,49.38,-66.91]],
## ID["ESRI",102003]]## [1] "+proj=aea +lat_0=37.5 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"- Convert to sp object as ‘dismo’ only works on sp objects. Plot the Voronoi tessellation:
ozone.sp = as(ozone.sf,
"Spatial")
ozone.voronoi = voronoi(ozone.sp,
ext = st_bbox(ca.sf))
class(ozone.voronoi)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"## Simple feature collection with 193 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -2354901 ymin: -364877.4 xmax: -1646015 ymax: 845703.7
## Projected CRS: +proj=aea +lat_0=37.5 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
## First 10 features:
## LATITUDE LONGITUDE ELEVATION OZONE geometry
## 1 39.1447 -123.2065 194 0.04650 POLYGON ((-2260182 535605.8...
## 2 39.4030 -123.3491 420 0.04969 POLYGON ((-2269084 625328.9...
## 3 37.7661 -122.3978 5 0.05000 POLYGON ((-2268422 356688.4...
## 4 37.9508 -122.3569 23 0.05799 POLYGON ((-2250035 365300.5...
## 5 36.6986 -121.6354 36 0.05860 POLYGON ((-2231077 228321.9...
## 6 35.3668 -120.8450 18 0.05889 POLYGON ((-2210740 71650.62...
## 7 37.8014 -122.2672 7 0.05990 POLYGON ((-2252019 359332.6...
## 8 37.0120 -122.1929 0 0.06060 POLYGON ((-2271587 276242.1...
## 9 38.4436 -122.7092 49 0.06350 POLYGON ((-2250589 432733.5...
## 10 41.7293 -122.6354 802 0.06369 POLYGON ((-1701646 845703.7...tm_shape(ozone.voronoi.sf) +
tm_polygons(col = "OZONE",
palette = "Set3") +
tm_layout(main.title = "Voronoi tessellation for California Ozone Layer",
legend.outside = T)
2 Inverse Distance Weighting Interpolation
idw function works with sp/sf objects (a CRS must be set) as well as a set of unsampled locations to be interpolated.
Load the data and convert it to an sf data and set the crs:
This example also utilises data from the fulmar dataset.
- Derive the idw estimate:
## Where k = 1
(fulmar.idw1 = idw(formula = fulmar ~ 1,
locations = fulmar.sf.1999,
newdata = ncp.grid.sf,
idp = 1))## [inverse distance weighted interpolation]## Simple feature collection with 2297 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 466500 ymin: 5699000 xmax: 736500 ymax: 6129000
## Projected CRS: WGS 84 / UTM zone 31N
## First 10 features:
## var1.pred var1.var geometry
## 1 1.551593 NA POINT (511500 6129000)
## 2 1.530795 NA POINT (516500 6129000)
## 3 1.480545 NA POINT (521500 6129000)
## 4 1.436883 NA POINT (526500 6129000)
## 5 1.554843 NA POINT (531500 6129000)
## 6 1.629213 NA POINT (536500 6129000)
## 7 1.657564 NA POINT (541500 6129000)
## 8 1.672028 NA POINT (546500 6129000)
## 9 1.739322 NA POINT (551500 6129000)
## 10 1.841795 NA POINT (556500 6129000)## Where k = 2
(fulmar.idw2 = idw(formula = fulmar ~ 1,
locations = fulmar.sf.1999,
newdata = ncp.grid.sf,
idp = 2))## [inverse distance weighted interpolation]## Simple feature collection with 2297 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 466500 ymin: 5699000 xmax: 736500 ymax: 6129000
## Projected CRS: WGS 84 / UTM zone 31N
## First 10 features:
## var1.pred var1.var geometry
## 1 1.3392497 NA POINT (511500 6129000)
## 2 1.0882239 NA POINT (516500 6129000)
## 3 0.6618765 NA POINT (521500 6129000)
## 4 0.4245152 NA POINT (526500 6129000)
## 5 1.0926721 NA POINT (531500 6129000)
## 6 1.6610271 NA POINT (536500 6129000)
## 7 1.6246147 NA POINT (541500 6129000)
## 8 1.4646454 NA POINT (546500 6129000)
## 9 1.6321421 NA POINT (551500 6129000)
## 10 3.0618283 NA POINT (556500 6129000)- Plot the two idw plots:
fulmar.buff.1 = st_buffer(fulmar.idw1,
dist = 3500)
fulmar.buff.2 = st_buffer(fulmar.idw2,
dist = 3500)
fulmar.idw.1.map = tm_shape(fulmar.buff.1) +
tm_fill(col = "var1.pred",
palette = "YlGnBu") +
tm_layout(main.title = "Fulmar density estimation, idp = 1",
main.title.size = 0.9)
fulmar.idw.2.map = tm_shape(fulmar.buff.2) +
tm_fill(col = "var1.pred",
palette = "YlGnBu",
breaks = c(0, 1, 2, 3, 4, 5, 6, Inf)) +
tm_layout(main.title = "Fulmar density estimation, idp = 2",
main.title.size = 0.9)
tmap_arrange(fulmar.idw.1.map,
fulmar.idw.2.map)
3 Kriging
Load the data and convert it to an sf data and set the crs:
This example also utilises data from the fulmar dataset.
data(fulmar)
fulmar.1999 = fulmar[fulmar$year == 1999, ]
fulmar.1999.sf = st_as_sf(fulmar.1999,
coords = c("x", "y"),
crs = 32631)
st_crs(fulmar.1999.sf); st_crs(fulmar.1999.sf)$proj; st_crs(fulmar.1999.sf)$proj4string## Coordinate Reference System:
## User input: EPSG:32631
## wkt:
## PROJCRS["WGS 84 / UTM zone 31N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 31N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",3,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["World - N hemisphere - 0°E to 6°E - by country"],
## BBOX[0,0,84,6]],
## ID["EPSG",32631]]## [1] "utm"## [1] "+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs"3.1 Construct the semi-variogram
- Estimate the average squared differences between observations for distance bands:
Giving boundaries is options but gives better control. Boundaries provide the distance bands. Given that the maximum distance in the “northings” direction is approximately 261.5km, we will choose a distance of 260km in bands of 5km. The number of lags will then be lags = (260/5) + 1 = 53.
- Determine the appropriate semi-variogram through a trial and error method and the automatic method:
Trial and error approach:
- Have to identify the value of sill, range and nugget from the initial plot.
## short long
## 1 Nug Nug (nugget)
## 2 Exp Exp (exponential)
## 3 Sph Sph (spherical)
## 4 Gau Gau (gaussian)
## 5 Exc Exclass (Exponential class/stable)
## 6 Mat Mat (Matern)
## 7 Ste Mat (Matern, M. Stein's parameterization)
## 8 Cir Cir (circular)
## 9 Lin Lin (linear)
## 10 Bes Bes (bessel)
## 11 Pen Pen (pentaspherical)
## 12 Per Per (periodic)
## 13 Wav Wav (wave)
## 14 Hol Hol (hole)
## 15 Log Log (logarithmic)
## 16 Pow Pow (power)
## 17 Spl Spl (spline)
## 18 Leg Leg (Legendre)
## 19 Err Err (Measurement error)
## 20 Int Int (Intercept)
fulmar.fvgm.exp = fit.variogram(fulmar.evgm,
vgm(model = "Ste",
psill = 10,
range = 100000,
nugget = 2))
plot(fulmar.evgm,
model = fulmar.fvgm.exp)
Automatic approach:
Automap currently only accepts projected sp objects.
fulmar.1999.sp = as(fulmar.1999.sf,
"Spatial")
(fulmar.autovgm = autofitVariogram(fulmar ~ 1,
input_data = fulmar.1999.sp))## $exp_var
## np dist gamma dir.hor dir.ver id
## 1 393 1784.391 2.901021 0 0 var1
## 2 1099 5870.347 1.786798 0 0 var1
## 3 1269 8731.925 4.297250 0 0 var1
## 4 2514 13875.243 3.645847 0 0 var1
## 5 2976 19231.834 3.515168 0 0 var1
## 6 3950 24431.539 5.504788 0 0 var1
## 7 16264 36637.791 6.044911 0 0 var1
## 8 20485 54703.306 7.875939 0 0 var1
## 9 37198 77231.559 12.431239 0 0 var1
## 10 36291 103969.539 14.249815 0 0 var1
## 11 32997 131163.796 11.983459 0 0 var1
## 12 37915 162590.448 12.434498 0 0 var1
##
## $var_model
## model psill range
## 1 Nug 2.798431 0.00
## 2 Gau 10.745451 59275.79
##
## $sserr
## [1] 0.000100295
##
## attr(,"class")
## [1] "autofitVariogram" "list"
We will use the Gaussian model as suggested by automap.
3.2 Ordinary Kriging Estimate
Load and prepare the data:
- Derive the Kriging estimate:
If no variogram is specificed to the argument “model”, the function defaults to IDW interpolation wth idp = 2.
ncp.grid.krig.est = krige(formula = fulmar ~ 1,
locations = fulmar.1999.sf,
newdata = ncp.grid.sf,
model = fulmar.fvgm.exp) ## How to use autovgm?## [using ordinary kriging]## var1.pred var1.var geometry
## Min. :-0.08568 Min. :1.948 POINT :2297
## 1st Qu.: 0.05974 1st Qu.:2.615 epsg:32631 : 0
## Median : 0.94513 Median :3.026 +proj=utm ...: 0
## Mean : 1.94352 Mean :3.330
## 3rd Qu.: 2.73492 3rd Qu.:3.809
## Max. :26.72930 Max. :8.550- Plot the Kriging estimates:
- Negative frequencies must be “eliminated” in the plot.
ncp.grid.krig.buff = st_buffer(ncp.grid.krig.est,
dist = 3500)
tm_shape(ncp.grid.krig.buff) +
tm_fill(col = "var1.pred",
palette = "YlGnBu",
breaks = c(-Inf, 2, 4, 6, 8, 10, 12, 14, Inf)) +
tm_layout(main.title = "Fulmar density estimation, Ordinary Kriging",
main.title.size = 0.8)
tm_shape(ncp.grid.krig.buff) +
tm_fill(col = "var1.var",
palette = "Blues") +
tm_layout(main.title = "Fulmar variance of Ordinary Kriging Estimate",
main.title.size = 0.8)
3.3 Simple Kriging Estimate
The Simple Kriging approach will derive the same output as the Ordinary Kriging Approach when \(\beta_i\) is specified in the argument.
- Run the OLS test without specifying the value of \(\beta_i\) in the argument and extracting the intercept:
##
## Call:
## lm(formula = fulmar ~ 1, data = fulmar.1999.sf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.116 -1.116 -1.116 -1.116 45.370
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1165 0.1247 8.951 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.368 on 728 degrees of freedom- Derive the Kriging estimate:
ncp.grid.s.krig.est = krige(formula = fulmar ~ 1,
locations = fulmar.1999.sf,
newdata = ncp.grid.sf,
beta = ncp.grid.s.ols.int,
model = fulmar.fvgm.exp)## [using simple kriging]- Plot the Kriging estimates:
- Negative frequencies must be “eliminated” in the plot.
ncp.grid.s.krig.buff = st_buffer(ncp.grid.s.krig.est,
dist = 3500)
tm_shape(ncp.grid.s.krig.buff) +
tm_fill(col = "var1.pred",
palette = "YlGnBu",
breaks = c(-Inf, 2, 4, 6, 8, 10, 12, 14, Inf)) +
tm_layout(main.title = "Fulmar density estimation, Simple Kriging",
main.title.size = 0.8)
There are similarity between the two estimates as seen by the plots. This reiterates the fact that all kriging techniques honour data characteristics, preserves the mean (trend) and preserve the spatial autocorrelation structure. Ordinary Kriging is a safe bet in many situations.
3.4 Universal Kriging Estimate in R
If we specify a regression equation in the krige function, then the output is the universal kriging estimate.
- Run the AIC test to see which model is a better fit for the dataset:
fulmar.mdl1.aic = AIC(lm(formula = fulmar ~ depth + coast,
data = fulmar.1999.sf))
fulmar.mdl2.aic = AIC(lm(formula = fulmar ~ depth,
data = fulmar.1999.sf))
fulmar.mdl3.aic = AIC(lm(formula = fulmar ~ coast,
data = fulmar.1999.sf))
fulmar.mdl4.aic = AIC(lm(formula = fulmar ~ 1,
data = fulmar.1999.sf))Of the four models, Model 1 has the smallest AIC (3703.98), therefore it is selected.
- Derive the Kriging estimate:
ncp.grid.u.krig.est = krige(formula = fulmar ~ depth + coast,
locations = fulmar.1999.sf,
newdata = ncp.grid.sf,
model = fulmar.fvgm.exp)## [using universal kriging]- Plot the Kriging estimates:
- Negative frequencies must be “eliminated” in the plot.
ncp.grid.u.krig.buff = st_buffer(ncp.grid.u.krig.est,
dist = 3500)
tm_shape(ncp.grid.u.krig.buff) +
tm_fill(col = "var1.pred",
palette = "YlGnBu",
breaks = c(-Inf, 2, 4, 6, 8, 10, 12, 14, Inf)) +
tm_layout(main.title = "Fulmar density estimation, Universal Kriging",
main.title.size = 0.8)
Applying the Universal Kriging approach, the output is still similar to the Ordinary Kriging approach.