M9_Lab8_Task2_Harold
Samuel Harold
2022-03-16
setwd("C:/Users/Sam/Documents/DSCI_605/Module_9/")
library(geoR)
library(rgeos)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)# help(parana)
names(parana)## $coords
## [1] "east" "north"
##
## $data
## [1] "data"
##
## $borders
## [1] "borders"
##
## $other
## [1] "loci.paper"parana$coords## east north
## [1,] 402.9529 164.52841
## [2,] 501.7049 428.77100
## [3,] 556.3262 445.27065
## [4,] 573.4043 447.04177
## [5,] 702.4228 272.29590
## [6,] 668.5442 261.67070
## [7,] 435.8477 286.54044
## [8,] 434.0125 317.90506
## [9,] 432.4622 288.37001
## [10,] 162.1458 286.29984
## [11,] 169.4766 334.53158
## [12,] 717.8468 181.56741
## [13,] 713.3433 214.88336
## [14,] 726.6630 207.27952
## [15,] 618.4688 127.43265
## [16,] 661.8808 128.81455
## [17,] 506.7167 212.86691
## [18,] 561.7874 138.89773
## [19,] 585.6395 212.59809
## [20,] 531.7106 131.62153
## [21,] 308.9718 181.97682
## [22,] 306.1047 148.69468
## [23,] 150.1220 162.12624
## [24,] 155.3399 154.86652
## [25,] 620.0018 112.64927
## [26,] 475.0577 83.64484
## [27,] 366.8719 92.21224
## [28,] 479.4953 461.96806
## [29,] 615.7235 400.63795
## [30,] 599.0729 456.12037
## [31,] 566.5676 445.22857
## [32,] 484.6613 423.22788
## [33,] 445.5438 400.99330
## [34,] 498.2953 426.92592
## [35,] 477.9447 356.79285
## [36,] 353.1834 446.50166
## [37,] 338.5246 381.74695
## [38,] 267.2320 373.39454
## [39,] 676.8785 254.18138
## [40,] 538.8900 308.77851
## [41,] 340.9647 271.01967
## [42,] 203.9369 309.41998
## [43,] 254.2625 240.19632
## [44,] 719.4624 177.84707
## [45,] 768.5087 199.12352
## [46,] 721.6547 209.20941
## [47,] 687.7536 187.56295
## [48,] 599.0604 210.66083
## [49,] 541.8933 185.11924
## [50,] 449.6865 196.16437
## [51,] 357.5152 188.11980
## [52,] 222.3395 154.59123
## [53,] 297.3230 176.26524
## [54,] 184.1727 211.03701
## [55,] 401.9894 70.37282
## [56,] 365.3621 77.42565
## [57,] 331.6807 110.26775
## [58,] 294.9501 113.43465
## [59,] 674.3159 185.89772
## [60,] 518.6908 378.94032
## [61,] 562.7977 362.20880
## [62,] 399.9123 354.61628
## [63,] 309.6951 375.85998
## [64,] 308.2159 359.22457
## [65,] 397.5989 443.18080
## [66,] 206.2979 359.35560
## [67,] 289.3800 370.04023
## [68,] 246.8030 374.89923
## [69,] 279.5530 345.88871
## [70,] 231.8296 450.35979
## [71,] 245.5946 445.05926
## [72,] 410.8090 242.10715
## [73,] 400.0153 339.85273
## [74,] 405.1119 338.04194
## [75,] 422.4753 262.48003
## [76,] 457.7338 306.92340
## [77,] 388.4718 299.16275
## [78,] 299.0782 285.24334
## [79,] 320.4695 329.84407
## [80,] 317.7274 279.96009
## [81,] 328.1532 256.09127
## [82,] 341.9647 271.01967
## [83,] 346.0594 313.53010
## [84,] 371.7581 278.71589
## [85,] 377.0997 247.38436
## [86,] 395.4668 267.83959
## [87,] 244.7930 297.28745
## [88,] 232.0861 343.23287
## [89,] 212.3630 311.44174
## [90,] 250.9284 238.28897
## [91,] 249.1147 339.84732
## [92,] 262.9887 321.61570
## [93,] 264.8356 312.41223
## [94,] 277.1810 281.21777
## [95,] 267.4829 255.20257
## [96,] 287.4249 273.99100
## [97,] 291.2186 246.34864
## [98,] 296.5940 340.60144
## [99,] 295.1108 325.80878
## [100,] 211.6614 263.38789
## [101,] 189.3152 281.38550
## [102,] 165.9700 267.90686
## [103,] 164.2818 267.86614
## [104,] 184.9871 323.79426
## [105,] 711.8998 124.40728
## [106,] 749.5943 175.47115
## [107,] 702.7330 179.95796
## [108,] 650.7635 178.79807
## [109,] 640.7307 180.75360
## [110,] 600.6570 199.57391
## [111,] 421.2830 177.57037
## [112,] 451.4439 174.02255
## [113,] 419.3326 223.70296
## [114,] 377.7171 179.08881
## [115,] 357.4957 189.96577
## [116,] 372.2726 223.34289
## [117,] 287.3583 170.56964
## [118,] 285.6835 170.54282
## [119,] 224.3215 139.84991
## [120,] 218.9926 154.52067
## [121,] 262.6645 146.13850
## [122,] 217.2403 158.17924
## [123,] 212.6301 217.21353
## [124,] 249.3585 142.19763
## [125,] 245.9086 147.67415
## [126,] 251.2398 131.14844
## [127,] 197.7439 205.80201
## [128,] 172.8530 192.28300
## [129,] 187.8333 198.18379
## [130,] 568.3420 114.86877
## [131,] 541.6892 120.51877
## [132,] 493.3411 98.43402
## [133,] 471.6876 103.94229
## [134,] 445.0324 105.70264
## [135,] 479.9921 118.72419
## [136,] 333.2294 119.52053
## [137,] 400.3281 70.36000
## [138,] 343.6607 84.56600
## [139,] 320.2432 93.49421
## [140,] 331.8482 97.34366
## [141,] 235.3269 92.03407
## [142,] 226.5333 114.02601
## [143,] 692.5446 170.87504ggplot(data.frame(cbind(parana$coords, Rainfall = parana$data))) +
geom_point(aes(east, north, color = Rainfall), size = 2) +
coord_fixed(ratio = 1) +
scale_color_gradient(low = "blue", high = "orange") +
geom_path(data = data.frame(parana$border), aes(east, north)) +
theme_bw()
world <- ne_countries(scale = "medium", returnclass = "sf") # get natural world count
class(world)## [1] "sf" "data.frame"names(world)## [1] "scalerank" "featurecla" "labelrank" "sovereignt" "sov_a3"
## [6] "adm0_dif" "level" "type" "admin" "adm0_a3"
## [11] "geou_dif" "geounit" "gu_a3" "su_dif" "subunit"
## [16] "su_a3" "brk_diff" "name" "name_long" "brk_a3"
## [21] "brk_name" "brk_group" "abbrev" "postal" "formal_en"
## [26] "formal_fr" "note_adm0" "note_brk" "name_sort" "name_alt"
## [31] "mapcolor7" "mapcolor8" "mapcolor9" "mapcolor13" "pop_est"
## [36] "gdp_md_est" "pop_year" "lastcensus" "gdp_year" "economy"
## [41] "income_grp" "wikipedia" "fips_10" "iso_a2" "iso_a3"
## [46] "iso_n3" "un_a3" "wb_a2" "wb_a3" "woe_id"
## [51] "adm0_a3_is" "adm0_a3_us" "adm0_a3_un" "adm0_a3_wb" "continent"
## [56] "region_un" "subregion" "region_wb" "name_len" "long_len"
## [61] "abbrev_len" "tiny" "homepart" "geometry"ggplot(data = world) +
# geom_sf() +
geom_sf(color = "black", fill = "lightgreen") +
xlab("Longitude") + ylab("Latitude") +
ggtitle("World map", subtitle = paste0("(", length(unique(world$name_sort)), " countries)"))
ggplot(data = world) +
geom_sf(aes(fill = pop_est)) +
scale_fill_viridis_c(option = "plasma", trans="sqrt")
ggplot(data = world) +
geom_sf() +
coord_sf(crs = "+proj=laea + lat_0=52 +lon_0=-10 +x_0=432000 +y_0=3210000 +ellps=GRS80 + units=m +no_defs ")
ggplot(data = world) +
geom_sf() +
coord_sf(crs = "+init=epsg:3035")
ggplot(data=world) +
geom_sf() +
coord_sf(crs = st_crs(3035)) #st_crs, get vector crs, 3035 in EPSG code
ggplot(data = world) +
geom_sf() +
coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97), expand = FALSE)
library(ggspatial)ggplot(data = world) +
geom_sf() +
annotation_scale(location = "bl", width_hint = 0.5) +
annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97))## Scale on map varies by more than 10%, scale bar may be inaccurate
## Scale on map varies by more than 10% scale bar may be inaccurate
library(sf)
sf::sf_use_s2(FALSE)## Spherical geometry (s2) switched off# sf::sf_use_s2(TRUE)world_points <- st_centroid(world)
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))ggplot(data = world) +
geom_sf() +
geom_text(data = world_points, aes(x = X, y = Y, label = name),
color = "darkblue", fontface = "bold", check_overlap = FALSE) +
annotate(geom = "text", x = -90, y = 26, label = "Gulf of Mexico",
fontface = "italic", color = "grey22", size = 6) +
coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97), expand = FALSE)
ggsave("map.pdf")
ggsave("map_web.png", width = 6, height = 6, dpi = "screen")Spatial Visualization
# devtools:: install_github("robinlovelace/geocompr")
library(raster)
library(dplyr)
# install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source')
library(spData) # gives hints about data
library(tmap) # for static and interactive maps
library(leaflet) # interactive maps
library(devtools)tm_shape(nz) +
tm_fill()
Add border layer to nz shape, will display in a new plot
tm_shape(nz) +
tm_borders()
Add fill and border layers to nz shape
tm_shape(nz) +
tm_fill() +
tm_borders()
map_nz = tm_shape(nz) + tm_polygons() # if you cannot see the plots, you can close R or rm(list-)
class(map_nz)## [1] "tmap"library(spDataLarge)data("nz_elev")
map_nz1 = map_nz +
tm_shape(nz_elev) + tm_raster(alpha = 0.9)nz_water = st_union(nz) %>% st_buffer(22200) %>%
st_cast(to = "LINESTRING")map_nz2 = map_nz1 +
tm_shape(nz_water) + tm_lines()map_nz3 = map_nz2 +
tm_shape(nz_height) + tm_dots()tmap_arrange(map_nz1, map_nz2, map_nz3)
ma1 = tm_shape(nz) + tm_fill(col = 'red')
ma2 = tm_shape(nz) + tm_fill(col = "red", alpha = 0.5)
ma3 = tm_shape(nz) + tm_borders(col = "blue")
ma4 = tm_shape(nz) + tm_borders(lwd = 3)
ma5 = tm_shape(nz) + tm_borders(lty = 3)
ma6 = tm_shape(nz) + tm_fill(col = "red", alpha = 0.5) +
tm_borders(col = "blue", lwd = 3, lty = 2)
tmap_arrange(ma1, ma2, ma3, ma4, ma5, ma6)
plot(st_geometry(nz), col = nz$Land_area) # works
# tm_shape(nz) + tm_fill(col = nz$Land_area) # fails
#. Error: Fill argument neither colors nor valid variable name(s)tm_shape(nz) + tm_fill(col = "Land_area")
legend_title = expression("Area (km"^2*")")
map_nza = tm_shape(nz) +
tm_fill(col = "Land_area", title = legend_title) + tm_borders()tm_shape(nz) + tm_polygons(col = "Median_income")
breaks = c(0, 3, 4, 5) * 10000
tm_shape(nz) + tm_polygons(col = "Median_income", breaks = breaks)
tm_shape(nz) + tm_polygons(col = "Median_income", n = 10)
tm_shape(nz) + tm_polygons(col = "Median_income", palette = "BuGn")
map_nz +
tm_compass(type = "8star", position = c("left", "top")) +
tm_scale_bar(breaks = c(0, 100, 200), text.size = 1)
map_nz + tm_layout(title = "New Zealand")
map_nz + tm_layout(scale = 5)
map_nz + tm_layout(bg.color = "lightblue")
map_nz + tm_layout(frame = FALSE)
map_nz +
tm_compass(type = "8star", position = c("left", "top")) +
tm_scale_bar(breaks = c(0, 100, 200), text.size = 1) +
tm_layout(title = "New Zealand") +
tm_layout(bg.color = "lightblue") +
tm_layout(frame = FALSE)
urb_1970_2030 = urban_agglomerations %>%
filter(year %in% c(1970, 1990, 2010, 2030))tmap_options(check.and.fix = TRUE)tm_shape(world) +
tm_polygons() +
tm_shape(urb_1970_2030) +
tm_symbols(col = "black", border.col = "white", size = "population_millions") +
tm_facets(by = "year", nrow = 2, free.coords = FALSE)
Homework
Please use the data I provide for mapping the points and raster image in one plot area. Download the sample data Spatial point data in .csv file:
plot.locations_HARV <- read.csv("HARV_PlotLocations.csv")Use ggplot(), geom_point(), geom_raster() to plot and overlay the points and raster image.
library(rgdal) # for vector work; sp package should always load with rgdal
library (raster) # for metadata/attributes- vectors or rasters
df <- read.csv("HARV_PlotLocations.csv")
utm18nCRS <- CRS("+proj=utm +zone=18 +datum=WGS84")
dfSp_HARV <- SpatialPointsDataFrame(df[,1:2], # easting and northing position
df, #the R object to convert
proj4string = utm18nCRS) # assign a CRSplot(dfSp_HARV, main="Map of Plot Locations")
library(ggsn);
library(ggplot2);
library(sf)Read points and plot using ggplot
df2 <- read.csv('HARV_PlotLocations.csv')
df2 <-data.frame(cbind(easting = df2$easting, northing = df2$northing, Elevation = df2$elevation))ggplot(df2)+
geom_point(aes(easting, northing, color = Elevation), size = 2) +
coord_fixed(ratio = 1) +
scale_color_gradient(low = "blue", high = "orange") +
theme_bw()
Hint: Install the “raster” package and use raster() function to read the raster image.
Hint: Read the data and transfer into data frame
CHM_HARV <- raster("HARV_chmCrop.tif")
CHM_HARV_df <- as.data.frame(CHM_HARV, xy = TRUE)
names(CHM_HARV_df)[1] <- 'Easting'
names(CHM_HARV_df)[2] <- 'Northing'
names(CHM_HARV_df)[3] <- 'Canopy_Height_Model'Overlay the two data frames. note: df is the dataframe of point data
Try to make your graph more professional since this course is about data visualization.
ggplot(df2)+
geom_raster(data = CHM_HARV_df, aes(x = Easting, y = Northing, fill = Canopy_Height_Model)) +
scale_fill_gradient2(low = "green", high = "dark green") +
geom_point(aes(easting, northing, color = Elevation), size = 3) +
scale_color_gradient(low = "blue", high = "red") +
ggtitle('Northeast Worcester Elevation/Canopy Height Model')
I would add a little more, but I am still unsure of what the geom points are indicating. I was able to find online however, the CHM means Canopy Height Model, so I made this green to match that of trees.