diff --git a/revdep/checks.rds b/revdep/checks.rds index f046d0fa..69985c0d 100644 Binary files a/revdep/checks.rds and b/revdep/checks.rds differ diff --git a/vignettes/Specified_stations_for_a_range_of_years.Rmd b/vignettes/Specified_stations_for_a_range_of_years.Rmd index 7e140dc3..3542b8a1 100644 --- a/vignettes/Specified_stations_for_a_range_of_years.Rmd +++ b/vignettes/Specified_stations_for_a_range_of_years.Rmd @@ -30,38 +30,39 @@ include: First we retrieve data from GADM.org that will provide the provincial spatial data for the survey area. We will then use this to find the centroid, which will be used to find the nearest stations. -```{r, RP Data, message=FALSE, eval=FALSE} + +```r library(raster) library(rgdal) RP <- getData(country = "Philippines", level = 1) - ``` Select the provinces involved in the survey and make a new object called `Central_Luzon`. -```{R, Central_Luzon, eval=FALSE} +```r Central_Luzon <- RP[RP@data$NAME_1 == "Pampanga" | RP@data$NAME_1 == "Tarlac" | RP@data$NAME_1 == "Pangasinan" | RP@data$NAME_1 == "La Union" | RP@data$NAME_1 == "Nueva Ecija" | RP@data$NAME_1 == "Bulacan", ] - ``` ## Dissolve Polygons and Find Centroid of Loop Survey Area Now that we have the provincial data imported, we will dissolve the polygons and find the centroid. -```{r, dissolve, message=FALSE, eval=FALSE} + +```r library(rgeos) Central_Luzon <- gUnaryUnion(Central_Luzon) centroid <- gCentroid(Central_Luzon) ``` -```{r, nearest_stations, message=FALSE, eval=FALSE} + +```r library(GSODR) library(readr) # Fetch station list from NCDC @@ -82,7 +83,6 @@ loop_stations <- station_meta[station_meta$STNID %in% loop_stations, ] loop_stations <- loop_stations[loop_stations$BEGIN <= 19591231 & loop_stations$END >= 20151231, ] - ``` ## Using get_GSOD to Fetch the Requested Station Files @@ -95,11 +95,11 @@ function. Be aware that it may result in incomplete data and error from the server. If it does this, see the following option for using the `format_GSOD` function. -```{r, fetch_weather, message=FALSE, eval=FALSE} + +```r get_GSOD(station = eval(parse(text = loop_stations[, 12])), years = 1960:2015, CSV = TRUE, dsn = "~/", filename = "Loop_Survey_Weather_1960-1969") - ``` ## Another Option @@ -116,7 +116,8 @@ directory, "gsod-YYYY.tar" and untar them locally and then use R to list the available files and select only the files for the stations of interest. Lastly, write the data to disk as a CSV file for saving and later use. -```{r, list_local_files, eval=FALSE} + +```r years <- 1960:2015 loop_stations <- eval(parse(text = loop_stations[, 12])) @@ -134,4 +135,44 @@ loop_data <- reformat_GSOD(file_list = local_files) write.csv(loop_data, file = "Loop_Survey_Weather_1960-1969") ``` +# Notes + +## Sources + +#### CHELSA climate layers +CHELSA (climatic surfaces at 1 km resolution) is based on a quasi-mechanistical +statistical downscaling of the ERA interim global circulation model +(Karger et al. 2016). ESA's CCI-LC cloud probability monthly averages are based +on the MODIS snow products (MOD10A2). + +#### EarthEnv MODIS cloud fraction + + +#### ESA's CCI-LC cloud probability + + +#### Elevation Values + +90m hole-filled SRTM digital elevation (Jarvis *et al.* 2008) was used +to identify and correct/remove elevation errors in data for station +locations between -60˚ and 60˚ latitude. This applies to cases here +where elevation was missing in the reported values as well. In case the +station reported an elevation and the DEM does not, the station reported +is taken. For stations beyond -60˚ and 60˚ latitude, the values are +station reported values in every instance. See + +for more detail on the correction methods. + +## WMO Resolution 40. NOAA Policy + +*Users of these data should take into account the following (from the +[NCDC website](http://www7.ncdc.noaa.gov/CDO/cdoselect.cmd?datasetabbv=GSOD&countryabbv=&georegionabbv=)):* + +> "The following data and products may have conditions placed on their +international commercial use. They can be used within the U.S. or for +non-commercial international activities without restriction. The +non-U.S. data cannot be redistributed for commercial purposes. +Re-distribution of these data by others must provide this same +notification." [WMO Resolution 40. NOAA +Policy](http://www.wmo.int/pages/about/Resolution40.html) diff --git a/vignettes/Specified_stations_for_a_range_of_years.html b/vignettes/Specified_stations_for_a_range_of_years.html index f9c3ae9f..5867d9f6 100644 --- a/vignettes/Specified_stations_for_a_range_of_years.html +++ b/vignettes/Specified_stations_for_a_range_of_years.html @@ -335,6 +335,52 @@

