NetCDFpointExtraction.r

install.packages("ncdf4")
install.packages("purrr")

library(ncdf4)
library(purrr)
setwd("/Users/marketa/Library/Mobile Documents/com~apple~CloudDocs/School/Thesis")

#reading the NetCDF file and getting specific variables
f <- nc_open("20180527_1441_historic_SPI_6_month.nc")
time <- ncdf4::ncvar_get(f, varid = "time")
spi <- ncdf4::ncvar_get(f, varid = "spi")
lat <- ncdf4::ncvar_get(f, varid = "latitude")
lon <- ncdf4::ncvar_get(f, varid = "longitude")

#setting the date format from a number to a normal date
time_d <- as.Date(time, format="%j", origin=as.Date("1900-01-01"))
time_years <- format(time_d, "%Y")
time_months <- format(time_d, "%m")
time_year_months <- format(time_d, "%Y-%m")

#reading the coordinates of the 1000 points that need to be assigned SPI value
rnd_pts <- read.csv("RNDM_points_coordinates.csv")

#finding the closest coordinate value
#here needs to be a for loop changing the row in matrix below that can also change the 1 in point1
data_frames <- map(seq(500, 1000), function(i) {
  rndx <- rnd_pts[i,3]
  rndy <- rnd_pts[i,4]
  positionx <- which.min(abs(lon - rndx))
  positiony <- which.min(abs(lat - rndy))
  x <- lon[positionx]
  y <- lat[positiony]
  
  point_2000 <- spi[lon == x, lat == y, time_years == "2000"]
  point_2001 <- spi[lon == x, lat == y, time_years == "2001"]
  point_2002 <- spi[lon == x, lat == y, time_years == "2002"]
  point_2003 <- spi[lon == x, lat == y, time_years == "2003"]
  point_2004 <- spi[lon == x, lat == y, time_years == "2004"]
  point_2005 <- spi[lon == x, lat == y, time_years == "2005"]
  point_2006 <- spi[lon == x, lat == y, time_years == "2006"]
  point_2007 <- spi[lon == x, lat == y, time_years == "2007"]
  point_2008 <- spi[lon == x, lat == y, time_years == "2008"]
  point_2009 <- spi[lon == x, lat == y, time_years == "2009"]
  point_2010 <- spi[lon == x, lat == y, time_years == "2010"]
  point_2011 <- spi[lon == x, lat == y, time_years == "2011"]
  point_2012 <- spi[lon == x, lat == y, time_years == "2012"]
  point_2013 <- spi[lon == x, lat == y, time_years == "2013"]
  point_2014 <- spi[lon == x, lat == y, time_years == "2014"]
  point_2015 <- spi[lon == x, lat == y, time_years == "2015"]
  point_2016 <- spi[lon == x, lat == y, time_years == "2016"]
  
  return(data.frame(point_2000, point_2001, point_2002, point_2003, point_2004, point_2005, point_2006, point_2007, point_2008, point_2009, point_2010, point_2011, point_2012, point_2013, point_2014, point_2015, point_2016))
  
})

#would it be possible to remove all point1_XXXX to save some space but keep point1?
#end of for loop
#combine all data frames
write.csv(data_frames, file="/Users/marketa/Library/Mobile Documents/com~apple~CloudDocs/School/Thesis/spi6m_500_1000.csv")

install.packages("ncdf4")
install.packages("purrr")

library(ncdf4)
library(purrr)
setwd("C:/Users/mpodebradska2/Desktop/School/Thesis/Indices_analysis/EDDI")

#reading the NetCDF file and getting specific variables
f <- nc_open("EDDI_03wk_2016.nc")
time <- ncdf4::ncvar_get(f, varid = "time")
EDDI <- ncdf4::ncvar_get(f, varid = "eddi")
lat <- ncdf4::ncvar_get(f, varid = "latitude")
lon <- ncdf4::ncvar_get(f, varid = "longitude")

#setting the date format from a number to a normal date
time_d <- as.Date(time, format="%j", origin=as.Date("1900-01-01"))
time_years <- format(time_d, "%Y")
time_months <- format(time_d, "%m")
time_year_months <- format(time_d, "%Y-%m")

#reading the coordinates of the 1000 points that need to be assigned SPI value
rnd_pts <- read.csv("RNDM_points_coordinates.csv")

#finding the closest coordinate value
#here needs to be a for loop changing the row in matrix below that can also change the 1 in point1
data_frames <- map(seq(1, 1000), function(i) {
rndx <- rnd_pts[i,3]
rndy <- rnd_pts[i,4]
positionx <- which.min(abs(lon - rndx))
positiony <- which.min(abs(lat - rndy))
x <- lon[positionx]
y <- lat[positiony]

point <- EDDI[lon == x, lat == y, time_years == "2016"]

return(data.frame(point))

})

#combine all data frames
write.csv(data_frames, file="C:/Users/mpodebradska2/Desktop/School/Thesis/Indices_analysis/EDDI/EDDI2016.csv")