
#1. Prerequisites
install.packages("ncdf4") # to handle netcdf files
install.packages("rgdal") # to handle coordinate system conversion
library(ncdf4)
library(rgdal)


#Replace the source with the source of the downloaded slohu3h.nc file's location
setwd("E:/Google Drive/Tritium/database")

#File name
fname<-"AP3H_v1.nc"

#Import the file
obsdata=nc_open( fname, write=FALSE, readunlim=TRUE, verbose=FALSE, auto_GMT=TRUE, suppress_dimvals=FALSE )

#2. Provide coordinates
#Provide the Lon and Lat values to find the nearest grid cell
#These coordinates should be in (WGS 84 / Pseudo-Mercator - Spherical Mercator; EPSG:3857)
#If not, first use the following code to convert standard Lon-Lat (WGS 84 - WGS84 - World Geodetic System 1984; EPSG:4326) to convert your coordinates to the desired one.

#2.1 Coordinate conversion
#Example to convert Lat, Lon
orig_coords <- data.frame(lon=19.07291, lat=46.094612)
coordinates(orig_coords) <- c('lon', 'lat')

#Determine the projection of the lat-long coordinates, by default it is EPSG:4326
proj4string(orig_coords) <- CRS("+init=epsg:4326")
print(summary(orig_coords))

#Convert the coordinates to the used metric system (EPSG:3857)
Metric_coords<-spTransform(orig_coords,CRS("+init=epsg:3857"))
print(summary(Metric_coords))

#Extract the converted coordinates and use them in the exctraction of the 3H time series
coords<-Metric_coords@coords
X <- coords[1]  # longitude of location
Y <- coords[2]  # latitude of location

#Determine the first time point from which you want to extract the data
obsdatadates <- as.Date(obsdata$dim$time$vals, origin = '1976-1-1 00:00:00')

# look for closest long
# look for closest lat
obsoutput <- ncvar_get(obsdata, varid = 'slohu3h',
            start= c(which.min(abs(obsdata$dim$Y$vals - Y)),which.min(abs(obsdata$dim$X$vals - X)),  
            1), count = c(1,1,-1))

#Gives back the coordinates of the closest grid cell
Y_grid=obsdata$dim$Y$vals[which.min(abs(obsdata$dim$Y$vals - Y))]
X_grid=obsdata$dim$X$vals[which.min(abs(obsdata$dim$X$vals - X))]


#Output plot and table in csv
datafinal <- data.frame(dates= obsdatadates, obs = obsoutput)

pdf(paste(X,' & ',Y,'3Htime series.pdf', sep = ''), 11.69, 8.27)
plot(datafinal)
dev.off()

write.table(datafinal,file = "3H.csv", sep=";",row.names=F)