Another Option

write.csv(loop_data, file = "Loop_Survey_Weather_1960-1969") +

Notes

+ +

Sources

+ +

CHELSA climate layers

+ +

CHELSA (climatic surfaces at 1 km resolution) is based on a quasi-mechanistical +statistical downscaling of the ERA interim global circulation model +(Karger et al. 2016). ESA's CCI-LC cloud probability monthly averages are based +on the MODIS snow products (MOD10A2). http://chelsa-climate.org/

+ +

EarthEnv MODIS cloud fraction

+ +

http://www.earthenv.org/cloud

+ +

ESA's CCI-LC cloud probability

+ +

http://maps.elie.ucl.ac.be/CCI/viewer/index.php

+ +

Elevation Values

+ +

90m hole-filled SRTM digital elevation (Jarvis et al. 2008) was used +to identify and correct/remove elevation errors in data for station +locations between -60˚ and 60˚ latitude. This applies to cases here +where elevation was missing in the reported values as well. In case the +station reported an elevation and the DEM does not, the station reported +is taken. For stations beyond -60˚ and 60˚ latitude, the values are +station reported values in every instance. See +https://github.com/adamhsparks/GSODR/blob/devel/data-raw/fetch_isd-history.md +for more detail on the correction methods.

+ +

WMO Resolution 40. NOAA Policy

+ +

Users of these data should take into account the following (from the +NCDC website):

+ +
+

“The following data and products may have conditions placed on their +international commercial use. They can be used within the U.S. or for +non-commercial international activities without restriction. The +non-U.S. data cannot be redistributed for commercial purposes. +Re-distribution of these data by others must provide this same +notification.” WMO Resolution 40. NOAA +Policy

