Raster_Extraction.R

easypackages::packages(c("tidyverse", 'sf', 'raster','sp','dplyr','spatialEco',
                         'exactextractr','tibble', 'here', "readr", "margrittr","tydr",
                         "readxl", "foreign", "data.table", "spData","spDataLarge", ))

remotes::install_github("Nowosad/spDataLarge")
setwd(here("Input_data"))

# Raster to point extraction  ---------------------------------------------
# Making SPAM spatial explicit 
library(readr)
dat_spam=read_csv("C:/Users/molch/OneDrive - Wageningen University & Research/PhD_WJS/Academic/RQ1/Data_analysis/Raw_dat/Zonation/SPAM/2010/spam2010v2r0_global_yield.csv/spam2010V2r0_global_Y_TA.csv")
dat_spam=select(dat_spam, 'x','y', 'alloc_key')

# Make SPAM spatial 
# projcrs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
# proj_mercator = "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
proj = "+proj=longlat +datum=WGS84 +no_defs +units=m" # taken from raster file below 
proj_meter = "+datum=WGS84 +units=m +no_defs"
dat_cent <- st_as_sf(x = dat_spam,                         
                     coords = c("x", "y"),
                     crs = proj)

# dat_long = st_as_sf(coords = c("lon", "lat"))

# dat_cent = filter(dat_cent, iso3=="FIN")
# dat_cent=dat_cent[1:10,]


# import raster datasets 
dat_rast = raster("D:/Raw_dat/Raw_dat/Phosphorus/nitrogen-fertilizer-global-geotif/NFertilizer_global_geotif/nfertilizer_global.tif")
dat_rast = raster("D:/Raw_dat/Raw_dat/Phosphorus/nitrogen-fertilizer-global-geotif/NFertilizer_global_geotif/nfertilizer_global.tif")
dat_rast = raster("D:/Raw_dat/Raw_dat/Phosphorus/nitrogen-fertilizer-global-geotif/NFertilizer_global_geotif/nfertilizer_global.tif")
dat_rast_SOC = raster("D:/Raw_dat/Raw_dat/Soil maps/Opensoilmaps/sol_organic.carbon_usda.tif")
dat_rast_ISCRIC = raster("D:/Raw_dat/Raw_dat/Soil maps/wise_05min_v12/wise_05min_v12/Grid/smw5by5min/w001001.adf")
dat_clim=getData('worldclim', var='bio', res=5) 
dat_rast_meantemp = dat_clim[[1]]

plot(clim[[1:4]])
# Import Polygons

# Point extraction of raster values and add it to the SPAM point --------
funrast = function(r){                  #insert the raster dataset from which you would like to extract values 
  r = raster::extract(r,                # raster layer
                      dat_cent,         #  SF with centroids for buffer
                      buffer = 5000,    # buffer size, units depend on CRS. Here in meters.10x10km(~5km radius)
                      fun=mean,         # what value to extract
                      df=TRUE)          # creating a df TRUE
  r$ID = dat_cent$alloc_key # assigning the identifier of the original df to the extracted df
  dat_spam <- merge(dat_spam, r, by.x = 'alloc_key', by.y = 'ID')# merge the two files together 
  }

system.time(
  funrast(dat_rast)
  )

# Extraction
ex_Nfert_spam = funrast(dat_rast)
ex_soil_spam = funrast(dat_rast_ISCRIC)
ex_SOC_spam = funrast(dat_rast_SOC)
# start SOC extraction at 7.30, 19.th

filter(test, nfertilizer_global>50)
summary(test)


# Zonal statistics  -------------------------------------------------------
# Load polygons
wd=setwd("C:/Users/molch/OneDrive - Wageningen University & Research/PhD_WJS/Academic/RQ1/Data_analysis/Raw_dat/Zonation/GADM/gadm36_shp")

# Load and check shapefile 
shape <- read_sf(dsn = wd, layer = "gadm36") # multipolygons containing all admin levels globally
# class(shape)
# shape$GID
shape=select(shape, 'UID', 'geometry', 'GID_0')
# y=filter(x, GID_0 == 'ALG') 

# Load raster 
rast=dat_rast = raster("D:/Raw_dat/Raw_dat/Phosphorus/nitrogen-fertilizer-global-geotif/NFertilizer_global_geotif/nfertilizer_global.tif")


# Perform zonal statistics ------------------------------------------------
FunZone=function(rast, shape){ # raster and vector
  c=zonal.stats(shape, d, stats = c("mean")) 
  c$UID <- shape$UID  # assigning the identifier of the original df to the extracted df
  dat <- merge(shape, c, by.x = 'UID', by.y = 'UID') # merge the two files togethere
}


dat$mean.nfertilizer_global
summary(dat$mean.nfertilizer_global)



c=zonal.stats(y, dat_rast, stats = c("mean")) 
c$UID <- x$UID  # assigning the identifier of the original df to the extracted df
dat <- merge(shape, c, by.x = 'UID', by.y = 'UID') # merge the two files together 

# Point = overlay(rast, df) #does not work

cent_mean_FIN <- raster::extract(rast,        # raster layer
                                 ds,   #           #  SF with centroids for buffer
                                 buffer = 5000,   # buffer size, units depend on CRS. Here in meters.10x10km(~5km radius)
                                 fun=mean,         # what value to extract
                                 df=TRUE)          # creating a df TRUE

cent_mean_FIN$ID = ds$alloc_key # assigning the identifier of the original df to the extracted df
centroids <- merge(dat, cent_mean_FIN, by.x = 'alloc_key', by.y = 'ID') # merge the two files together