library(geoR)
## Warning: package 'geoR' was built under R version 4.4.3
## --------------------------------------------------------------
## 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)
## Warning: package 'gstat' was built under R version 4.4.3
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:
##
## coalash
plot(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 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 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: NAD83
plot(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] 4269
st_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] TRUE
ggplot() +
#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()
## Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
## : Ignoring unknown aesthetics: label
#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()
## 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`
#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
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)
#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% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 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%
#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))
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
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 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] 4326
st_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 at first 3 distances
head(dist.cville, n = 3)
## ...1 ...2 Station Year Temp Lat Lon
## 1 215239.0 1 Ambagan(Ctg) 2018 25.68333 22.35000 91.81700
## 2 113061.2 2 Barisal 2018 25.36667 22.75000 90.36667
## 3 123504.7 3 Bhola 2018 25.60833 22.68333 90.65000
## geometry
## 1 POINT (91.817 22.35)
## 2 POINT (90.36667 22.75)
## 3 POINT (90.65 22.68333)
# 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)
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 705
plot(variog(elevation, op = "cloud"), pch = 20)
## variog: computing omnidirectional variogram
estvar<-variog(elevation)
## variog: computing omnidirectional variogram
names(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
estvarm<-variog(elevation,estimator.type="modulus")
## variog: computing omnidirectional variogram
estvarm$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.4134
plot(estvar,pch=20)
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: optim
lines(fit_g)
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: optim
lines(fit_e, col = "red")
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: optim
lines(fit_s, col = "blue")
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: optim
lines(fit_c, col = "green")
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: optim
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(g.l$u, cov.model = g.l$cov.model, kappa = kappa, : NaN
## value in cov.spatial
lines(fit_m, col = "blue4")
## Warning in cov.spatial(x, cov.model = my.l$cov.model, kappa = my.l$kappa, :
## Infinity value in cov.spatial
## Warning in cov.spatial(x, cov.model = my.l$cov.model, kappa = my.l$kappa, : NA
## value in cov.spatial
## Warning in cov.spatial(x, cov.model = my.l$cov.model, kappa = my.l$kappa, : NaN
## value in cov.spatial
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")
## Warning: Unknown or uninitialised column: `Tem`.
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()
## Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
## : Ignoring unknown aesthetics: label
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 variogram
names(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.8170139
estvarm<-variog(temp.geo,estimator.type="modulus")
## variog: computing omnidirectional variogram
estvarm$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.15157734
plot(estvar,pch=20)
plot(estvarm,pch=20)