+
+ diff --git a/vignettes/Specified_stations_for_a_range_of_years.md b/vignettes/Specified_stations_for_a_range_of_years.md index c100996e..3542b8a1 100644 --- a/vignettes/Specified_stations_for_a_range_of_years.md +++ b/vignettes/Specified_stations_for_a_range_of_years.md @@ -135,4 +135,44 @@ loop_data <- reformat_GSOD(file_list = local_files) write.csv(loop_data, file = "Loop_Survey_Weather_1960-1969") ``` +# Notes + +## Sources + +#### CHELSA climate layers +CHELSA (climatic surfaces at 1 km resolution) is based on a quasi-mechanistical +statistical downscaling of the ERA interim global circulation model +(Karger et al. 2016). ESA's CCI-LC cloud probability monthly averages are based +on the MODIS snow products (MOD10A2). + +#### EarthEnv MODIS cloud fraction + + +#### ESA's CCI-LC cloud probability + + +#### Elevation Values + +90m hole-filled SRTM digital elevation (Jarvis *et al.* 2008) was used +to identify and correct/remove elevation errors in data for station +locations between -60˚ and 60˚ latitude. This applies to cases here +where elevation was missing in the reported values as well. In case the +station reported an elevation and the DEM does not, the station reported +is taken. For stations beyond -60˚ and 60˚ latitude, the values are +station reported values in every instance. See + +for more detail on the correction methods. + +## WMO Resolution 40. NOAA Policy + +*Users of these data should take into account the following (from the +[NCDC website](http://www7.ncdc.noaa.gov/CDO/cdoselect.cmd?datasetabbv=GSOD&countryabbv=&georegionabbv=)):* + +> "The following data and products may have conditions placed on their +international commercial use. They can be used within the U.S. or for +non-commercial international activities without restriction. The +non-U.S. data cannot be redistributed for commercial purposes. +Re-distribution of these data by others must provide this same +notification." [WMO Resolution 40. NOAA +Policy](http://www.wmo.int/pages/about/Resolution40.html) diff --git a/vignettes/Working_with_spatial_and_climate_data.R b/vignettes/Working_with_spatial_and_climate_data.R new file mode 100644 index 00000000..d9ac29b5 --- /dev/null +++ b/vignettes/Working_with_spatial_and_climate_data.R @@ -0,0 +1,63 @@ +## ----example_1, eval=TRUE, message=FALSE, results='hide'----------------- +library(GSODR) +get_GSOD(years = 2010, country = "Philippines", dsn = "~/", + filename = "PHL-2010", GPKG = TRUE, max_missing = 5) + +## ----example_1.1, eval=TRUE, echo=TRUE----------------------------------- +library(rgdal) +library(spacetime) +library(plotKML) + +layers <- ogrListLayers(dsn = path.expand("~/PHL-2010.gpkg")) +pnts <- readOGR(dsn = path.expand("~/PHL-2010.gpkg"), layers[1]) + +# Plot results in Google Earth as a spacetime object: +pnts$DATE = as.Date(paste(pnts$YEAR, pnts$MONTH, pnts$DAY, sep = "-")) +row.names(pnts) <- paste("point", 1:nrow(pnts), sep = "") + +tmp_ST <- STIDF(sp = as(pnts, "SpatialPoints"), + time = pnts$DATE - 0.5, + data = pnts@data[, c("TEMP", "STNID")], + endTime = pnts$DATE + 0.5) + +shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png" + +kml(tmp_ST, dtime = 24 * 3600, colour = TEMP, shape = shape, labels = TEMP, + file.name = "Temperatures_PHL_2010-2010.kml", folder.name = "TEMP") + +system("zip -m Temperatures_PHL_2010-2010.kmz Temperatures_PHL_2010-2010.kml") + +## ----example_1.2, eval=TRUE, fig.cap="Comparison of GSOD daily values and average monthly values with CHELSA climate monthly values", fig.width=8---- +library(dplyr) +library(ggplot2) +library(reshape2) + +data(GSOD_clim) +cnames <- paste0("CHELSA_temp_", 1:12, "_1979-2013") +clim_temp <- GSOD_clim[GSOD_clim$STNID %in% pnts$STNID, + paste(c("STNID", cnames))] +clim_temp_df <- data.frame(STNID = rep(clim_temp$STNID, 12), + MONTHC = as.vector(sapply(1:12, rep, + times = nrow(clim_temp))), + CHELSA_TEMP = as.vector(unlist(clim_temp[, cnames]))) + +pnts$MONTHC <- as.numeric(paste(pnts$MONTH)) +temp <- left_join(pnts@data, clim_temp_df, by = c("STNID", "MONTHC")) + +temp <- temp %>% + group_by(MONTH) %>% + mutate(AVG_DAILY_TEMP = round(mean(TEMP), 1)) + +df_melt <- na.omit(melt(temp[, c("STNID", "DATE", "CHELSA_TEMP", "TEMP", "AVG_DAILY_TEMP")], + id = c("DATE", "STNID"))) + +ggplot(df_melt, aes(x = DATE, y = value)) + + geom_point(aes(color = variable), alpha = 0.2) + + scale_x_date(date_labels = "%b") + + ylab("Temperature (C)") + + xlab("Month") + + labs(colour = "") + + scale_color_brewer(palette = "Dark2") + + facet_wrap( ~ STNID) + + diff --git a/vignettes/Working_with_spatial_and_climate_data.Rmd b/vignettes/Working_with_spatial_and_climate_data.Rmd index a0363acb..ff9eca64 100644 --- a/vignettes/Working_with_spatial_and_climate_data.Rmd +++ b/vignettes/Working_with_spatial_and_climate_data.Rmd @@ -22,20 +22,58 @@ CHELSA temperatures (1979-2013). Download data for Philippines for year 2010 and generate a spatial, year summary file, PHL-2010.gpkg, in the user's home directory. -```{r example_1, eval=TRUE, message=FALSE, results='hide'} + +```r library(GSODR) get_GSOD(years = 2010, country = "Philippines", dsn = "~/", filename = "PHL-2010", GPKG = TRUE, max_missing = 5) ``` -```{r example_1.1, eval=TRUE, echo=TRUE} + +```r library(rgdal) +``` + +``` +## Loading required package: sp +``` + +``` +## rgdal: version: 1.2-4, (SVN revision 643) +## Geospatial Data Abstraction Library extensions to R successfully loaded +## Loaded GDAL runtime: GDAL 1.11.5, released 2016/07/01 +## Path to GDAL shared files: /usr/local/Cellar/gdal/1.11.5_1/share/gdal +## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493] +## Path to PROJ.4 shared files: (autodetected) +## Linking to sp version: 1.2-3 +``` + +```r library(spacetime) library(plotKML) +``` + +``` +## plotKML version 0.5-6 (2016-05-02) +``` +``` +## URL: http://plotkml.r-forge.r-project.org/ +``` + +```r layers <- ogrListLayers(dsn = path.expand("~/PHL-2010.gpkg")) pnts <- readOGR(dsn = path.expand("~/PHL-2010.gpkg"), layers[1]) +``` + +``` +## OGR data source with driver: GPKG +## Source: "/Users/U8004755/PHL-2010.gpkg", layer: "GSOD" +## with 4703 features +## It has 46 fields +``` +```r # Plot results in Google Earth as a spacetime object: pnts$DATE = as.Date(paste(pnts$YEAR, pnts$MONTH, pnts$DAY, sep = "-")) row.names(pnts) <- paste("point", 1:nrow(pnts), sep = "") @@ -49,15 +87,50 @@ shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png" kml(tmp_ST, dtime = 24 * 3600, colour = TEMP, shape = shape, labels = TEMP, file.name = "Temperatures_PHL_2010-2010.kml", folder.name = "TEMP") +``` + +``` +## KML file opened for writing... +``` +``` +## Writing to KML... +``` + +``` +## Closing Temperatures_PHL_2010-2010.kml +``` + +```r system("zip -m Temperatures_PHL_2010-2010.kmz Temperatures_PHL_2010-2010.kml") ``` Compare the GSOD weather data from the Philippines with climatic data provided by the GSODR package in the `GSOD_clim` data set. -```{r example_1.2, eval=TRUE, fig.cap="Comparison of GSOD daily values and average monthly values with CHELSA climate monthly values", fig.width=8} + +```r library(dplyr) +``` + +``` +## +## Attaching package: 'dplyr' +``` + +``` +## The following objects are masked from 'package:stats': +## +## filter, lag +``` + +``` +## The following objects are masked from 'package:base': +## +## intersect, setdiff, setequal, union +``` + +```r library(ggplot2) library(reshape2) @@ -72,7 +145,14 @@ clim_temp_df <- data.frame(STNID = rep(clim_temp$STNID, 12), pnts$MONTHC <- as.numeric(paste(pnts$MONTH)) temp <- left_join(pnts@data, clim_temp_df, by = c("STNID", "MONTHC")) +``` + +``` +## Warning in left_join_impl(x, y, by$x, by$y, suffix$x, suffix$y): joining +## factors with different levels, coercing to character vector +``` +```r temp <- temp %>% group_by(MONTH) %>% mutate(AVG_DAILY_TEMP = round(mean(TEMP), 1)) @@ -88,15 +168,30 @@ ggplot(df_melt, aes(x = DATE, y = value)) + labs(colour = "") + scale_color_brewer(palette = "Dark2") + facet_wrap( ~ STNID) - ``` +![Comparison of GSOD daily values and average monthly values with CHELSA climate monthly values](figure/example_1.2-1.png) + # Notes -## Elevation Values +## Sources + +#### CHELSA climate layers +CHELSA (climatic surfaces at 1 km resolution) is based on a quasi-mechanistical +statistical downscaling of the ERA interim global circulation model +(Karger et al. 2016). ESA's CCI-LC cloud probability monthly averages are based +on the MODIS snow products (MOD10A2). + +#### EarthEnv MODIS cloud fraction + + +#### ESA's CCI-LC cloud probability + + +#### Elevation Values -90 metre (90m) hole-filled SRTM digital elevation (Jarvis *et al.* 2008) was -used to identify and correct/remove elevation errors in data for station +90m hole-filled SRTM digital elevation (Jarvis *et al.* 2008) was used +to identify and correct/remove elevation errors in data for station locations between -60˚ and 60˚ latitude. This applies to cases here where elevation was missing in the reported values as well. In case the station reported an elevation and the DEM does not, the station reported diff --git a/vignettes/Working_with_spatial_and_climate_data.html b/vignettes/Working_with_spatial_and_climate_data.html new file mode 100644 index 00000000..82cccd9b --- /dev/null +++ b/vignettes/Working_with_spatial_and_climate_data.html @@ -0,0 +1,408 @@ + + + + + +Introduction + + + + + + + + + + + + + + + + + + + + +

