AST 511
Humayara Turin
2025-09-30
1 Spatial Data
1.0.1 Libraries
library(sf) #For working with spatial data
library(tidyverse) #For data wrangling and ggplot2
library(ggspatial) #For spatial data with ggplot2
library(tigris) #For downloading spatial data directly into R1.1 Read in Virginia polygon
setwd("G:/My Drive/Turin/5th year MS/AST 511 [Enviromental and Spatial Statistics]/practical part/Spatial data/Data")
virginia <- st_read("G:/My Drive/Turin/5th year MS/AST 511 [Enviromental and Spatial Statistics]/practical part/Spatial data/Data/VirginiaBoundary.shp") ## Reading layer `VirginiaBoundary' from data source
## `G:\My Drive\Turin\5th year MS\AST 511 [Enviromental and Spatial Statistics]\practical part\Spatial data\Data\VirginiaBoundary.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -83.67541 ymin: 36.54074 xmax: -75.86704 ymax: 39.46601
## Geodetic CRS: NAD83plot(virginia)
as.character(st_geometry_type(virginia))[1]## [1] "POLYGON"plot(virginia$geometry)
1.1.1 Transform CRS
st_crs(virginia)$epsg## [1] 4269virginia.crs <- st_transform(virginia, crs = 26918)1.2 Data frame of cities
cities <- data.frame(
City = c("Richmond", "Virginia Beach", "Norfolk", "Chesapeake", "Arlington", "Charlottesville"),
Latitude = c(37.5413, 36.8529, 36.8508, 36.7682, 38.8797, 38.0293),
Longitude = c(-77.4348, -75.9780, -76.2859, -76.2875, -77.1080, -78.4767))
str(cities)## 'data.frame': 6 obs. of 3 variables:
## $ City : chr "Richmond" "Virginia Beach" "Norfolk" "Chesapeake" ...
## $ Latitude : num 37.5 36.9 36.9 36.8 38.9 ...
## $ Longitude: num -77.4 -76 -76.3 -76.3 -77.1 ...1.2.1 Convert to sf object
cities.sf <-
st_as_sf(cities, coords = c("Longitude", "Latitude"), remove = FALSE, crs = 4326)
plot(cities.sf$geometry)
1.2.2 Transform CRS
st_crs(cities.sf)$epsg## [1] 4326cities.crs <- st_transform(cities.sf, crs = 26918)1.2.3 compare crs and the plot
#Check that the layers have the same CRS
st_crs(virginia.crs) == st_crs(cities.crs)## [1] TRUEggplot() +
#Add Virginia polygon
geom_sf(data = virginia.crs) +
#Add cities layer
geom_sf(data = cities.crs, aes(color = City, label = City), size = 3) +
labs(x = "Longitude", y = "Latitude", title = "Virginia Cities") +
theme_bw()## Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
## : Ignoring unknown aesthetics: label
1.2.4 Get bounding box(the smallest rectangle that contains all geometries) spatial object cities.crs.
cities.extent <- st_bbox(cities.crs)
cities.extent## xmin ymin xmax ymax
## 194831.0 4069931.4 412812.9 4305538.81.2.5 Plot using the packages ggplot and ggspatial
ggplot() +
#Crop Virginia boundary to spatial extent of cities and add Virginia layer
geom_sf(data = st_crop(virginia.crs, cities.extent)) +
#Add cities layer
geom_sf(data = cities.crs, aes(color = City, label = City), size = 3) +
#Add scale bar to bottom left from ggspatial
annotation_scale(location = "bl", height = unit(.25, "cm"),
width = unit(1, "cm"), pad_x = unit(0.3, "in"),
pad_y = unit(0.3, "in")) +
#Add north arrow to bottom left from ggspatial
annotation_north_arrow(height = unit(1, "cm"), width = unit(1, "cm"),
which_north = "true", location = "bl",
pad_x = unit(0.4, "in"), pad_y = unit(0.5, "in")) +
labs(x = "Longitude", y = "Latitude", title = "Virginia Cities") +
theme_bw()## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries## Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
## : Ignoring unknown aesthetics: label## Warning in annotation_scale(location = "bl", height = unit(0.25, "cm"), :
## Ignoring unknown parameters: `width`
1.2.6 Distance from each city to the border of Virginia
dist.border <- st_distance(cities.crs, st_boundary(virginia.crs))
print(dist.border)## Units: [m]
## [,1]
## [1,] 85469.5113
## [2,] 121.7268
## [3,] 11906.0329
## [4,] 20091.7204
## [5,] 3054.4664
## [6,] 82565.95711.2.7 Calculate distance from Charlottesville to each city
dist.cville <-
st_distance(cities.crs, filter(cities.crs, City == "Charlottesville")) %>%
bind_cols(cities.crs)## New names:
## • `` -> `...1`head(dist.cville, n = 3)## ...1 City Latitude Longitude geometry
## 1 106615.8 Richmond 37.5413 -77.4348 POINT (284889.6 4157709)
## 2 256808.3 Virginia Beach 36.8529 -75.9780 POINT (412812.9 4079001)
## 3 233902.9 Norfolk 36.8508 -76.2859 POINT (385359.8 4079093)1.3 county data
Get county data for Virginia from tigris package and transform the CRS
va.counties <- counties(state = "VA") %>%
st_transform(crs = 26918)## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%1.4 Road data
Get road data and transform the CRS
roads <- primary_secondary_roads(state = "VA", year = 2021) %>% st_transform(crs = 26918)## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%1.4.1 Plot with ggplot
ggplot() +
#Crop Virginia counties to cities extent and plot
geom_sf(data = st_crop(va.counties, cities.extent), color = "gray100", fill = "gray88", size = .25) +
#Crop roads to cities extent and plot
geom_sf(data = st_crop(roads, cities.extent), color = "black", size = .3) +
#Plot cities
geom_sf(data = cities.crs, aes(color = City), size = 4) +
labs(x = "Longitude", y = "Latitude", title = "Map of Virginia Cities") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
1.5 Read in Bangladesh Data
bd <- st_read("G:/My Drive/Turin/5th year MS/AST 511 [Enviromental and Spatial Statistics]/practical part/Spatial data/Data/BGD_adm0.shp")## Reading layer `BGD_adm0' from data source
## `G:\My Drive\Turin\5th year MS\AST 511 [Enviromental and Spatial Statistics]\practical part\Spatial data\Data\BGD_adm0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 67 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 88.01057 ymin: 20.74111 xmax: 92.67366 ymax: 26.6344
## Geodetic CRS: WGS 84plot(bd$geometry)
plot(bd)
as.character(st_geometry_type(bd))## [1] "MULTIPOLYGON"1.6 Read in Temperature Data
tem <- read_csv("G:/My Drive/Turin/5th year MS/AST 511 [Enviromental and Spatial Statistics]/practical part/Spatial data/Data/Average_temp_2018.csv")
tem <- tem[,c(2,5,6)]
tem <- st_as_sf(tem, coords = c("Lon", "Lat"), remove = FALSE, crs = 4326)
plot(tem$geometry)
plot(tem)
1.6.1 compare crs and the plot
st_crs(bd)$epsg## [1] 4326st_crs(tem)$epsg## [1] 4326ggplot()+
geom_sf(data = bd)+
geom_sf(data=tem, aes(color=Station, label=Station), size=3)+
labs(x = "Longitude", y = "Latitude", title = "BD Temperature") +
theme_bw()## Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
## : Ignoring unknown aesthetics: label
1.6.2 set labels in BD
ggplot()+
geom_sf(data = bd)+
geom_sf(data=tem, aes(color=Station))+
annotation_scale(location = "bl", height = unit(.25, "cm"),
width = unit(1, "cm"), pad_x = unit(0.3, "in"),
pad_y = unit(0.3, "in")) +
annotation_north_arrow(height = unit(1, "cm"), width = unit(1, "cm"),
which_north = "true", location = "bl",
pad_x = unit(0.4, "in"), pad_y = unit(0.5, "in")) +
labs(x = "Longitude", y = "Latitude", title = "Bangladesh's envermental station") +
theme_bw()
1.6.3 Plot using the packages ggplot and ggspatial
station.extent <- st_bbox(tem)
ggplot() +
geom_sf(data = st_crop(bd, station.extent)) +
geom_sf(data = tem, aes(color = Station, label = Station), size = 3) +
annotation_scale(location = "bl", height = unit(.25, "cm"),
width = unit(1, "cm"), pad_x = unit(0.3, "in"),
pad_y = unit(0.3, "in")) +
annotation_north_arrow(height = unit(1, "cm"), width = unit(1, "cm"),
which_north = "true", location = "bl",
pad_x = unit(0.4, "in"), pad_y = unit(0.5, "in")) +
labs(x = "Longitude", y = "Latitude", title = "BD Stations") +
theme_bw()## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries## Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
## : Ignoring unknown aesthetics: label## Warning in annotation_scale(location = "bl", height = unit(0.25, "cm"), :
## Ignoring unknown parameters: `width`
1.6.4 Calculate distance from boundary to each city
dist.dhaka <- st_distance(tem, st_boundary(bd))
head(dist.dhaka)## Units: [m]
## [,1]
## [1,] 6342.159
## [2,] 2491.445
## [3,] 4080.144
## [4,] 56254.037
## [5,] 8151.806
## [6,] 4298.5341.6.5 Calculate distance from Dhaka to each city
dist.dhaka <- st_distance(tem, filter(tem, Station == "Dhaka")) %>%
bind_cols(tem)
head(dist.dhaka, n=4)## ...1 Station Lat Lon geometry
## 1 215239.0 Ambagan(Ctg) 22.35000 91.81700 POINT (91.817 22.35)
## 2 113061.2 Barisal 22.75000 90.36667 POINT (90.36667 22.75)
## 3 123504.7 Bhola 22.68333 90.65000 POINT (90.65 22.68333)
## 4 158507.8 Bogra 24.85000 89.36667 POINT (89.36667 24.85)1.7 Bangladesh district boundaries
bd.districts <- st_read("G:/My Drive/Turin/5th year MS/AST 511 [Enviromental and Spatial Statistics]/practical part/Spatial data/Data/BGD_adm1.shp") %>%
st_transform(crs = 32646) # UTM Zone 46N## Reading layer `BGD_adm1' from data source
## `G:\My Drive\Turin\5th year MS\AST 511 [Enviromental and Spatial Statistics]\practical part\Spatial data\Data\BGD_adm1.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 88.01057 ymin: 20.74111 xmax: 92.67366 ymax: 26.6344
## Geodetic CRS: WGS 841.8 Load roads shapefile (e.g., from Geofabrik OSM)
bd.roads <- st_read("path/to/bangladesh-roads.shp") %>%
st_transform(crs = 32646)ggplot() +
# Crop and plot Bangladesh districts
geom_sf(data = st_crop(bd.districts, stations.extent),
color = "gray100", fill = "gray88", size = 0.25) +
# Crop and plot roads
geom_sf(data = st_crop(bd.roads, stations.extent),
color = "black", size = 0.3) +
# Plot stations (cities)
geom_sf(data = tem, aes(color = Station), size = 4) +
labs(x = "Longitude", y = "Latitude",
title = "Map of Bangladesh Monitoring Stations") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))2 Geostatistics
2.0.1 Loading Libraries
library(readr)
library(geoR)
library(gstat)
library(scatterplot3d)
library(sf)
library(ggplot2)2.1 Coalash Data
data("coalash")2.2 Large scale component
plot(coalash$x,coalash$y,pch=4,cex=0.75,xlab="x",ylab="y")
scatterplot3d(coalash$x,coalash$y, coalash$coalash, xlab = "x", ylab = "y", zlab = "coal")
The 3D scatter plot illustrates how coal ash values change across the two spatial dimensions, x and y. Each point represents an observation, with its position on the horizontal axes showing location and its height indicating the coal ash value. If the points appear to align along a plane or slope, it suggests a spatial trend in coal ash distribution, whereas a more random scatter indicates little to no spatial relationship.
2.2.1 Text instead of points
plot(coalash$x,coalash$y,"n",xlab="x",ylab="y")
text(coalash$x,coalash$y,label=coalash$coalash,cex=0.5)
“n” means no points are drawn in the plot
par(mfrow=c(1,1))
plot(unique(coalash$x),tapply(coalash$coalash,coalash$x,mean),xlab="x",
ylab="mean coal",pch=20)
plot(unique(coalash$y),tapply(coalash$coalash,coalash$y,mean),xlab="y",
ylab="mean coal",pch=20)
2.2.2 - tapply(): computes the mean coalash value for each unique x
2.2.3 - unique(): gives the distinct x values
2.2.4 - plot(): makes a scatter plot of mean coalash vs. x
2.2.5 - pch=20: plots solid circle points
2.2.6 - Overall: the plot shows how the average coalash changes across x
m <- median(coalash$coalash)
plot(coalash$x[coalash$coalash <=m], coalash$y[coalash$coalash <=m], pch =5, xlab = "x", ylab = "y")
points(coalash$x[coalash$coalash >m], coalash$y[coalash$coalash >m], pch =20) 
2.3 Average temperature data
Try for Average temperature data. In which region temperature greater or less than mean…
dat <- read_csv("G:/My Drive/Turin/5th year MS/AST 511 [Enviromental and Spatial Statistics]/practical part/Spatial data/Data/Average_temp_2018.csv")
plot(dat$Lon, dat$Lat, pch=4,cex=0.75,xlab="Longitude",ylab="Latitude")
scatterplot3d(dat$Lon, dat$Lat, dat$Temp, xlab="Longitude",ylab="Latitude",
zlab = "Temperature")
plot(dat$Lon, dat$Lat, "n", pch=4,cex=0.75,xlab="Longitude",ylab="Latitude")
text(dat$Lon, dat$Lat, labels = dat$Temp, cex=0.75)
par(mfrow=c(1,1))
plot(unique(dat$Lon), tapply(dat$Temp, dat$Lon, mean), xlab="Longitude",
ylab="mean Temp",pch=20)
plot(unique(dat$Lat), tapply(dat$Temp, dat$Lat, mean), xlab="Latitude",
ylab="mean Temp",pch=20)
m1<- median(dat$Temp)
plot(dat$Lon[dat$Temp<=m1], dat$Lat[dat$Temp<=m1], pch=5, xlab = "Longitude", ylab = "Latitude")
points(dat$Lon[dat$Temp>m1], dat$Lat[dat$Temp>m1], pch=20)
2.3.1 For Geographical map
# Convert to sf
dat.sf <- st_as_sf(dat, coords = c("Lon", "Lat"), remove = FALSE, crs = 4326)
bd <- st_read("G:/My Drive/Turin/5th year MS/AST 511 [Enviromental and Spatial Statistics]/practical part/Spatial data/Data/BGD_adm0.shp")## Reading layer `BGD_adm0' from data source
## `G:\My Drive\Turin\5th year MS\AST 511 [Enviromental and Spatial Statistics]\practical part\Spatial data\Data\BGD_adm0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 67 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 88.01057 ymin: 20.74111 xmax: 92.67366 ymax: 26.6344
## Geodetic CRS: WGS 84# Transform to proper CRS for Bangladesh (UTM zone 46N)
dat.crs <- st_transform(dat.sf, crs = 32646)
bd.crs <- st_transform(bd, crs = 32646)
# Plot
ggplot() +
geom_sf(data = bd.crs) +
geom_sf(data = dat.crs, aes(color = as.factor(Station)), size = 2) +
labs(x = "Longitude", y = "Latitude", title = "Median size") +
theme_bw()
2.4 VARIOGRAM CLOUD
data(elevation)
head(elevation)## $coords
## x y
## 1 0.3 6.1
## 2 1.4 6.2
## 3 2.4 6.1
## 4 3.6 6.2
## 5 5.7 6.2
## 6 1.6 5.2
## 7 2.9 5.1
## 8 3.4 5.3
## 9 3.4 5.7
## 10 4.8 5.6
## 11 5.3 5.0
## 12 6.2 5.2
## 13 0.2 4.3
## 14 0.9 4.2
## 15 2.3 4.8
## 16 2.5 4.5
## 17 3.0 4.5
## 18 3.5 4.5
## 19 4.1 4.6
## 20 4.9 4.2
## 21 6.3 4.3
## 22 0.9 3.2
## 23 1.7 3.8
## 24 2.4 3.8
## 25 3.7 3.5
## 26 4.5 3.2
## 27 5.2 3.2
## 28 6.3 3.4
## 29 0.3 2.4
## 30 2.0 2.7
## 31 3.8 2.3
## 32 6.3 2.2
## 33 0.6 1.7
## 34 1.5 1.8
## 35 2.1 1.8
## 36 2.1 1.1
## 37 3.1 1.1
## 38 4.5 1.8
## 39 5.5 1.7
## 40 5.7 1.0
## 41 6.2 1.0
## 42 0.4 0.5
## 43 1.4 0.6
## 44 1.4 0.1
## 45 2.1 0.7
## 46 2.3 0.3
## 47 3.1 0.0
## 48 4.1 0.8
## 49 5.4 0.4
## 50 6.0 0.1
## 51 5.7 3.0
## 52 3.6 6.0
##
## $data
## [1] 870 793 755 690 800 800 730 728 710 780 804 855 830 813 762 765 740 765 760
## [20] 790 820 855 812 773 812 827 805 840 890 820 873 875 873 865 841 862 908 855
## [39] 850 882 910 940 915 890 880 870 880 960 890 860 830 705class(elevation)## [1] "geodata"2.4.1 Empirical Variograms
plot(variog(elevation, op = "cloud"), pch = 20)## variog: computing omnidirectional variogram
Spatial autocorrelation is decreasing with the increase of distance.
2.4.2 For Temperature Data
class(dat)## [1] "spec_tbl_df" "tbl_df" "tbl" "data.frame"names(dat)## [1] "...1" "Station" "Year" "Temp" "Lat" "Lon"dat.gdata<-as.geodata(dat, coords.col = c("Lat", "Lon"), data.col = "Temp")
#dat.gdata1<-as.geodata(dat, coords.col = c(5,6), data.col = 4)
plot(variog(dat.gdata, op = "cloud"), pch = 20)## variog: computing omnidirectional variogram
2.5 ESTIMATED SEMI-VARIOGRAM
Two different estimation methods (classical vs. modulus)
# Classical
estvar<-variog(elevation)## variog: computing omnidirectional variogramnames(estvar)## [1] "u" "v" "n" "sd"
## [5] "bins.lim" "ind.bin" "var.mark" "beta.ols"
## [9] "output.type" "max.dist" "estimator.type" "n.data"
## [13] "lambda" "trend" "pairs.min" "nugget.tolerance"
## [17] "direction" "tolerance" "uvec" "call"estvar$v## [1] 168.6250 695.7431 1195.4654 2206.3057 3073.5943 4037.8634 4605.0765
## [8] 6316.2616 6641.9603 6319.5886 5444.7903 3417.1000 3261.1250# Modulus
estvarm<-variog(elevation, estimator.type="modulus")## variog: computing omnidirectional variogramestvarm$v## [1] 216.0084 875.2365 1523.2942 2743.6665 3779.0148 4868.4651 5562.8754
## [8] 7817.4419 7075.8987 6754.4119 6959.8737 3681.7172 3033.4134plot(estvar,pch=20)
plot(estvarm,pch=20)
2.6 Fit different variogram models
EXPERIMENTAL VARIOGRAM FIITING ON EXPERIMENTAL PLOT
plot(estvar,pch=20)
# Gaussian model
fit_g<-variofit(estvar,ini.cov.pars=c(5000,5)
,cov.model="gaussian", fix.nugget = FALSE, nugget = 1)## variofit: covariance model used is gaussian
## variofit: weights used: npairs
## variofit: minimisation function used: optimlines(fit_g)
# Exponential model
fit_e<-variofit(estvar,ini.cov.pars=c(5000,5)
,cov.model="exponential", fix.nugget = FALSE, nugget = 1)## variofit: covariance model used is exponential
## variofit: weights used: npairs
## variofit: minimisation function used: optimlines(fit_e, col = "red")
# Spherical model
fit_s<-variofit(estvar,ini.cov.pars=c(5000,5)
,cov.model="spherical", fix.nugget = FALSE, nugget = 1)## variofit: covariance model used is spherical
## variofit: weights used: npairs
## variofit: minimisation function used: optimlines(fit_s, col = "blue")
# Circular model
fit_c<-variofit(estvar,ini.cov.pars=c(5000,5)
,cov.model="circular", fix.nugget = FALSE, nugget = 1)## variofit: covariance model used is circular
## variofit: weights used: npairs
## variofit: minimisation function used: optimlines(fit_c, col = "green")
# Matérn model
fit_m<-variofit(estvar,ini.cov.pars=c(5000,5)
,cov.model="matern", fix.nugget = FALSE, nugget = 0, fix.kappa = FALSE, kappa = 0.3)## variofit: covariance model used is matern
## variofit: weights used: npairs
## variofit: minimisation function used: optimlines(fit_m, col = "blue4")
2.6.1 For Temperature Data
dat.gdata<-as.geodata(dat, coords.col = c("Lat", "Lon"), data.col = "Temp")
estvar1 <- variog(dat.gdata, estimator.type = "modulus")## variog: computing omnidirectional variogramestvarm1 <- variog(dat.gdata)## variog: computing omnidirectional variogramplot(estvar1, pch=20)
plot(estvarm1, pch=20)
fit_g1<-variofit(estvar1, ini.cov.pars=c(5000,5)
,cov.model="gaussian", fix.nugget = FALSE, nugget = 1)## variofit: covariance model used is gaussian
## variofit: weights used: npairs
## variofit: minimisation function used: optim## Warning in variofit(estvar1, ini.cov.pars = c(5000, 5), cov.model = "gaussian",
## : unreasonable initial value for sigmasq (too high)## Warning in variofit(estvar1, ini.cov.pars = c(5000, 5), cov.model = "gaussian",
## : unreasonable initial value for sigmasq + nugget (too high)lines(fit_g1)