Merge pull request #74 from bczernecki/dev

climate 1.0.4

.Rbuildignore

@@ -9,5 +9,5 @@
 ^\.httr-oauth$
 ^data-raw$
 ^vignettes/articles$
-
 ^\.github$
+^\.Rhistory$

.github/workflows/R-CMD-check.yaml

@@ -31,9 +31,9 @@ jobs:
       matrix:
         config:
           - {os: macOS-latest,   r: 'release'}
-          - {os: windows-latest, r: '3.6'}
-          - {os: ubuntu-16.04,   r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"}
-          - {os: ubuntu-16.04,   r: '3.5',     rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"}
+          - {os: windows-latest, r: 'release'}
+          - {os: ubuntu-18.04,   r: 'release'}
+          - {os: ubuntu-18.04,   r: '3.6'}
     env:
       R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
       RSPM: ${{ matrix.config.rspm }}
@@ -67,10 +67,14 @@ jobs:
       - name: Install system dependencies
         if: runner.os == 'Linux'
         run: |
+          ls
+          pwd
+          echo -e 'LC_CTYPE = "en_US.UTF-8"\nLC_NUMERIC="C"\nLC_TIME="en_US.UTF-8"\nLC_COLLATE="C"\nLC_MONETARY="en_US.UTF-8"\nLC_MESSAGES="en_US.UTF-8"\nLC_PAPER="en_US.UTF-8"\nLC_NAME="C"\nLC_ADDRESS="C"\nLC_TELEPHONE="C"\nLC_MEASUREMENT="en_US.UTF-8"\nLC_IDENTIFICATION="C"' > .Renviron
+
           while read -r cmd
           do
             eval sudo $cmd
-          done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "16.04"))')
+          done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "18.04"))')
 
       - name: Install dependencies
         run: |

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: climate
 Title: Interface to Download Meteorological (and Hydrological) Datasets
