Viewshed Maya

Load in your datasets. This should all be linked ot files on the GIS Server Can take a little while

#load in the DEM and resmaple to 3m resolution
dem3 <- rast('//gis/gisprojects/2024_spring/maya/psx_1_2_dem_1m.tif') %>% 
          terra::aggregate(fact=3,fun=mean)

#load land use DEM and resample to extent of the dem
lu <- rast('//gis/gisprojects/2024_spring/charlie_stuff/land_cover_maya.tif')
ext(lu)<-ext(dem3)
terra::crs(lu)<-terra::crs(dem3)
lu<-terra::resample(lu,dem3)

#load square structures
squares <- st_read('//gis/gisprojects/2024_spring/charlie_stuff/psx_1_2_squares.shp') %>% 
  st_transform(terra::crs(dem3)) %>% 
  mutate(id = row_number())

#load sites
sites <- st_read('//gis/gisprojects/2024_spring/charlie_stuff/psx_1_2_sites.shp') %>% 
  st_transform(terra::crs(dem3))

#create grids at the extent of lu,
grids<-st_make_grid(lu,
                    cellsize=125,
                    # n=c(175,175),
                    square=F) %>% 
                    st_as_sf()

---

Run the whole chunk below to build the function. Don’t modify anything in here! **********

```r
#the following is a function to add random field to hexes, remove buffers around square platforms, add heights to a given percentage of forest hexagons and to all bajo areas

randFor <- function (struct,perc,obs) {

  #add random column to grids and convert to vector. set.seed means you'll get a different result every time
  set.seed( as.integer((as.double(Sys.time())*1000+Sys.getpid()) %% 2^31) )
  grids<-grids %>% mutate(rand=round(runif(nrow(.))*100,0)) %>% 
                      vect()
  
  #rasterize the grids
  forPercRast <- rasterize(x=grids,y=dem3,field=\rand\,fun=max)
  
  #rasterize buffers around square platforms
  gridsBuff <- rasterize(x=squares %>% st_buffer(200),y=dem3,fun=max)
  
  #get rid of non forested pixels
  forPercRast[lu!=7]<-NA
  
  #get rid of areas outside of the buffer
  forPercRast[gridsBuff==1]<-NA
  remove(gridsBuff)
  
  forPercRast[is.na(forPercRast)]<-0 # set non forest to 0 offset
  forPercRast[forPercRast<perc]<-0 #set unchosen forest hexes to 0 offset
  forPercRast[forPercRast>=perc]<-15.52 #set forest height in random hexes
  forPercRast[lu==10]<-4.8 #add bajo tree height
  
  dem_model<-dem3+forPercRast
  remove(forPercRast)
  
  #get the coordinates of the structure that will be your observer point
  startSquareCoords <- squares[struct,] %>% st_coordinates()
  
  test <- terra::viewshed(x=dem_model,loc=squares[struct,] %>% st_coordinates(),observer=obs)
  remove(dem_model)
return (test)
}

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->



<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuXG50bV9zaGFwZShzcXVhcmVzKSt0bV9kb3RzKCkrXG50bV9zaGFwZShzcXVhcmVzKSt0bV90ZXh0KCdpZCcpXG5gYGAifQ== -->

```r

tm_shape(squares)+tm_dots()+
tm_shape(squares)+tm_text('id')

use this chunk only to test out a single observer at one square platform

#run the viewshed model for one point. takes about 10 seconds to run one

#struct is the id of the structures

#obs is observer height in meters

#setting perc to 100 means that 0% of pixels will be forested
#setting perc to 75 means that 25% of pixels will be forested
#setting perc to 50 means that 50% of pixels will be forested

viewTest <- randFor(struct=10,perc=50,obs=50)

plot(viewTest)

writeRaster(viewTest,'./data/singleview_struc10_50p_50m.tif',overwrite=TRUE)

This chunk below is what you’ll run for your all 17 square viewshed

modify parameters, run the code and then save the output with a different name

```r
#this chunk runs through all the structures and adds them together. this is exactly like what ArcGIS does when you have multiple observer points

#takes about 4 minutes

#this is what you'll do for each of the nine scenarios

allView<-sapply(squares$id,
                   FUN=randFor,
                   perc=75,
                   obs=50)  %>% rast() %>% sum()

#plot your results
plot(allView)

# export for use in ArcGIS - 
# change the output name to match model parameters from above
writeRaster(allView,'./data/view_25p_50m.tif',overwrite=TRUE)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->





***Let's talk about this one. I think it's pretty neat but we might want to save it for after presentation***


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuI3RoaXMgaXMgc29tZXRoaW5nIEkganVzdCB0ZXN0ZWQuIFNpbmNlIHRoZSByZXN1bHRzIGFyZSByYW5kb20sIHlvdSdsbCBnZXQgYSBkaWZmZXJlbnQgbWFwIGV2ZXJ5IHRpbWUuIGJ5IHJ1bm5pbmcgdGhlIHNhbWUgcGFyYW1ldGVycyAxMCB0aW1lcywgeW91IGNhbiBnZXQgYW4gaWRlYSBvZiB3aGljaCBwaXhlbHMgY2FuIGJlIHNlZW4gbm8gbWF0dGVyIHRoZSBmb3Jlc3QgY29uZmlndXJhdGlvblxuXG4jcnVuIHRoaXMgZm91ciB0aW1lczogcGVyYyA9IDc1ICgyNSUpIGFuZCBwZXJjID0gNTAgKDUwJSkgYXQgYm90aCBoZWlnaHRzICg1IGFuZCA1MCBtKVxuXG5mb3IgKGkgaW4gMToxNykge1xuI3J1biAxMCB0aW1lcyBmb3IgZWFjaCBzaXRlXG5yZXBlYXRTaW5nbGU8LXNhcHBseShyZXAoaSwxMCksXG4gICAgICAgICAgICAgICAgICAgIEZVTj1yYW5kRm9yLFxuICAgICAgICAgICAgICAgICAgICBwZXJjPTUwLFxuICAgICAgICAgICAgICAgICAgICBvYnM9NTApICU+JSByYXN0KCkgJT4lIHN1bSgpXG5cbiNjcmVhdGUgYSBuYW1lIHRvIHBhc3RlIGludG8gd3JpdGVSYXN0ZXJcbm5hbWUgPC0gcGFzdGUoJy4vZGF0YS9zaXRlJyxpLCdfMTB4XzUwcF81MG0udGlmJyxzZXA9JycpXG5cbiN3cml0ZSB0aGUgcmFzdGVyXG53cml0ZVJhc3RlcihyZXBlYXRTaW5nbGUsbmFtZSxvdmVyd3JpdGU9VFJVRSlcblxuI3JlbW92ZSB0aGUgZmlsZSBhbmQgdGhlbiBzdGFydCBhZ2FpblxucmVtb3ZlKHJlcGVhdFNpbmdsZSlcbn1cblxuYGBgIn0= -->

```r
#this is something I just tested. Since the results are random, you'll get a different map every time. by running the same parameters 10 times, you can get an idea of which pixels can be seen no matter the forest configuration

#run this four times: perc = 75 (25%) and perc = 50 (50%) at both heights (5 and 50 m)

for (i in 1:17) {
#run 10 times for each site
repeatSingle<-sapply(rep(i,10),
                    FUN=randFor,
                    perc=50,
                    obs=50) %>% rast() %>% sum()

#create a name to paste into writeRaster
name <- paste('./data/site',i,'_10x_50p_50m.tif',sep='')

#write the raster
writeRaster(repeatSingle,name,overwrite=TRUE)

#remove the file and then start again
remove(repeatSingle)
}