Introduction

+ +

The GSODR package provides the ability to interact with GSOD data using +spatial methods. The get_GSOD function allows for the data to be saved as +a GeoPackage file which can be read by most GIS software packages or in R using +R's GIS capabilities with contributed packages as well.

+ +

Following is an example of how you might download and save GSOD annual data for +a given country, Philippines in this example, and convert it into a KML file for +viewing in GoogleEarth. The second portion uses the same GeoPackage file to +import the data back into R and combine the GSOD data with the included CHELSA +data and plot the station temperatures for daily GSOD, average monthly GSOD and +CHELSA temperatures (1979-2013).

+ +

Example - Download and plot data for a single country

+ +

Download data for Philippines for year 2010 and generate a spatial, year summary +file, PHL-2010.gpkg, in the user's home directory.

+ +
library(GSODR)
+get_GSOD(years = 2010, country = "Philippines", dsn = "~/",
+         filename = "PHL-2010", GPKG = TRUE, max_missing = 5)
+
+ +
library(rgdal)
+
+ +
## Loading required package: sp
+
+ +
## rgdal: version: 1.2-4, (SVN revision 643)
+##  Geospatial Data Abstraction Library extensions to R successfully loaded
+##  Loaded GDAL runtime: GDAL 1.11.5, released 2016/07/01
+##  Path to GDAL shared files: /usr/local/Cellar/gdal/1.11.5_1/share/gdal
+##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
+##  Path to PROJ.4 shared files: (autodetected)
+##  Linking to sp version: 1.2-3
+
+ +
library(spacetime)
+library(plotKML)
+
+ +
## plotKML version 0.5-6 (2016-05-02)
+
+ +
## URL: http://plotkml.r-forge.r-project.org/
+
+ +
layers <- ogrListLayers(dsn = path.expand("~/PHL-2010.gpkg"))
+pnts <- readOGR(dsn = path.expand("~/PHL-2010.gpkg"), layers[1])
+
+ +
## OGR data source with driver: GPKG 
+## Source: "/Users/U8004755/PHL-2010.gpkg", layer: "GSOD"
+## with 4703 features
+## It has 46 fields
+
+ +
# Plot results in Google Earth as a spacetime object:
+pnts$DATE = as.Date(paste(pnts$YEAR, pnts$MONTH, pnts$DAY, sep = "-"))
+row.names(pnts) <- paste("point", 1:nrow(pnts), sep = "")
+
+tmp_ST <- STIDF(sp = as(pnts, "SpatialPoints"),
+                time = pnts$DATE - 0.5,
+                data = pnts@data[, c("TEMP", "STNID")],
+                endTime = pnts$DATE + 0.5)
+
+shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png"
+
+kml(tmp_ST, dtime = 24 * 3600, colour = TEMP, shape = shape, labels = TEMP,
+    file.name = "Temperatures_PHL_2010-2010.kml", folder.name = "TEMP")
+
+ +
## KML file opened for writing...
+
+ +
## Writing to KML...
+
+ +
## Closing  Temperatures_PHL_2010-2010.kml
+
+ +
system("zip -m Temperatures_PHL_2010-2010.kmz Temperatures_PHL_2010-2010.kml")
+
+ +

