pacman::p_load( "exactextractr","raster","ggplot2","unikn","mapview",
"gridExtra","rgdal","fields",
"RColorBrewer","ncdf4","rasterVis", "moments", "tictoc","rgeos",
"snow","future.apply","sp","tigris","sf","tidyr","rcartocolor","parallel","snow","future.apply")
require(maptools)
data(wrld_simpl)
nc = nc_open("C:/Users/tamis/Downloads/FDSI.nc")
time_var= ncvar_get(nc,'Time')
Latitude= ncvar_get(nc,'Latitude')
Longitude= ncvar_get(nc,'Longitude')
FDSI= ncvar_get(nc,'FDSI')
nc_close(nc)
FDSIbrk= brick("C:/Users/tamis/Downloads/FDSI.nc")
IPCC_shp = rgdal::readOGR("C:/Users/tamis/Downloads/SampleData-master/SampleData-master/CMIP_land/CMIP_land.shp",verbose = FALSE)
head(IPCC_shp@data)
global_crop1 = crop(coastlines,extent(IPCC_shp))
raster::plot(IPCC_shp, add=TRUE) # Add IPCC land regions in blue color
raster::plot(coastlines, add=TRUE) # Add global coastline
#~ Extract feature data as data frame
In order to make the operation faster we have to extract the desired layers.
Subset the raster data to the dates of interest
nl=nlayers(FDSIbrk)
start_date <- which(getZ(FDSIbrk) == '20210821')
end_date <- which(getZ(FDSIbrk) == '20210822')
FDSI_brk_sub <- FDSIbrk[[start_date:end_date]]
featureData = exact_extract(FDSIbrk, # r brick
IPCC_shp, # Convert shapefile to sf (simple feature)
force_df = FALSE, # Output as a data.frame
include_xy = FALSE, # Incluse cell lat-long in output?
fun = 'mean', # Specify any function to apply for each feature extracted data
progress = TRUE) # Progressbar
names(featureData)
Plot mean fdsi map for IPCC
#create a folder where it can be save.
setwd("G:/.shortcut-targets-by-id/1_munBiOkEi1EyelAxqX1FaQ56WRqt47Z/FLASH_shared/Ashely Johnson")
dir.create("Img_Assignment_2")
setwd("G:/.shortcut-targets-by-id/1_munBiOkEi1EyelAxqX1FaQ56WRqt47Z/FLASH_shared/Ashely Johnson/Img_Assignment_2")
Create the preffered colors.
mypal= brewer.pal(11,"Spectral") # Max colors available in the palette=11
mypal= colorRampPalette(mypal)(20) # Expand color palette to 20plot(venshp)
for(i in 2678:nlayers(FDSI_brk_sub)) {
mean_map=ggplot() +
geom_sf(data = (IPCC_shp), # CONUS shp as sf object (simple feature)
aes(fill = featureData[[i]])) + # Plot fill color= mean soil moisture
scale_fill_gradientn(colors=rev(mypal), # Use user-defined colormap
name = "Mean FDSI", # Name of the colorbar
na.value = "transparent", # transparent NA cells
labels = (c("0", "0.2", "0.4", "0.6", "0.8")), # Labels of colorbar
breaks = seq(0,0.8,by =0.2), # Set breaks of colorbar
limits = c(0,0.8))+ # To invert color, use -1
coord_sf(crs = 2163)+ # Reprojecting polygon 4326 or 3083
theme_void() + # Plot theme. Try: theme_bw
theme(legend.position = c(0.2, 0.1),
legend.direction = "horizontal",
legend.key.width = unit(5, "mm"),
legend.key.height = unit(4, "mm")) + ggtitle(names(featureData[i]))
Saving the plot to disk:
ggsave(plot = mean_map, # ggplot object to be exported
filename =paste(names(featureData[i]), "png", sep="."), # Export file name
bg = "white", # Background color
units="in", dpi=300, # Export dim unit & Dots per inch (resolution)
width=10, height=7) # Height, width
}
mean_map