-Version: 1.0.3
+Version: 1.0.4
 Authors@R: c(person(given = "Bartosz",
            family = "Czernecki",
            role = c("aut", "cre"),
@@ -27,7 +27,7 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.2
 Depends: 
     R (>= 3.1)
 Imports: 

NAMESPACE

@@ -18,6 +18,7 @@ export(nearest_stations_ogimet)
 export(ogimet_daily)
 export(ogimet_hourly)
 export(sounding_wyoming)
+export(spheroid_dist)
 export(stations_ogimet)
 export(test_url)
 import(httr)

NEWS.md

@@ -1,3 +1,10 @@
+# climate 1.0.4
+
+* Function `spheroid_dist` added to improve accuracy of calculations between points, but also avoid installing GIS dependencies (thanks to @kadyb)
+* Function `nearest_stations_imgw()` now uses the Vincenty's formula in `spheroid_dist` to calculate the distance between points on a spheroid, not the Euclidean distance (previously results were inaccurate for some specific cases)
+* minor bugs fixes and improvements
+
+
 # climate 1.0.3
 
 * Adding possibility to download BUFR vertical sounding dataset from `http://weather.uwyo.edu/upperair/sounding.html`; extra information with supporting example added to the `sounding_wyoming`'s documentation

R/check_locale.R

@@ -6,16 +6,14 @@
 
 check_locale = function(){
 
-  if(Sys.getlocale("LC_CTYPE") %in% c("C.UTF-8", "en_US.iso885915")){
+  if (Sys.getlocale("LC_CTYPE") %in% c("C.UTF-8", "en_US.iso885915")) {
     locale = Sys.getlocale("LC_CTYPE")
     message(paste0("    Your system locale is: " , locale," which may cause trouble. 
     Please consider changing it manually while working with climate, e.g.: 
     Sys.setlocale(category = 'LC_ALL', locale = 'en_US.UTF-8') "))
     Sys.sleep(4)
     return(1)
-
   } else {
-
     return(0)
   }
 

R/meteo_imgw_daily.R

@@ -91,15 +91,15 @@ meteo_imgw_daily = function(rank = "synop", year, status = FALSE, coords = FALSE
         colnames(data1) = meta[[1]]$parameters
 
         file2 = paste(temp2, dir(temp2), sep = "/")[2]
-        if(translit){
+        if (translit) {
           data2 = as.data.frame(data.table::fread(cmd = paste("iconv -f CP1250 -t ASCII//TRANSLIT", file2)))
         } else {
           data2 = read.csv(file2, header = FALSE, stringsAsFactors = FALSE, fileEncoding = "CP1250")
         }
         colnames(data2) = meta[[2]]$parameters
 
         # usuwa statusy
-        if(status == FALSE){
+        if (status == FALSE) {
           data1[grep("^Status", colnames(data1))] = NULL
           data2[grep("^Status", colnames(data2))] = NULL
         }
@@ -123,7 +123,7 @@ meteo_imgw_daily = function(rank = "synop", year, status = FALSE, coords = FALSE
 
     ######################
     ###### KLIMAT: #######
-    if(rank == "climate") {
+    if (rank == "climate") {
       address = paste0(base_url, "dane_meteorologiczne/dobowe/klimat",
                         "/", catalog, "/")
       #folder_contents = getURL(address, ftp.use.epsv = FALSE, dirlistonly = FALSE) # zawartosc folderu dla wybranego yearu
@@ -138,29 +138,29 @@ meteo_imgw_daily = function(rank = "synop", year, status = FALSE, coords = FALSE
       # w tym miejscu trzeba przemyslec fragment kodu do dodania dla pojedynczej stacji jesli tak sobie zazyczy uzytkownik:
       # na podstawie zawartosci obiektu files
 
-      for(j in seq_along(addresses_to_download)){
+      for (j in seq_along(addresses_to_download)) {
         temp = tempfile()
         temp2 = tempfile()
         test_url(addresses_to_download[j], temp)
         unzip(zipfile = temp, exdir = temp2)
         file1 = paste(temp2, dir(temp2), sep = "/")[1]
-        if(translit){
+        if (translit) {
           data1 = as.data.frame(data.table::fread(cmd = paste("iconv -f CP1250 -t ASCII//TRANSLIT", file1)))
         } else {
           data1 = read.csv(file1, header = FALSE, stringsAsFactors = FALSE, fileEncoding = "CP1250")
         }
         colnames(data1) = meta[[1]]$parameters
 
         file2 = paste(temp2, dir(temp2), sep = "/")[2]
-        if(translit){
+        if (translit) {
           data2 = as.data.frame(data.table::fread(cmd = paste("iconv -f CP1250 -t ASCII//TRANSLIT", file2)))
         } else {
           data2 = read.csv(file2, header = FALSE, stringsAsFactors = FALSE, fileEncoding = "CP1250")
         }
         colnames(data2) = meta[[2]]$parameters
 
         # usuwa statusy
-        if(status == FALSE){
+        if (status == FALSE) {
           data1[grep("^Status", colnames(data1))] = NULL
           data2[grep("^Status", colnames(data2))] = NULL
         }
@@ -176,7 +176,7 @@ meteo_imgw_daily = function(rank = "synop", year, status = FALSE, coords = FALSE
 
     ######################
     ######## OPAD: #######
-    if(rank == "precip") {
+    if (rank == "precip") {
       address = paste0(base_url, "dane_meteorologiczne/dobowe/opad",
                         "/", catalog, "/")
       #folder_contents = getURL(address, ftp.use.epsv = FALSE, dirlistonly = FALSE) # zawartosc folderu dla wybranego yearu
@@ -189,21 +189,21 @@ meteo_imgw_daily = function(rank = "synop", year, status = FALSE, coords = FALSE
       files = as.character(readHTMLTable(folder_contents)[[1]]$Name[ind])
       addresses_to_download = paste0(address, files)
 
-      for(j in seq_along(addresses_to_download)){
+      for (j in seq_along(addresses_to_download)) {
         temp = tempfile()
         temp2 = tempfile()
         test_url(addresses_to_download[j], temp)
         unzip(zipfile = temp, exdir = temp2)
         file1 = paste(temp2, dir(temp2), sep = "/")[1]
-        if(translit){
+        if (translit) {
           data1 = as.data.frame(data.table::fread(cmd = paste("iconv -f CP1250 -t ASCII//TRANSLIT", file1)))
         } else {
           data1 = read.csv(file1, header = FALSE, stringsAsFactors = FALSE, fileEncoding = "CP1250")
         }
 
         colnames(data1) = meta[[1]]$parameters
         # usuwa statusy
-        if(status == FALSE){
+        if (status == FALSE) {
           data1[grep("^Status", colnames(data1))] = NULL
         }
 
@@ -250,7 +250,7 @@ meteo_imgw_daily = function(rank = "synop", year, status = FALSE, coords = FALSE
   }
 
   # sortowanie w zaleznosci od nazw kolumn - raz jest "kod stacji", raz "id"
-  if(sum(grepl(x = colnames(all_data), pattern = "Kod stacji"))){
+  if (sum(grepl(x = colnames(all_data), pattern = "Kod stacji"))) {
     all_data = all_data[order(all_data$`Kod stacji`, all_data$Rok, all_data$Miesiac, all_data$Dzien), ]
   } else {
     all_data = all_data[order(all_data$id, all_data$Rok, all_data$Miesiac, all_data$Dzien), ]

R/meteo_imgw_monthly.R

@@ -26,8 +26,8 @@
 #'   
 #'   # please note that station names may change over time 
 #'   # and thus 2 names are required in some cases:
-#'   df = meteo_imgw_monthly(rank = 'synop', year = 1991:2000, 
-#'           coords = TRUE, station = c("POZNAŃ","POZNAŃ-ŁAWICA")) 
+#'   # df = meteo_imgw_monthly(rank = 'synop', year = 1991:2000, 
+#'   #         coords = TRUE, station = c("POZNAŃ","POZNAŃ-ŁAWICA")) 
 #' }
 #'
 meteo_imgw_monthly = function(rank = "synop", year, status = FALSE, coords = FALSE, station = NULL, col_names = "short", ...){

R/nearest_stations_imgw.R

@@ -3,89 +3,102 @@
 #' Returns a data frame of meteorological or hydrological stations with their coordinates in particular year. 
 #' The returned object is valid only for a given year and type of stations (e.g. "synop", "climate" or "precip"). If `add_map = TRUE` additional map of downloaded data is added. 
 #'
-#' @param type data name;"meteo" (default), "hydro" 
-#' @param rank rank of the stations: "synop" (default), "climate", or "precip"; Only valid if type = "meteo
-#' @param year select year for serching nearest station
+#' @param type data name; "meteo" (default), "hydro" 
+#' @param rank rank of the stations: "synop" (default), "climate", or "precip"; Only valid if type = "meteo"
+#' @param year select year for searching nearest station
 #' @param add_map logical - whether to draw a map for a returned data frame (requires maps/mapdata packages)
 #' @param point a vector of two coordinates (longitude, latitude) for a point we want to find nearest stations to (e.g. c(15, 53)); If not provided calculated as a mean longitude and latitude for the entire dataset
 #' @param no_of_stations how many nearest stations will be returned from the given geographical coordinates. 50 used by default
 #' @param ... extra arguments to be provided to the [graphics::plot()] function (only if add_map = TRUE)
-#' @importFrom XML readHTMLTable
 #' @export
-#' @return A data.frame with a list of nearest stations. Each row represents metadata for station which collected measurements in a given year. Particular columns contain stations metadata (e.g. station ID, geographical coordinates, official name, distance from a given coordinates). 
+#' @return A data.frame with a list of nearest stations. Each row represents metadata for station which collected measurements in a given year. Particular columns contain stations metadata (e.g. station ID, geographical coordinates, official name, distance in kilometers from a given coordinates). 
 #'  
 #' @examples 
 #' \donttest{
 #'   nearest_stations_imgw(type = "hydro", 
-#'   rank="synop",
-#'   year=2018,
+#'   rank = "synop",
+#'   year = 2018,
 #'   point = c(17, 52),
 #'   add_map = TRUE, 
 #'   no_of_stations = 4)
 #' }
 #'
 
-nearest_stations_imgw = function(type = "meteo", 
-                                  rank = "synop",
-                                  year = 2018,
-                                  add_map = TRUE, 
-                                  point = NULL, 
-                                  no_of_stations = 50, ...){
-  if (length(point)>2) {
-    stop("Too many points for the distance calculations. Please provide just one pair of coordinates (e.g. point = c(17,53))")
-  } else if (length(point)<2 | length(point) == 0) {
-    message("
-            The point argument should have two coordinates.
-            We will provide nearest stations for mean location of all available stations.
-            To change it please change the `point` argument c(LON,LAT)" )
-    Sys.sleep(2)
+nearest_stations_imgw = function(type = "meteo",
+                                 rank = "synop",
+                                 year = 2018,
+                                 add_map = TRUE,
+                                 point = NULL,
+                                 no_of_stations = 50,
+                                 ...) {
+  if (length(point) > 2) {
+    stop(paste("Too many points for the distance calculations.",
+               "Please provide just one pair of coordinates (e.g. point = c(17,53))"))
+  } else if (length(point) < 2) {
+    msg = paste("The point argument should have two coordinates.",
+                 "We will provide nearest stations for mean location of all available stations.",
+                 "To change it please change the `point` argument c(LON,LAT)",
+                 sep = "\n")
+    message(msg)
+    Sys.sleep(2) # display message for 2 sec
   }
 
-
-  if (max(year)>=as.integer(substr(Sys.Date(),1,4))-1) {
+  if (!is.null(point) && point[1] > 180) {
+    stop("x should be longitude")
+  } else if (!is.null(point) && point[2] > 90) {
+    stop("y should be latitude")
+  }
+
+  if (max(year) >= as.integer(substr(Sys.Date(), 1, 4)) - 1) {
     message("Data cannot be provided for this repository. Please check the available records at: \n
             https://danepubliczne.imgw.pl/data/dane_pomiarowo_obserwacyjne/")
+    Sys.sleep(2)
   }
 
-  if (type == "meteo"){
-    result = unique(meteo_imgw_monthly(rank = rank, year = year, coords = T)[,c(2:5)])
-  } else if (type=="hydro"){
-    result = unique(hydro_imgw_annual(year = year, coords = T)[,c(1:4)])
+  if (type == "meteo") {
+    result = unique(meteo_imgw_monthly(rank = rank, year = year, coords = TRUE)[, c(2:5)])
+  } else if (type == "hydro") {
+    result = unique(hydro_imgw_annual(year = year, coords = TRUE)[, c(1:4)])
   } else {
     stop("You've provided wrong type argument; please use: \"meteo\", or \"hydro\"")
   }
 
-  if (dim(result)[1]==0) {
+  if (dim(result)[1] == 0) {
     stop("Propobly there is no data in the downloaded object. Please check available records:  
         https://danepubliczne.imgw.pl/data/dane_pomiarowo_obserwacyjne/")
-  } 
-  if (is.null(point)){
-  # workaround for different column names:
-    if(any(colnames(result)=="LON")) point = c(round(mean(result$LON,na.rm=T),2),round(mean(result$LAT,na.rm=T),2))
-    if(any(colnames(result)=="X")) point = c(round(mean(result$X,na.rm=T),2),round(mean(result$Y,na.rm=T),2))
   }
 
+  if (is.null(point)) {
+  # workaround for different column names:
+    if ("LON" %in% colnames(result))
+      point = c(round(mean(result$LON, na.rm = TRUE), 2),
+                round(mean(result$LAT, na.rm = TRUE), 2))
+    if ("X" %in% colnames(result))
+      point = c(round(mean(result$X, na.rm = TRUE), 2),
+                round(mean(result$Y, na.rm = TRUE), 2))
+  }
 
-  point = as.data.frame(t(point))
-  names(point) = c("X", "Y")
-  distmatrix = rbind(point,result[, 2:3])
-  distance_points = stats::dist(distmatrix, method = "euclidean")[1:dim(result)[1]]
-  result["distance [km]"] = round(distance_points * 112.196672, 3)
-  orderd_distance = result[order(result$distance), ]
-  result = orderd_distance[1:no_of_stations, ]
+  dist_vec = double(nrow(result))
+  for (i in seq_along(dist_vec)) {
+    dist_vec[i] = spheroid_dist(point, c(result$X[i], result$Y[i]))
+  }
 
-  message('Currently "climate" package use approximate distance calculation\nIn order to get accurate results please use GIS-based solutions')
+  result["distance"] = round(dist_vec, 3)
+  result = result[order(result$distance), ]
+  result = result[1:no_of_stations, ]
 
+
   # removing rows with all NA records from the obtained dataset;
   # otherwise there might be problems with plotting infinite xlim, ylim, etc..
-  result = result[!apply(is.na(result), 1, sum) == ncol(result),]
+  result = result[!apply(is.na(result), 1, sum) == ncol(result), ]
 
-  if(add_map == TRUE){
-    if (!requireNamespace("maps", quietly = TRUE)){
+  if (add_map == TRUE) {
+    if (!requireNamespace("maps", quietly = TRUE)) {
       stop("package maps required, please install it first")
     }
     # plot a little bit more:
-    addfactor = as.numeric(diff(stats::quantile(result$Y, na.rm = TRUE, c(0.48, 0.51)))) #lat Y
+    addfactor = stats::quantile(result$Y, c(0.48, 0.51), na.rm = TRUE) #lat Y
+    addfactor = as.numeric(diff(addfactor))
     addfactor = ifelse(addfactor > 0.2, 0.2, addfactor)
     addfactor = ifelse(addfactor < 0.05, 0.05, addfactor)
 
@@ -96,8 +109,10 @@ nearest_stations_imgw = function(type = "meteo",
       pch = 19,
       xlab = "longitude",
       ylab = "latitude",
-      xlim = (c(min(c(result$X, point$X), na.rm = T) - 1, max(c(result$X, point$X),na.rm = T) + 1)),
-      ylim = (c(min(c(result$Y, point$Y), na.rm = T) - 1, max(c(result$Y, point$Y),na.rm = T) + 1)),
+      xlim = c(min(c(result[["X"]], point[1]), na.rm = TRUE) - 1,
+               max(c(result[["X"]], point[1]), na.rm = TRUE) + 1),
+      ylim = c(min(c(result[["Y"]], point[2]), na.rm = TRUE) - 1,
+               max(c(result[["Y"]], point[2]), na.rm = TRUE) + 1),
       ...
     )
     graphics::points(
@@ -117,9 +132,10 @@ nearest_stations_imgw = function(type = "meteo",
     maps::map(add = TRUE)
 
   }
-  if (length(year)>1) {
-    message("Please provide only one year. For more years station's metadata may change (name, location or station may stop collecting data)")
+
+  if (length(year) > 1) {
+    message(paste0("Please provide only one year. For more years station's metadata",
+                   "may change (name, location or station may stop collecting data)"))
   }
   return(result)
 }
-