Compare the GSOD weather data from the Philippines with climatic data provided +by the GSODR package in the GSOD_clim data set.

+ +
library(dplyr)
+
+ +
## 
+## Attaching package: 'dplyr'
+
+ +
## The following objects are masked from 'package:stats':
+## 
+##     filter, lag
+
+ +
## The following objects are masked from 'package:base':
+## 
+##     intersect, setdiff, setequal, union
+
+ +
library(ggplot2)
+library(reshape2)
+
+data(GSOD_clim)
+cnames <- paste0("CHELSA_temp_", 1:12, "_1979-2013")
+clim_temp <- GSOD_clim[GSOD_clim$STNID %in% pnts$STNID,
+                       paste(c("STNID", cnames))]
+clim_temp_df <- data.frame(STNID = rep(clim_temp$STNID, 12),
+                           MONTHC = as.vector(sapply(1:12, rep,
+                                                    times = nrow(clim_temp))), 
+                           CHELSA_TEMP = as.vector(unlist(clim_temp[, cnames])))
+
+pnts$MONTHC <- as.numeric(paste(pnts$MONTH))
+temp <- left_join(pnts@data, clim_temp_df, by = c("STNID", "MONTHC"))
+
+ +
## Warning in left_join_impl(x, y, by$x, by$y, suffix$x, suffix$y): joining
+## factors with different levels, coercing to character vector
+
+ +
temp <- temp %>% 
+  group_by(MONTH) %>% 
+  mutate(AVG_DAILY_TEMP = round(mean(TEMP), 1))
+
+df_melt <- na.omit(melt(temp[, c("STNID", "DATE", "CHELSA_TEMP", "TEMP", "AVG_DAILY_TEMP")],
+                        id = c("DATE", "STNID")))
+
+ggplot(df_melt, aes(x = DATE, y = value)) +
+  geom_point(aes(color = variable), alpha = 0.2) +
+  scale_x_date(date_labels = "%b") +
+  ylab("Temperature (C)") +
+  xlab("Month") +
+  labs(colour = "") +
+  scale_color_brewer(palette = "Dark2") +
+  facet_wrap( ~ STNID)
+
+ +

