AST 511
Md Monowar Hossain
2025-08-28
1st incourse
library(geoR)## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.9-6 (built on 2025-08-29) is now loaded
## --------------------------------------------------------------library(gstat)
data("coalash")
plot(coalash$x,coalash$y)
plot(coalash$x,coalash$y,pch=4)
library(scatterplot3d)
scatterplot3d(coalash$x,coalash$y,coalash$coalash)
attach(coalash)## The following object is masked _by_ .GlobalEnv:
##
## coalashplot(x,y)
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 R
#setwd("~/AST 511")
setwd("D:/Masters/AST 531 Statistical Computing I/AST 531 (LAB 511)")
#Read in Virginia polygon
virginia <- st_read("Data 511/VirginiaBoundary.shp") ## Reading layer `VirginiaBoundary' from data source
## `D:\Masters\AST 531 Statistical Computing I\AST 531 (LAB 511)\Data 511\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)
#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))
#Examine structure of data frame
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 ...#Convert to sf object
cities.sf <-
st_as_sf(cities, coords = c("Longitude", "Latitude"), remove = FALSE, crs = 4326)
#Plot the cities
plot(cities.sf$geometry)
st_crs(virginia)$epsg## [1] 4269st_crs(cities.sf)$epsg## [1] 4326#Transform CRS
virginia.crs <- st_transform(virginia, crs = 26918)
cities.crs <- st_transform(cities.sf, crs = 26918)
#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) +
#Add titles
labs(x = "Longitude", y = "Latitude", title = "Virginia Cities") +
#Add theme
theme_bw()
#Get bounding box of the cities
cities.extent <- st_bbox(cities.crs)
#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")) +
#Add titles
labs(x = "Longitude", y = "Latitude", title = "Virginia Cities") +
#Add theme
theme_bw()
#Distance from each city to the border of Virginia
dist.border <- st_distance(cities.crs, st_boundary(virginia.crs))
#See results
print(dist.border)## Units: [m]
## [,1]
## [1,] 85469.5113
## [2,] 121.7268
## [3,] 11906.0329
## [4,] 20091.7204
## [5,] 3054.4664
## [6,] 82565.9571#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`#Look at first 3 distances
dist.cville## ...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)
## 4 239121.9 Chesapeake 36.7682 -76.2875 POINT (385093.7 4069931)
## 5 152313.8 Arlington 38.8797 -77.1080 POINT (317147 4305539)
## 6 0.0 Charlottesville 38.0293 -78.4767 POINT (194831 4214774)#Get county data for Virginia from tigris package and transform the CRS
va.counties <- counties(state = "VA") %>%
st_transform(crs = 26918)## Retrieving data for the year 2024## | | | 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% | |========================== | 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%#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% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 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% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 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% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%cities.extent <- st_bbox(cities.crs)
#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) +
#Add titles
labs(x = "Longitude", y = "Latitude", title = "Map of Virginia Cities") +
#Adjust theme
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
library(geoR)
library(gstat)
data("coalash")
############# Large scale component ##############
plot(coalash$x,coalash$y,pch=4,cex=0.75,xlab="x",ylab="y")
library(scatterplot3d)
scatterplot3d(coalash$x,coalash$y, coalash$coalash, xlab = "x", ylab = "y", zlab = "coal")
plot(coalash$x,coalash$y,"n",xlab="x",ylab="y")
text(coalash$x,coalash$y,label=coalash$coalash,cex=0.5)
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)
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)
##################################################
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 705plot(variog(elevation, op = "cloud"), pch = 20)## variog: computing omnidirectional variogram
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.1250estvarm<-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.4134fit_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: optimfit_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: optimfit_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: optimfit_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: optimfit_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: optimplot(estvar,pch=20)
lines(fit_g)
lines(fit_e, col = "red")
lines(fit_s, col = "blue")
lines(fit_c, col = "green")
lines(fit_m, col = "blue4")
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 R
#Read in Bangladesh
bgd <- st_read("Data 511/BGD_adm0.shp") ## Reading layer `BGD_adm0' from data source
## `D:\Masters\AST 531 Statistical Computing I\AST 531 (LAB 511)\Data 511\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#plot(bgd)
as.character(st_geometry_type(bgd))[1]## [1] "MULTIPOLYGON"plot(bgd$geometry)
library(readr)
temp <- read_csv("Data 511/Average_temp_2018.csv")## New names:
## Rows: 35 Columns: 6
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (1): Station dbl (5): ...1, Year, Temp, Lat, Lon
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`#Examine structure of data frame
str(temp)## spc_tbl_ [35 × 6] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ ...1 : num [1:35] 1 2 3 4 5 6 7 8 9 10 ...
## $ Station: chr [1:35] "Ambagan(Ctg)" "Barisal" "Bhola" "Bogra" ...
## $ Year : num [1:35] 2018 2018 2018 2018 2018 ...
## $ Temp : num [1:35] 25.7 25.4 25.6 25.3 26.2 ...
## $ Lat : num [1:35] 22.4 22.8 22.7 24.9 23.3 ...
## $ Lon : num [1:35] 91.8 90.4 90.7 89.4 90.7 ...
## - attr(*, "spec")=
## .. cols(
## .. ...1 = col_double(),
## .. Station = col_character(),
## .. Year = col_double(),
## .. Temp = col_double(),
## .. Lat = col_double(),
## .. Lon = col_double()
## .. )
## - attr(*, "problems")=<externalptr>tem.sf<-st_as_sf(x = temp,coords = c("Lon", "Lat"), remove = FALSE, crs = 4326)
#Plot the cities
plot(tem.sf$geometry)
st_crs(bgd)$epsg## [1] 4326st_crs(tem.sf)$epsg## [1] 4326#Transform CRS
# virginia.crs <- st_transform(virginia, crs = 26918)
# cities.crs <- st_transform(cities.sf, crs = 26918)
#Check that the layers have the same CRS
# st_crs(virginia.crs) == st_crs(cities.crs)
ggplot() +
#Add Virginia polygon
geom_sf(data = bgd) +
#Add cities layer
geom_sf(data = tem.sf, aes(color = Station, label = Station), size = 3) +
#Add titles
labs(x = "Longitude", y = "Latitude", title = "Bangladesh Cities") +
#Add theme
theme_bw()
#Get bounding box of the cities
tem.extent <- st_bbox(tem.sf)
#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(bgd, tem.extent)) +
#Add cities layer
geom_sf(data = tem.sf, aes(color = Station, label = Station), 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")) +
#Add titles
labs(x = "Longitude", y = "Latitude", title = "Virginia Cities") +
#Add theme
theme_bw()
# Distance from each Station to the border of Bangladesh
dist.border <- st_distance(tem.sf, st_boundary(bgd))
# See results
print(dist.border)## Units: [m]
## [,1]
## [1,] 6342.1591
## [2,] 2491.4450
## [3,] 4080.1436
## [4,] 56254.0368
## [5,] 8151.8056
## [6,] 4298.5340
## [7,] 17206.1147
## [8,] 7989.6323
## [9,] 331.0868
## [10,] 65087.2476
## [11,] 16381.1277
## [12,] 86839.0141
## [13,] 2333.5789
## [14,] 2033.8490
## [15,] 32486.1776
## [16,] 17869.9803
## [17,] 532.3991
## [18,] 3361.3528
## [19,] 120.0074
## [20,] 17683.7922
## [21,] 28625.1446
## [22,] 320.5012
## [23,] 47611.9419
## [24,] 2088.5302
## [25,] 6028.5595
## [26,] 34190.6594
## [27,] 33044.4426
## [28,] 259.4593
## [29,] 13045.2760
## [30,] 5370.7863
## [31,] 5968.2587
## [32,] 27715.3041
## [33,] 30938.5015
## [34,] 110094.9736
## [35,] 875.9095# Calculate distance from Dhaka to each Station
dist.cville <-
st_distance(tem.sf, filter(tem.sf, Station == "Dhaka")) %>%
bind_cols(tem.sf)## New names:
## • `` -> `...1`
## • `...1` -> `...2`# look all
dist.cville## ...1 ...2 Station Year Temp Lat Lon
## 1 215239.00 1 Ambagan(Ctg) 2018 25.68333 22.35000 91.81700
## 2 113061.17 2 Barisal 2018 25.36667 22.75000 90.36667
## 3 123504.73 3 Bhola 2018 25.60833 22.68333 90.65000
## 4 158507.77 4 Bogra 2018 25.29167 24.85000 89.36667
## 5 64292.63 5 Chandpur 2018 26.20833 23.26667 90.70000
## 6 222063.76 6 Chittagong 2018 26.05833 22.27000 91.82000
## 7 160029.26 7 chuadanga 2018 24.80000 23.65000 88.81667
## 8 89546.83 8 Comilla 2018 25.17500 23.43333 91.18333
## 9 304579.90 9 Cox's Bazar 2018 25.89167 21.45000 91.96667
## 10 0.00 10 Dhaka 2018 26.08333 23.76667 90.38333
## 11 270814.57 11 Dinajpur 2018 24.64167 25.65000 88.68333
## 12 57384.34 12 Faridpur 2018 25.49167 23.60000 89.85000
## 13 133300.42 13 Feni 2018 25.05000 23.03333 91.41667
## 14 165389.34 14 Hatiya 2018 25.53333 22.43333 91.10000
## 15 141495.45 15 Ishurdi 2018 24.86667 24.13333 89.05000
## 16 140019.23 16 Jessore 2018 25.65833 23.18333 89.16667
## 17 198892.44 17 Khepupara 2018 25.77500 21.98333 90.23333
## 18 139620.00 18 Khulna 2018 26.15000 22.78333 89.53333
## 19 263851.73 19 Kutubdia 2018 25.85833 21.81667 91.85000
## 20 123977.81 20 M.court 2018 25.87500 22.86667 91.10000
## 21 69766.09 21 Madaripur 2018 25.40833 23.16667 90.18333
## 22 165265.79 22 Mongla 2018 26.02500 22.46667 89.60000
## 23 105756.92 23 Mymensingh 2018 24.72500 24.71667 90.43333
## 24 159461.72 24 Patuakhali 2018 25.96667 22.33333 90.33333
## 25 183465.78 25 Rajshahi 2018 25.00833 24.36667 88.70000
## 26 230876.46 26 Rangamati 2018 24.70833 22.53333 92.20000
## 27 247600.64 27 Rangpur 2018 24.71667 25.73333 89.23333
## 28 178582.47 28 Sandwip 2018 25.70000 22.48333 91.43333
## 29 176840.99 29 Satkhira 2018 25.97500 22.71667 89.08333
## 30 188222.34 30 Sitakunda 2018 25.03333 22.58333 91.70000
## 31 149375.30 31 Srimangal 2018 24.14167 24.30000 91.73333
## 32 265640.86 32 sydpur 2018 24.90000 25.75000 88.91667
## 33 197424.44 33 Sylhet 2018 24.80000 24.90000 91.88333
## 34 70553.00 34 Tangail 2018 25.16667 24.25000 89.93333
## 35 377946.50 35 Teknaf 2018 26.08333 20.86667 92.30000
## geometry
## 1 POINT (91.817 22.35)
## 2 POINT (90.36667 22.75)
## 3 POINT (90.65 22.68333)
## 4 POINT (89.36667 24.85)
## 5 POINT (90.7 23.26667)
## 6 POINT (91.82 22.27)
## 7 POINT (88.81667 23.65)
## 8 POINT (91.18333 23.43333)
## 9 POINT (91.96667 21.45)
## 10 POINT (90.38333 23.76667)
## 11 POINT (88.68333 25.65)
## 12 POINT (89.85 23.6)
## 13 POINT (91.41667 23.03333)
## 14 POINT (91.1 22.43333)
## 15 POINT (89.05 24.13333)
## 16 POINT (89.16667 23.18333)
## 17 POINT (90.23333 21.98333)
## 18 POINT (89.53333 22.78333)
## 19 POINT (91.85 21.81667)
## 20 POINT (91.1 22.86667)
## 21 POINT (90.18333 23.16667)
## 22 POINT (89.6 22.46667)
## 23 POINT (90.43333 24.71667)
## 24 POINT (90.33333 22.33333)
## 25 POINT (88.7 24.36667)
## 26 POINT (92.2 22.53333)
## 27 POINT (89.23333 25.73333)
## 28 POINT (91.43333 22.48333)
## 29 POINT (89.08333 22.71667)
## 30 POINT (91.7 22.58333)
## 31 POINT (91.73333 24.3)
## 32 POINT (88.91667 25.75)
## 33 POINT (91.88333 24.9)
## 34 POINT (89.93333 24.25)
## 35 POINT (92.3 20.86667)library(geoR)
library(gstat)
library(readr)
tem <- read_csv("Data 511/Average_temp_2018.csv")## New names:
## Rows: 35 Columns: 6
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (1): Station dbl (5): ...1, Year, Temp, Lat, Lon
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`############# Large scale component ##############
plot(tem$Lat,tem$Lon,pch=4,cex=0.75,xlab="x",ylab="y")
library(scatterplot3d)
scatterplot3d(tem$Lat,tem$Lon, tem$Tem, xlab = "x", ylab = "y", zlab = "coal")
plot(tem$Lat,tem$Lon,"n")
text(tem$Lat,tem$Lon,label=round(tem$Temp,1),cex=0.5)
m <- mean(tem$Temp)
plot(tem$Lat[tem$Temp <=m], tem$Lon[tem$Temp <=m], pch =5, xlab = "x", ylab = "y")
points(tem$Lat[tem$Temp >m], tem$Lon[tem$Temp >m], pch =20) 
tem.sf<-st_as_sf(x = temp,coords = c("Lon", "Lat"), remove = FALSE, crs = 4326)
#Plot the cities
plot(tem.sf$geometry)
tem.sf<-tem.sf %>%
mutate(ind= ifelse(Temp> m,1,0))
ggplot() +
#Add Virginia polygon
geom_sf(data = bgd) +
#Add cities layer
geom_sf(data = tem.sf,aes(colour=ind,label=ind)) +
#Add titles
labs(x = "Longitude", y = "Latitude", title = "Bangladesh Cities") +
#Add theme
theme_bw()
temp.geo<-as.geodata(obj = tem,coords.col = c("Lat","Lon"),data.col = "Temp")
plot(variog(temp.geo, op = "cloud"), pch = 20)## variog: computing omnidirectional variogram
estvar<-variog(temp.geo)## 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] 0.1101736 0.1566892 0.1736328 0.2630166 0.3666166 0.3439492 0.4121331
## [8] 0.3871172 0.3632668 0.4757031 0.4574074 0.6544213 0.8170139estvarm<-variog(temp.geo,estimator.type="modulus")## variog: computing omnidirectional variogramestvarm$v## [1] 0.09313272 0.17600819 0.20014569 0.23310726 0.42897885 0.42226529
## [7] 0.40777491 0.44315910 0.56044560 0.59703263 0.57959215 1.03836586
## [13] 1.15157734plot(estvar,pch=20)
plot(estvarm,pch=20)
### for variogram function need sf
v <- variogram(Temp ~ 1, tem.sf, cloud = TRUE)
plot(v)
after 1st incourse
Grid
library("sp")
library("sf")
library("stars") ## Loading required package: abinddata("meuse")
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
bbox <- st_bbox(meuse)
bbox## xmin ymin xmax ymax
## 178605 329714 181390 333611cell_size <- 40
x <- seq(bbox$xmin, bbox$xmax, by=cell_size)
y <- seq(bbox$ymin, bbox$ymax, by=cell_size)
meuse_grid <- expand.grid(x=x, y=y)
plot(meuse_grid$x, meuse_grid$y, cex=0.1)
meuse_grid$tmp <- 1
meuse_grid <- st_as_stars(meuse_grid, crs=st_crs(meuse))
st_crs(meuse_grid) <- st_crs(meuse) Kriging
library(sp)
library(gstat)
library(sf)
library(mapview)
data("meuse")
data(meuse.grid)
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
plot(meuse)
plot(meuse$geometry)
mapview(meuse, zcol = "zinc", map.types = "CartoDB.Voyager")meuse.grid <- st_as_sf(meuse.grid, coords = c("x", "y"),
crs = 28992)
mapview(meuse.grid, zcol="dist", map.types = "CartoDB.Voyager")ggplot()+
geom_sf(data=meuse)
ggplot()+
geom_sf(data=meuse.grid)
hist(meuse$zinc)
hist(log(meuse$zinc))
measure the variogram
vc <- variogram(log(zinc) ~ 1, meuse, cloud = TRUE)
plot(vc)
v <- variogram(log(zinc) ~ 1, data = meuse)
plot(v)
vgm()## 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)fit the variogram model
vinitial <- vgm(psill = 0.5, model = "Sph",
range = 900, nugget = 0.1)
plot(v, vinitial, cutoff = 1000, cex = 1.5)
fv <- fit.variogram(object = v,
model = vinitial)
fv## model psill range
## 1 Nug 0.05066017 0.0000
## 2 Sph 0.59060556 897.0044plot(v, fv, cex = 1.5)
krining/ prediction
library(ggplot2)
library(viridis)## Loading required package: viridisLitek <- gstat(formula = log(zinc) ~ 1, data = meuse, model = fv)
kpred <- predict(k, meuse.grid)## [using ordinary kriging]ggplot() +
geom_sf(data = kpred, aes(color = var1.pred))+
geom_sf(data = meuse) +
scale_color_viridis(name = "log(zinc)") +
theme_bw()
ggplot() +
geom_sf(data = kpred, aes(color = var1.var)) +
geom_sf(data = meuse) +
scale_color_viridis(name = "variance") +
theme_bw()
temp_data
library(readr)
data <- read_csv("Data 511/Average_temp_2018.csv")## New names:
## Rows: 35 Columns: 6
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (1): Station dbl (5): ...1, Year, Temp, Lat, Lon
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`head(data)## # A tibble: 6 × 6
## ...1 Station Year Temp Lat Lon
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1 Ambagan(Ctg) 2018 25.7 22.4 91.8
## 2 2 Barisal 2018 25.4 22.8 90.4
## 3 3 Bhola 2018 25.6 22.7 90.6
## 4 4 Bogra 2018 25.3 24.8 89.4
## 5 5 Chandpur 2018 26.2 23.3 90.7
## 6 6 Chittagong 2018 26.1 22.3 91.8data2<- st_as_sf(data,coords = c("Lon","Lat"),crs= 28992)
bbox<-st_bbox(data2)
bbox## xmin ymin xmax ymax
## 88.68333 20.86667 92.30000 25.75000cell_si<-40
#x<-seq(bbox$xmin,bbox$ymin, by =40)
#x
grid<-st_make_grid(x = data2,cellsize = 0.01,what = "centers")
grid<- st_sf(grid)
grid<- st_sf(geometry=grid)
grid## Simple feature collection with 177018 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 88.68833 ymin: 20.87167 xmax: 92.29833 ymax: 25.75167
## Projected CRS: Amersfoort / RD New
## First 10 features:
## grid
## 1 POINT (88.68833 20.87167)
## 2 POINT (88.69833 20.87167)
## 3 POINT (88.70833 20.87167)
## 4 POINT (88.71833 20.87167)
## 5 POINT (88.72833 20.87167)
## 6 POINT (88.73833 20.87167)
## 7 POINT (88.74833 20.87167)
## 8 POINT (88.75833 20.87167)
## 9 POINT (88.76833 20.87167)
## 10 POINT (88.77833 20.87167)library("gstat")
vc<-variogram(object = Temp~1,data = data2,cloud=T)
plot(vc)
v<-variogram(object = Temp~1,data = data2)
plot(v)
# Sphericalmodel
v_sph<-fit.variogram(v,model= vgm(psill=0.5,model="Sph",range=900,nugget=0.09))## Warning in fit.variogram(v, model = vgm(psill = 0.5, model = "Sph", range =
## 900, : No convergence after 200 iterations: try different initial values?# Exponentialmodel
v_exp<-fit.variogram(v,model= vgm(psill=0.4,model="Exp",range=850,nugget=0.08))## Warning in fit.variogram(v, model = vgm(psill = 0.4, model = "Exp", range =
## 850, : No convergence after 200 iterations: try different initial values?# Gaussianmodel
v_gau<-fit.variogram(v,model= vgm(psill=0.45,model="Gau",range=950,nugget=0.10))
# Maternmodel
v_mat<-fit.variogram(v,model= vgm(psill=0.48,model="Mat",range=900,nugget=0.09,kappa=0.5))## Warning in fit.variogram(v, model = vgm(psill = 0.48, model = "Mat", range =
## 900, : No convergence after 200 iterations: try different initial values?# Linearmodel
v_lin<-fit.variogram(v,model= vgm(psill=0.48,model="Lin",range=900,nugget=0.09,kappa=0.5))## Warning in fit.variogram(v, model = vgm(psill = 0.48, model = "Lin", range =
## 900, : singular model in variogram fitplot(v,v_sph,main= 'Sphericalmodel')
plot(v,v_exp,main= 'Exponentialmodel')
plot(v,v_gau,main= 'Gaussianmodel')
plot(v,v_mat,main= 'Maternmodel')
plot(v,v_lin,main= 'Linearmodel')
kriging/ prediction
library(viridis)
k<-gstat(formula = Temp ~1,data = data2,model = v_gau)
kpre<- predict(k,grid)## [using ordinary kriging]ggplot()+
geom_sf(data = kpre,aes(color=var1.pred))+
geom_sf(data = data2)+
scale_color_viridis(name="Temp")
ggplot()+
geom_sf(data = kpre,aes(color=var1.var))+
geom_sf(data = data2)+
scale_color_viridis()