Github-R-intro/Spatial test.Rmd at main · arsenyerlan/Github-R-intro · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
---
title: "Spatial test"
output: html_document
date: "2023-01-20"
---

```{r}
source("setup.R")
#Import colorado counties with tigris

counties<-counties(state="CO")
#Import roads for Larimer county
roads<-roads(state="CO", county="Larimer")

# set tmap mode to interactive
tmap_mode("view")


qtm(counties)+
  qtm(roads)

tm_shape(counties)+
  tm_polygons()
#Look at the class of counties
class(counties)
```

```{r}
# point data
poudre_points <- data.frame(name = c("Mishawaka", "Rustic", "Blue Lake Trailhead"),
                            long = c(-105.35634, -105.58159, -105.85563),
                            lat = c(40.68752, 40.69687, 40.57960))
#convert to spatial
poudre_points_sf<-st_as_sf(poudre_points, coords=c("long", "lat"), crs= 4326)
class(poudre_points_sf)
poudre_points_sf
```


```{r}
#raster data
elevation<-get_elev_raster(counties, z= 7)
elevation
qtm(elevation)
tm_shape(elevation)+
  tm_raster(style="cont", title="Elevation (m)")
```


```{r}
#the terra package
elevation<-rast(elevation)
elevation
names(elevation)<-"Elevation"
elevation
#check projections
st_crs(counties)

crs(counties)==crs(elevation)

#project elevation layer
elevation_prj<-terra::project(elevation, counties)
#crop elevation to counties extent
elevation_crop<-crop(elevation, ext(counties))
qtm(elevation_crop)
```


```{r}
#Exercises
##1 Filter out the counties data set to only include Larimer, Denver, and Pueblo counties.

counties_fl<-counties %>%
  filter(NAME %in% c('Larimer', 'Denver', 'Pueblo'))
qtm(counties_fl)
##2 Make a map of the counties data colored by county area. Make a second map of counties colored by their total area of water.
tm_shape(counties_fl) +
  tm_polygons(col = 'ALAND', )
tm_shape(counties_fl) +
  tm_polygons(col = 'AWATER')
?tm_polygons
## Make a barplot comparing the elevation of your 3 points in the Poudre Canyon


roads %>%
  filter(grepl('Poudre',FULLNAME, ignore.case = T))

roads %>%
  filter(str_detect(FULLNAME, "Poudre"))

### The code above contains use code for finding information in a vector that you do not know the exact name of.
pdrHwy <- roads %>%
  filter(FULLNAME == 'Poudre Canyon Hwy')

elevation <- get_elev_raster(counties, z = 7)
elevation <- rast(elevation)
names(elevation) <- 'Elevation'
elevationCrop <- crop(elevation, ext(pdrHwy))
writeRaster(elevationCrop, 'data/poudreElevation.tif')
pdrPointElevation <- extract(x= elevationCrop,
                             y= pdrPntsProject,
                             xy= TRUE
                             )
pdrPointElevation %>%
  ggplot() +
    geom_col(mapping = aes(x=ID, y=Elevation))

?extract()
```


?extract()

#read and write spatial data

#safe sf/vector data
write_sf(counties, "data/counties.shp")

#save raster data
writeRaster(elevation_crop,"data/elevation.tif")

#save R data
save(counties, roads, file= "data/spatial_objects.RData")