Comparison of GSOD daily values and average monthly values with CHELSA climate monthly values

+ +

Notes

+ +

Sources

+ +

CHELSA climate layers

+ +

CHELSA (climatic surfaces at 1 km resolution) is based on a quasi-mechanistical +statistical downscaling of the ERA interim global circulation model +(Karger et al. 2016). ESA's CCI-LC cloud probability monthly averages are based +on the MODIS snow products (MOD10A2). http://chelsa-climate.org/

+ +

EarthEnv MODIS cloud fraction

+ +

http://www.earthenv.org/cloud

+ +

ESA's CCI-LC cloud probability

+ +

http://maps.elie.ucl.ac.be/CCI/viewer/index.php

+ +

Elevation Values

+ +

90m hole-filled SRTM digital elevation (Jarvis et al. 2008) was used +to identify and correct/remove elevation errors in data for station +locations between -60˚ and 60˚ latitude. This applies to cases here +where elevation was missing in the reported values as well. In case the +station reported an elevation and the DEM does not, the station reported +is taken. For stations beyond -60˚ and 60˚ latitude, the values are +station reported values in every instance. See +https://github.com/adamhsparks/GSODR/blob/devel/data-raw/fetch_isd-history.md +for more detail on the correction methods.

+ +

WMO Resolution 40. NOAA Policy

