Introduction to NetCDF handling: A test case using CM-SAF data
Maria Karypidou
24 October 2018
This is an introductory tutorial on some basic handling of a NetCDF file in R. For this purpose we will use some CM-SAF data [1], which are freely available online. In order to download your data from the website, you will first need to set up a user account here [2]. Then, you can start your quiry and request the product you want.
In R, there are two ways one can read in a NetCDF file. We will go through both and then you decide which one you will adopt.
So, let’s start, by loading in the packages/libraries we will work with! Remember, you must first “install.packages” in you R environment and then load the libraries. You install the libraries once, but you load them in each R session anew.
library("raster")## Loading required package: splibrary("RColorBrewer")
library("viridis")## Loading required package: viridisLitelibrary("maptools")## Checking rgeos availability: TRUElibrary("rgdal")## rgdal: version: 1.3-4, (SVN revision 766)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/user/Documents/R/win-library/3.5/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/user/Documents/R/win-library/3.5/rgdal/proj
## Linking to sp version: 1.3-1library("maps")
library("mapdata")
library("rgeos")## rgeos version: 0.3-28, (SVN revision 572)
## GEOS runtime version: 3.6.1-CAPI-1.10.1 r0
## Linking to sp version: 1.3-1
## Polygon checking: TRUElibrary("ggplot2")
library("latticeExtra")## Loading required package: lattice##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layerlibrary("grid")
library("colorspace")##
## Attaching package: 'colorspace'## The following object is masked from 'package:raster':
##
## RGBlibrary("rasterVis")
library("ncdf4")
library("colorRamps")
library("mapproj")
library("fields")## Loading required package: spam## Loading required package: dotCall64## Spam version 2.2-0 (2018-06-19) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## See www.image.ucar.edu/~nychka/Fields for
## a vignette and other supplements.library("ncdf4.helpers")
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:rgeos':
##
## intersect, setdiff, union## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("ggmap")
library("viridis")
library("weathermetrics")# We need to define where our data sit. You must change the path here and put your own path accordingly!
ncpath = "C:/Users/user/Documents/Tutorials/Tutorial2/"
#
# and how are they called
ncname = "SIS_sarah_daily_2007_2008_greece_monmean"
#
# and define what type of files are they?
# ".nc" is how NetCDF files are represented
# The "paste" command in R is, I think, a pure blessing! Find more about "paste" here [3]. It is expremely simple # and extremely useful!
# Now, let's wrap them all together
ncfname = paste(ncpath, ncname, ".nc", sep="")
ncfname## [1] "C:/Users/user/Documents/Tutorials/Tutorial2/SIS_sarah_daily_2007_2008_greece_monmean.nc"#
# The name of the variable we will work with
dname = "SIS" Let’s open our file and see what it contains. When we hit “print” we will see all the attributes thats our NetCDF file containes. Before proceeding with any sort of analysis, it is VERY IMPORTANT that we know the attributes of our file.
# By "nc_open" we ask R to open the netcdf file on the exact location we defined in the previous lines
ncin = nc_open(ncfname)
#
# and we will print "ncin", so that we see what it contains.
#
print(ncin)## File C:/Users/user/Documents/Tutorials/Tutorial2/SIS_sarah_daily_2007_2008_greece_monmean.nc (NC_FORMAT_NETCDF4_CLASSIC):
##
## 2 variables (excluding dimension variables):
## double time_bnds[bnds,time]
## short SIS[lon,lat,time]
## standard_name: surface_downwelling_shortwave_flux_in_air
## long_name: Surface Downwelling Shortwave Radiation
## units: W m-2
## _FillValue: -999
## missing_value: -999
##
## 4 dimensions:
## lon Size:313
## standard_name: longitude
## long_name: longitude
## units: degrees_east
## axis: X
## lat Size:168
## standard_name: latitude
## long_name: latitude
## units: degrees_north
## axis: Y
## time Size:24 *** is unlimited ***
## standard_name: time
## bounds: time_bnds
## units: days since 1983-1-1 00:00:00
## calendar: proleptic_gregorian
## axis: T
## bnds Size:2
##
## 38 global attributes:
## CDI: Climate Data Interface version 1.7.1 (http://mpimet.mpg.de/cdi)
## history: Wed Oct 24 18:38:02 2018: cdo monmean SIS_sarah_daily_2007_2008_greece.nc SIS_sarah_daily_2007_2008_greece_monmean.nc
## Wed Oct 24 18:30:05 2018: cdo sellonlatbox,19,34.6,34.6,43 SIS_sarah_daily_2007_2008.nc SIS_sarah_daily_2007_2008_greece.nc
## Wed Oct 24 17:50:56 2018: cdo selvar,SIS sis_sarah_daily_2007_2008.nc SIS_sarah_daily_2007_2008.nc
## Wed Oct 24 17:15:14 2018: cdo selyear,2007/2008 sis_sarah_daily_1990_2008.nc sis_sarah_daily_2007_2008.nc
## Fri Apr 20 10:10:41 2018: cdo mergetime cmSIS_1990.nc cmSIS_1991.nc cmSIS_1992.nc cmSIS_1993.nc cmSIS_1994.nc cmSIS_1995.nc cmSIS_1996.nc cmSIS_1997.nc cmSIS_1998.nc cmSIS_1999.nc cmSIS_2000.nc cmSIS_2001.nc cmSIS_2002.nc cmSIS_2003.nc cmSIS_2004.nc cmSIS_2005.nc cmSIS_2006.nc cmSIS_2007.nc cmSIS_2008.nc sis_sarah_daily_1990_2008.nc
## Fri Apr 20 03:53:32 2018: cdo mergetime SISdm200801010000003231000101MA.nc SISdm200801020000003231000101MA.nc SISdm200801030000003231000101MA.nc SISdm200801040000003231000101MA.nc SISdm200801050000003231000101MA.nc SISdm200801060000003231000101MA.nc SISdm200801070000003231000101MA.nc SISdm200801080000003231000101MA.nc SISdm200801090000003231000101MA.nc SISdm200801100000003231000101MA.nc SISdm200801110000003231000101MA.nc SISdm200801120000003231000101MA.nc SISdm200801130000003231000101MA.nc SISdm200801140000003231000101MA.nc SISdm200801150000003231000101MA.nc SISdm200801160000003231000101MA.nc SISdm200801170000003231000101MA.nc SISdm200801180000003231000101MA.nc SISdm200801190000003231000101MA.nc SISdm200801200000003231000101MA.nc SISdm200801210000003231000101MA.nc SISdm200801220000003231000101MA.nc SISdm200801230000003231000101MA.nc SISdm200801240000003231000101MA.nc SISdm200801250000003231000101MA.nc SISdm200801260000003231000101MA.nc SISdm200801270000003231000101MA.nc SISdm200801280000003231000101MA.nc
## institution: EUMETSAT/CMSAF
## Conventions: CF-1.6, ACDD-1.3
## title: Surface Solar Radiation Climate Data Record using Heliosat, Edition 2 (SARAH-2); CM SAF
## summary: This file contains data from the CM SAF Surface Solar Radiation Climate Data Record using Heliosat, Edition 2 (SARAH-2).
## id: DOI:10.5676/EUM_SAF_CM/SARAH/V002
## creator_name: DE/DWD
## creator_email: contact.cmsaf@dwd.de
## creator_url: http://www.cmsaf.eu/
## creator_type: institution
## publisher_name: EUMETSAT/CMSAF
## publisher_email: contact.cmsaf@dwd.de
## publisher_url: http://www.cmsaf.eu/
## publisher_type: institution
## references: http://dx.doi.org/10.5676/EUM_SAF_CM/SARAH/V002
## keywords_vocabulary: GCMD Science Keywords, Version 8.1
## keywords: EARTH SCIENCE > ATMOSPHERE > ATMOSPHERIC RADIATION > INCOMING SOLAR RADIATION
## project: Satellite Application Facility on Climate Monitoring (CM SAF)
## standard_name_vocabulary: Standard Name Table (v28, 07 January 2015)
## date_created: 2016-11-03T21:45:46Z
## geospatial_lat_units: degrees_north
## geospatial_lat_min: -65
## geospatial_lat_max: 65
## geospatial_lon_units: degrees_east
## geospatial_lon_min: -65
## geospatial_lon_max: 65
## geospatial_lat_resolution: 0.05 degree
## geospatial_lon_resolution: 0.05 degree
## time_coverage_duration: P1D
## time_coverage_resolution: P1D
## platform_vocabulary: GCMD Platforms, Version 8.1
## platform: Earth Observation Satellites > METEOSAT
## instrument_vocabulary: GCMD Instruments, Version 8.1
## instrument: MVIRI > Meteosat Visible Infra-Red Imager
## product_version: 2.0
## frequency: mon
## CDO: Climate Data Operators version 1.7.1 (http://mpimet.mpg.de/cdo)Our NetCDf file contains latitude, longitude and time as separate variables. Let’s make R read those and let’s store them to separate variables.
lon = ncvar_get(ncin,"lon")
nlon = dim(lon)
head(lon)## [1] 19.00 19.05 19.10 19.15 19.20 19.25summary(lon)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 19.0 22.9 26.8 26.8 30.7 34.6#
lat = ncvar_get(ncin,"lat")
nlat = dim(lat)
head(lat)## [1] 34.65 34.70 34.75 34.80 34.85 34.90summary(lat)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 34.65 36.74 38.83 38.83 40.91 43.00#
print(c(nlon,nlat))## [1] 313 168#
#
time = ncvar_get(ncin,"time")
time## [1] 8781.0 8810.5 8840.0 8870.5 8901.0 8931.5 8962.0 8993.0 9023.5 9054.0
## [11] 9084.5 9115.0 9146.0 9176.0 9206.0 9236.5 9267.0 9297.5 9328.0 9359.0
## [21] 9389.5 9420.0 9450.5 9481.0Now let’s get our SIS variable
tmp_array = ncvar_get(ncin,dname)
dlname = ncatt_get(ncin,dname,"long_name")
dunits = ncatt_get(ncin,dname,"units")
fillvalue = ncatt_get(ncin,dname,"_FillValue")
dim(tmp_array)## [1] 313 168 24From the above line we see that R says that our “ncin” file has 24 time steps. This is correct as we have loaded monthly data for two years.
Now let’s get our SIS variable
greece = getData('GADM', country='GR', level=0)
turkey = getData('GADM', country='TR', level=0)
bulgaria = getData('GADM', country='BG', level=0)
albania = getData('GADM', country='AL', level=0)
#
# This is how our shapefile looks like
plot(greece)
m = 1
tmp_slice = tmp_array[,,m]
image.plot(lon,lat,tmp_slice, col=rev(rainbow(100)))
plot(greece, add = TRUE)
plot(turkey, add = TRUE)
plot(bulgaria, add = TRUE)
plot(albania, add = TRUE)
Now, let’s assume that we do not want to plot the map of the whole area, but we only want to get the values for a selected location, Thessaloniki for instance. We know that the latitude of Thessaloniki is 40.736851 and it’s Longitude is 22.920227. Somehow, we have to search for this location over our whole NetCDF file and select the grid box that is placed closer to the coordinates of Thessaloniki.
sis_time = nc.get.time.series(ncin, v = "SIS",
time.dim.name = "time")
#
lon_index = which.min(abs(lon - 22.920227))
lat_index = which.min(abs(lat - 40.736851))
sis_salonica = nc.get.var.subset.by.axes(ncin, "SIS",
axis.indices = list(X = lon_index,
Y = lat_index))
data_frame(time = sis_time,
sis_salonica = as.vector(sis_salonica)) %>%
mutate(time = as.Date(format(time, "%Y-%m-%d"))) %>%
ggplot(aes(x = time, y = sis_salonica)) +
geom_line() +
xlab("Date") + ylab("Shortwave Solar Radiation @Surface (SIS) (W/m^2)") +
ggtitle("CM-SAF SIS data for 2007-2008 (mean monthly values)",
subtitle = "At Thessaloniki") +
theme_classic()
However, there is a much quicker way of working things around, using raster stacks.
my_stack = stack("C:/Users/user/Documents/Tutorials/Tutorial2/SIS_sarah_daily_2007_2008_greece_monmean.nc")
#
# Raster stacks have layers
nlayers(my_stack)## [1] 24# So, you can define which layer you want to plot
plot(my_stack[[1]], main="My first R map")
plot(greece, add = TRUE)
plot(turkey, add = TRUE)
plot(bulgaria, add = TRUE)
plot(albania, add = TRUE)
#
summary(my_stack)## X2007.01.16 X2007.02.14 X2007.03.16 X2007.04.15 X2007.05.16
## Min. 37 45 107 133 169
## 1st Qu. 84 96 154 234 255
## Median 102 109 181 248 275
## 3rd Qu. 114 128 198 260 287
## Max. 136 151 227 285 314
## NA's 0 0 0 0 0
## X2007.06.15 X2007.07.16 X2007.08.16 X2007.09.15 X2007.10.16
## Min. 227 269 201 144 89
## 1st Qu. 303 327 273 204 129
## Median 316 335 295 238 154
## 3rd Qu. 330 340 304 255 170
## Max. 357 366 332 277 208
## NA's 0 0 0 0 0
## X2007.11.15 X2007.12.16 X2008.01.16 X2008.02.15 X2008.03.16
## Min. 34 25 25 60 79
## 1st Qu. 80 62 77 117 155
## Median 96 77 94 132 175
## 3rd Qu. 109 88 108 142 199
## Max. 145 113 133 170 231
## NA's 0 0 0 0 0
## X2008.04.15 X2008.05.16 X2008.06.15 X2008.07.16 X2008.08.16
## Min. 135 170 195 217 226
## 1st Qu. 192 271 306 308 285
## Median 221 289 328 332 295
## 3rd Qu. 245 306 340 340 303
## Max. 286 340 367 363 327
## NA's 0 0 0 0 0
## X2008.09.15 X2008.10.16 X2008.11.15 X2008.12.16
## Min. 136 104 57 30
## 1st Qu. 183 144 84 62
## Median 210 163 104 72
## 3rd Qu. 226 175 118 87
## Max. 255 202 149 118
## NA's 0 0 0 0min(my_stack)## class : RasterLayer
## dimensions : 168, 313, 52584 (nrow, ncol, ncell)
## resolution : 0.05, 0.04999999 (x, y)
## extent : 18.975, 34.625, 34.625, 43.025 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 25, 113 (min, max)max(my_stack)## class : RasterLayer
## dimensions : 168, 313, 52584 (nrow, ncol, ncell)
## resolution : 0.05, 0.04999999 (x, y)
## extent : 18.975, 34.625, 34.625, 43.025 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 269, 367 (min, max)#
wor = readOGR("F:/Shapefiles/zipped/neighbors/neighbors.shp")## OGR data source with driver: ESRI Shapefile
## Source: "F:\Shapefiles\zipped\neighbors\neighbors.shp", layer: "neighbors"
## with 8 features
## It has 11 fields
## Integer64 fields read as strings: POP2005akto = readOGR("F:/Shapefiles/zipped/akto_wgs84/akto_wgs84.shp")## OGR data source with driver: ESRI Shapefile
## Source: "F:\Shapefiles\zipped\akto_wgs84\akto_wgs84.shp", layer: "akto_wgs84"
## with 5452 features
## It has 2 fieldsgr = readOGR("F:/Shapefiles/zipped/nomoi_wgs84/nomoi_wgs84.shp")## OGR data source with driver: ESRI Shapefile
## Source: "F:\Shapefiles\zipped\nomoi_wgs84\nomoi_wgs84.shp", layer: "nomoi_wgs84"
## with 55 features
## It has 8 fields
## Integer64 fields read as strings: POP#
r2 = crop(my_stack, extent(gr))
r3 = mask(r2, gr)
my_stack = stack(r3)
min(my_stack)## class : RasterLayer
## dimensions : 138, 205, 28290 (nrow, ncol, ncell)
## resolution : 0.05, 0.04999999 (x, y)
## extent : 19.375, 29.625, 34.825, 41.725 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 40, 104 (min, max)max(my_stack)## class : RasterLayer
## dimensions : 138, 205, 28290 (nrow, ncol, ncell)
## resolution : 0.05, 0.04999999 (x, y)
## extent : 19.375, 29.625, 34.825, 41.725 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 296, 367 (min, max)min(my_stack[[1]])## Warning in min(new("RasterLayer", file = new(".RasterFile", name = "",
## datanotation = "FLT4S", : Nothing to summarize if you provide a single
## RasterLayer; see cellStats## class : RasterLayer
## dimensions : 138, 205, 28290 (nrow, ncol, ncell)
## resolution : 0.05, 0.04999999 (x, y)
## extent : 19.375, 29.625, 34.825, 41.725 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : X2007.01.16
## values : 65, 136 (min, max)max(my_stack[[1]])## Warning in max(new("RasterLayer", file = new(".RasterFile", name = "",
## datanotation = "FLT4S", : Nothing to summarize if you provide a single
## RasterLayer; see cellStats## class : RasterLayer
## dimensions : 138, 205, 28290 (nrow, ncol, ncell)
## resolution : 0.05, 0.04999999 (x, y)
## extent : 19.375, 29.625, 34.825, 41.725 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : X2007.01.16
## values : 65, 136 (min, max)#
plot(my_stack[[1]], col=rainbow(50), axes = F, legend=T, zlim=c(60, 140),
las=1, main = "SIS from CM-SAF for 1/1/2007", cex.main=1, legend.shrink=1, legend.width=1)
plot(wor, col = "grey", add=T)
plot(akto, add=T)
lon_saloniki = 22.920227
lat_soliniki = 40.736851
ts <- extract(my_stack, cbind(lon_saloniki, lat_soliniki))
ts_saloniki = as.data.frame(ts)
class(ts_saloniki)## [1] "data.frame"head(ts_saloniki)## X2007.01.16 X2007.02.14 X2007.03.16 X2007.04.15 X2007.05.16 X2007.06.15
## 1 105 92 168 246 252 304
## X2007.07.16 X2007.08.16 X2007.09.15 X2007.10.16 X2007.11.15 X2007.12.16
## 1 331 279 206 129 85 74
## X2008.01.16 X2008.02.15 X2008.03.16 X2008.04.15 X2008.05.16 X2008.06.15
## 1 80 136 183 207 282 295
## X2008.07.16 X2008.08.16 X2008.09.15 X2008.10.16 X2008.11.15 X2008.12.16
## 1 316 288 191 148 76 69ts_saloniki = t(ts_saloniki)
print(ts_saloniki)## [,1]
## X2007.01.16 105
## X2007.02.14 92
## X2007.03.16 168
## X2007.04.15 246
## X2007.05.16 252
## X2007.06.15 304
## X2007.07.16 331
## X2007.08.16 279
## X2007.09.15 206
## X2007.10.16 129
## X2007.11.15 85
## X2007.12.16 74
## X2008.01.16 80
## X2008.02.15 136
## X2008.03.16 183
## X2008.04.15 207
## X2008.05.16 282
## X2008.06.15 295
## X2008.07.16 316
## X2008.08.16 288
## X2008.09.15 191
## X2008.10.16 148
## X2008.11.15 76
## X2008.12.16 69ts_saloniki = as.data.frame(ts_saloniki)
print(ts_saloniki)## V1
## X2007.01.16 105
## X2007.02.14 92
## X2007.03.16 168
## X2007.04.15 246
## X2007.05.16 252
## X2007.06.15 304
## X2007.07.16 331
## X2007.08.16 279
## X2007.09.15 206
## X2007.10.16 129
## X2007.11.15 85
## X2007.12.16 74
## X2008.01.16 80
## X2008.02.15 136
## X2008.03.16 183
## X2008.04.15 207
## X2008.05.16 282
## X2008.06.15 295
## X2008.07.16 316
## X2008.08.16 288
## X2008.09.15 191
## X2008.10.16 148
## X2008.11.15 76
## X2008.12.16 69ts_saloniki = as.data.frame(ts_saloniki)
print(ts_saloniki)## V1
## X2007.01.16 105
## X2007.02.14 92
## X2007.03.16 168
## X2007.04.15 246
## X2007.05.16 252
## X2007.06.15 304
## X2007.07.16 331
## X2007.08.16 279
## X2007.09.15 206
## X2007.10.16 129
## X2007.11.15 85
## X2007.12.16 74
## X2008.01.16 80
## X2008.02.15 136
## X2008.03.16 183
## X2008.04.15 207
## X2008.05.16 282
## X2008.06.15 295
## X2008.07.16 316
## X2008.08.16 288
## X2008.09.15 191
## X2008.10.16 148
## X2008.11.15 76
## X2008.12.16 69ts_saloniki$DATE = seq.Date(min(as.Date("2007-1-1")), max(as.Date("2008-12-31")), by="month")
ts_saloniki## V1 DATE
## X2007.01.16 105 2007-01-01
## X2007.02.14 92 2007-02-01
## X2007.03.16 168 2007-03-01
## X2007.04.15 246 2007-04-01
## X2007.05.16 252 2007-05-01
## X2007.06.15 304 2007-06-01
## X2007.07.16 331 2007-07-01
## X2007.08.16 279 2007-08-01
## X2007.09.15 206 2007-09-01
## X2007.10.16 129 2007-10-01
## X2007.11.15 85 2007-11-01
## X2007.12.16 74 2007-12-01
## X2008.01.16 80 2008-01-01
## X2008.02.15 136 2008-02-01
## X2008.03.16 183 2008-03-01
## X2008.04.15 207 2008-04-01
## X2008.05.16 282 2008-05-01
## X2008.06.15 295 2008-06-01
## X2008.07.16 316 2008-07-01
## X2008.08.16 288 2008-08-01
## X2008.09.15 191 2008-09-01
## X2008.10.16 148 2008-10-01
## X2008.11.15 76 2008-11-01
## X2008.12.16 69 2008-12-01min(ts_saloniki$V1)## [1] 69max(ts_saloniki$V1)## [1] 331colnames(ts_saloniki)[colnames(ts_saloniki)=="V1"] = "SIS"
head(ts_saloniki)## SIS DATE
## X2007.01.16 105 2007-01-01
## X2007.02.14 92 2007-02-01
## X2007.03.16 168 2007-03-01
## X2007.04.15 246 2007-04-01
## X2007.05.16 252 2007-05-01
## X2007.06.15 304 2007-06-01#
ggplot(ts_saloniki, aes(DATE))+
#
geom_line(aes(y=SIS, col="SIS"), size=0.8)+
xlab("Date") + ylab("Shortwave Solar Radiation @Surface (SIS) (W/m^2)") +
ggtitle("CM-SAF SIS data for 2007-2008 (mean monthly values)",
subtitle = "At Thessaloniki")
print("The End")## [1] "The End"I hope you enjoyed it! Please send me your feedback at karypidou@geo.auth.gr
[1] https://www.cmsaf.eu/EN/Home/home_node.html [2] https://wui.cmsaf.eu/safira/action/viewLogin;jsessionid=E609AF65C4BAF27E102129D198889D70.ku_1?menuName=NUTZER_HOME [3] https://www.rdocumentation.org/packages/base/versions/3.5.1/topics/paste