+ +

Users of these data should take into account the following (from the +NCDC website):

+ +
+

“The following data and products may have conditions placed on their +international commercial use. They can be used within the U.S. or for +non-commercial international activities without restriction. The +non-U.S. data cannot be redistributed for commercial purposes. +Re-distribution of these data by others must provide this same +notification.” WMO Resolution 40. NOAA +Policy

+
+ +

References

+ +

Stachelek, J. 2016. Using the Geopackage Format with R. +URL: https://jsta.github.io/2016/07/14/geopackage-r.html

+ + + + diff --git a/vignettes/Working_with_spatial_and_climate_data.md b/vignettes/Working_with_spatial_and_climate_data.md new file mode 100644 index 00000000..ff9eca64 --- /dev/null +++ b/vignettes/Working_with_spatial_and_climate_data.md @@ -0,0 +1,218 @@ + + +# Introduction + +The GSODR package provides the ability to interact with GSOD data using +spatial methods. The `get_GSOD` function allows for the data to be saved as +a GeoPackage file which can be read by most GIS software packages or in R using +R's GIS capabilities with contributed packages as well. + +Following is an example of how you might download and save GSOD annual data for +a given country, Philippines in this example, and convert it into a KML file for +viewing in GoogleEarth. The second portion uses the same GeoPackage file to +import the data back into R and combine the GSOD data with the included CHELSA +data and plot the station temperatures for daily GSOD, average monthly GSOD and +CHELSA temperatures (1979-2013). + +## Example - Download and plot data for a single country +Download data for Philippines for year 2010 and generate a spatial, year summary +file, PHL-2010.gpkg, in the user's home directory. + + +```r +library(GSODR) +get_GSOD(years = 2010, country = "Philippines", dsn = "~/", + filename = "PHL-2010", GPKG = TRUE, max_missing = 5) +``` + + +```r +library(rgdal) +``` + +``` +## Loading required package: sp +``` + +``` +## rgdal: version: 1.2-4, (SVN revision 643) +## Geospatial Data Abstraction Library extensions to R successfully loaded +## Loaded GDAL runtime: GDAL 1.11.5, released 2016/07/01 +## Path to GDAL shared files: /usr/local/Cellar/gdal/1.11.5_1/share/gdal +## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493] +## Path to PROJ.4 shared files: (autodetected) +## Linking to sp version: 1.2-3 +``` + +```r +library(spacetime) +library(plotKML) +``` + +``` +## plotKML version 0.5-6 (2016-05-02) +``` + +``` +## URL: http://plotkml.r-forge.r-project.org/ +``` + +```r +layers <- ogrListLayers(dsn = path.expand("~/PHL-2010.gpkg")) +pnts <- readOGR(dsn = path.expand("~/PHL-2010.gpkg"), layers[1]) +``` + +``` +## OGR data source with driver: GPKG +## Source: "/Users/U8004755/PHL-2010.gpkg", layer: "GSOD" +## with 4703 features +## It has 46 fields +``` + +```r +# Plot results in Google Earth as a spacetime object: +pnts$DATE = as.Date(paste(pnts$YEAR, pnts$MONTH, pnts$DAY, sep = "-")) +row.names(pnts) <- paste("point", 1:nrow(pnts), sep = "") + +tmp_ST <- STIDF(sp = as(pnts, "SpatialPoints"), + time = pnts$DATE - 0.5, + data = pnts@data[, c("TEMP", "STNID")], + endTime = pnts$DATE + 0.5) + +shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png" + +kml(tmp_ST, dtime = 24 * 3600, colour = TEMP, shape = shape, labels = TEMP, + file.name = "Temperatures_PHL_2010-2010.kml", folder.name = "TEMP") +``` + +``` +## KML file opened for writing... +``` + +``` +## Writing to KML... +``` + +``` +## Closing Temperatures_PHL_2010-2010.kml +``` + +```r +system("zip -m Temperatures_PHL_2010-2010.kmz Temperatures_PHL_2010-2010.kml") +``` + +Compare the GSOD weather data from the Philippines with climatic data provided +by the GSODR package in the `GSOD_clim` data set. + + +```r +library(dplyr) +``` + +``` +## +## Attaching package: 'dplyr' +``` + +``` +## The following objects are masked from 'package:stats': +## +## filter, lag +``` + +``` +## The following objects are masked from 'package:base': +## +## intersect, setdiff, setequal, union +``` + +```r +library(ggplot2) +library(reshape2) + +data(GSOD_clim) +cnames <- paste0("CHELSA_temp_", 1:12, "_1979-2013") +clim_temp <- GSOD_clim[GSOD_clim$STNID %in% pnts$STNID, + paste(c("STNID", cnames))] +clim_temp_df <- data.frame(STNID = rep(clim_temp$STNID, 12), + MONTHC = as.vector(sapply(1:12, rep, + times = nrow(clim_temp))), + CHELSA_TEMP = as.vector(unlist(clim_temp[, cnames]))) + +pnts$MONTHC <- as.numeric(paste(pnts$MONTH)) +temp <- left_join(pnts@data, clim_temp_df, by = c("STNID", "MONTHC")) +``` + +``` +## Warning in left_join_impl(x, y, by$x, by$y, suffix$x, suffix$y): joining +## factors with different levels, coercing to character vector +``` + +```r +temp <- temp %>% + group_by(MONTH) %>% + mutate(AVG_DAILY_TEMP = round(mean(TEMP), 1)) + +df_melt <- na.omit(melt(temp[, c("STNID", "DATE", "CHELSA_TEMP", "TEMP", "AVG_DAILY_TEMP")], + id = c("DATE", "STNID"))) + +ggplot(df_melt, aes(x = DATE, y = value)) + + geom_point(aes(color = variable), alpha = 0.2) + + scale_x_date(date_labels = "%b") + + ylab("Temperature (C)") + + xlab("Month") + + labs(colour = "") + + scale_color_brewer(palette = "Dark2") + + facet_wrap( ~ STNID) +``` + +![Comparison of GSOD daily values and average monthly values with CHELSA climate monthly values](figure/example_1.2-1.png) + +# Notes + +## Sources + +#### CHELSA climate layers +CHELSA (climatic surfaces at 1 km resolution) is based on a quasi-mechanistical +statistical downscaling of the ERA interim global circulation model +(Karger et al. 2016). ESA's CCI-LC cloud probability monthly averages are based +on the MODIS snow products (MOD10A2). + +#### EarthEnv MODIS cloud fraction + + +#### ESA's CCI-LC cloud probability + + +#### Elevation Values + +90m hole-filled SRTM digital elevation (Jarvis *et al.* 2008) was used +to identify and correct/remove elevation errors in data for station +locations between -60˚ and 60˚ latitude. This applies to cases here +where elevation was missing in the reported values as well. In case the +station reported an elevation and the DEM does not, the station reported +is taken. For stations beyond -60˚ and 60˚ latitude, the values are +station reported values in every instance. See + +for more detail on the correction methods. + +## WMO Resolution 40. NOAA Policy + +*Users of these data should take into account the following (from the +[NCDC website](http://www7.ncdc.noaa.gov/CDO/cdoselect.cmd?datasetabbv=GSOD&countryabbv=&georegionabbv=)):* + +> "The following data and products may have conditions placed on their +international commercial use. They can be used within the U.S. or for +non-commercial international activities without restriction. The +non-U.S. data cannot be redistributed for commercial purposes. +Re-distribution of these data by others must provide this same +notification." [WMO Resolution 40. NOAA +Policy](http://www.wmo.int/pages/about/Resolution40.html) + +# References +Stachelek, J. 2016. Using the Geopackage Format with R. +URL: https://jsta.github.io/2016/07/14/geopackage-r.html diff --git a/vignettes/figure/example_1.2-1.png b/vignettes/figure/example_1.2-1.png new file mode 100644 index 00000000..eac980c8 Binary files /dev/null and b/vignettes/figure/example_1.2-1.png differ