/
rxtracto_3D.R
executable file
·206 lines (184 loc) · 8.6 KB
/
rxtracto_3D.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
#' Extract environmental data in a 3-dimensional box from an 'ERDDAP' server using 'rerddap'.
#'
#' \code{rxtracto_3D} uses the R program 'rerddap' to extract environmental data
#' from an 'ERDDAP' server in an (x,y,z, time) bounding box.
#' The same call could be made directly in rerddap,
#' but function is maintained as it is used in the polygon routine.
#' @export
#' @param dataInfo - the return from an 'rerddap:info' call to an 'ERDDAP' server
#' @param parameter - character string containing the name of the parameter to extract
#' @param xcoord - a real array with the x-coordinates of the trajectory (if longitude in #' decimal degrees East, either 0-360 or -180 to 180)
#' @param ycoord - a real array with the y-coordinate of the trajectory (if latitude in
#' decimal degrees N; -90 to 90)
#' @param zcoord - a real array with the z-coordinate (usually altitude or depth)
#' @param tcoord - a character array with the times of the trajectory in
#' "YYYY-MM-DD" - for now restricted to be time.
#' @param xName - character string with name of the xcoord in the 'ERDDAP' dataset (default "longitude")
#' @param yName - character string with name of the ycoord in the 'ERDDAP' dataset (default "latitude")
#' @param zName - character string with name of the zcoord in the 'ERDDAP' dataset (default "altitude")
#' @param tName - character string with name of the tcoord in the 'ERDDAP' dataset (default "time")
#' @param verbose - logical variable (default FALSE) if the the URL request should be verbose
#' @param cache_remove - logical variable (default TRUE) whether to delete 'rerddap' cache
#' @return If successful a structure with data and dimensions:
#' \itemize{
#' \item extract$data - the data array dimensioned (lon,lat,time)
#' \item extract$varname - the name of the parameter extracted
#' \item extract$datasetname - ERDDAP dataset name
#' \item extract$longitude - the longitudes on some scale as request
#' \item extract$latitude - the latitudes always going south to north
#' \item extract$time - the times of the extracts
#' }
#' else an error string
#' @examples
#' ## toy example to show use
#' ## and keep execution time low
#' ##
#' # dataInfo <- rerddap::info('erdHadISST')
#' parameter <- 'sst'
#' xcoord <- c(-130.5, -130.5)
#' ycoord <- c(40.5, 40.5)
#' tcoord <- c('2006-01-16', '2006-01-16')
#' # extract <- rxtracto_3D(dataInfo, parameter, xcoord = xcoord, ycoord = ycoord,
#' # tcoord = tcoord)
#'
#' ## bathymetry example
#' ## 2-D example getting bathymetry
#' dataInfo <- rerddap::info('etopo360')
#' parameter <- 'altitude'
#' # extract <- rxtracto_3D(dataInfo, parameter, xcoord = xcoord, ycoord = ycoord)
rxtracto_3D <- function(dataInfo, parameter = NULL, xcoord = NULL,
ycoord = NULL, zcoord = NULL, tcoord = NULL,
xName = 'longitude', yName = 'latitude',
zName = 'altitude', tName = 'time',
verbose=FALSE, cache_remove = TRUE) {
# Check Passed Info -------------------------------------------------------
rerddap::cache_setup(temp_dir = TRUE)
callDims <- list(xcoord, ycoord, zcoord, tcoord)
names(callDims) <- c(xName, yName, zName, tName)
dataInfo1 <- dataInfo
urlbase <- dataInfo1$base_url
urlbase <- checkInput(dataInfo1, parameter, urlbase, callDims)
if (is.numeric(urlbase)) {
if (urlbase == -999) {
return("error in inputs")
} else {
return('url is not a valid erddap server')
}
}
# Check and readjust coordinate variables ---------------------------------
# get the actual coordinate values for the dataset
allCoords <- dimvars(dataInfo1)
dataCoordList <- getfileCoords(attr(dataInfo1, "datasetid"), allCoords, urlbase)
if (is.numeric(dataCoordList) ) {
return("Error retrieving coordinate variable")
}
working_coords <- remapCoords(dataInfo1, callDims, dataCoordList, urlbase)
dataInfo1 <- working_coords$dataInfo1
cross_dateline_180 <- working_coords$cross_dateline_180
# Check request is within dataset bounds ----------------------------------
#get limits over new coordinates
xcoordLim <- working_coords$xcoord1
#if (working_coords$latSouth) {
# ycoordLim <- c(min(working_coords$ycoord1), max(working_coords$ycoord1))
#} else {
# ycoordLim <- c(max(working_coords$ycoord1), min(working_coords$ycoord1))
#}
ycoordLim <- c(min(working_coords$ycoord1), max(working_coords$ycoord1))
zcoordLim <- NULL
if (!is.null(working_coords$zcoord1)) {
zcoordLim <- working_coords$zcoord1
if (length(zcoordLim) == 1) {
zcoordLim <- c(zcoordLim, zcoordLim)
}
}
tcoordLim <- NULL
if (!is.null(working_coords$tcoord1)) {
# check for last in time, and convert
isoTime <- dataCoordList$time
udtTime <- parsedate::parse_iso_8601(isoTime)
tcoord1 <- removeLast(isoTime, working_coords$tcoord1)
tcoord1 <- parsedate::parse_iso_8601(tcoord1)
tcoordLim <- tcoord1
}
dimargs <- list(xcoordLim, ycoordLim, zcoordLim, tcoordLim)
names(dimargs) <- c(xName, yName, zName, tName)
dimargs <- Filter(Negate(is.null), dimargs)
#check that coordinate bounds are contained in the dataset
bound_check <- checkBounds(dataCoordList, dimargs, cross_dateline_180)
if (bound_check != 0){
return( 'error in given bounds')
}
# Find dataset coordinates closest to requested coordinates ---------------
erddapList <- findERDDAPcoord(dataCoordList, isoTime, udtTime,
xcoordLim, ycoordLim, tcoordLim, zcoordLim,
xName, yName, tName, zName, cross_dateline_180)
erddapCoords <- erddapList$erddapCoords
# Construct the griddap() command from the input --------------------------
if (cross_dateline_180) {
# upper_bound <- round(max(dataCoordList$longitude), 3)
upper_bound <- max(dataCoordList$longitude) - 0.0001
xcoord_temp <- c(erddapCoords$erddapXcoord[1], upper_bound)
extract1 <- data_extract_read(dataInfo1, callDims, urlbase,
xName, yName, zName, tName, parameter,
xcoord_temp, erddapCoords$erddapYcoord,
erddapCoords$erddapTcoord, erddapCoords$erddapZcoord,
verbose, cache_remove )
if (!is.list(extract1)) {
text1 <- "There was an error in the url call, perhaps a time out."
text2 <- "See message on screen and URL called"
print(paste(text1, text2))
return("URL cal error")
}
# lower_bound <- round(min(dataCoordList$longitude), 3)
lower_bound <- min(dataCoordList$longitude) + 0.0001
xcoord_temp <- c(lower_bound, erddapCoords$erddapXcoord[2])
extract2 <- data_extract_read(dataInfo1, callDims, urlbase,
xName, yName, zName, tName, parameter,
xcoord_temp, erddapCoords$erddapYcoord,
erddapCoords$erddapTcoord, erddapCoords$erddapZcoord,
verbose, cache_remove )
if (!is.list(extract2)) {
text1 <- "There was an error in the url call, perhaps a time out."
text2 <- "See message on screen and URL called"
print(paste(text1, text2))
return("URL cal error")
}
extract2$longitude = make360(extract2$longitude)
# extract <- list(NA, NA, NA, NA, NA, NA)
extract <- vector("list", 6)
lat_len <- length(extract1$latitude)
time_len <- length(extract1$time)
lon_len <- length(extract1$longitude) + length(extract2$longitude)
temp_array <- array(NA_real_, dim = c(lon_len, lat_len, time_len))
temp_array[1:length(extract1$longitude), ,] <- extract1[[1]]
temp_array[(length(extract1$longitude) + 1):lon_len, ,] <- extract2[[1]]
names(extract) <- names(extract1)
extract[[1]] <- temp_array
extract$datasetname <- extract1$datasetname
extract$latitude <- extract1$latitude
extract$altitude <- ifelse(is.null(extract1$altitude), NA, extract1$altitude)
if (is.null(extract1$time)) {
extract$time <- NA
} else{
extract$time <- extract1$time
}
extract$longitude <- c(extract1$longitude, extract2$longitude)
}else {
extract <- data_extract_read(dataInfo1, callDims, urlbase,
xName, yName, zName, tName, parameter,
erddapCoords$erddapXcoord, erddapCoords$erddapYcoord,
erddapCoords$erddapTcoord, erddapCoords$erddapZcoord,
verbose, cache_remove )
}
if (!is.list(extract)) {
text1 <- "There was an error in the url call, perhaps a time out."
text2 <- "See message on screen and URL called"
print(paste(text1, text2))
return("URL call error")
}
extract <- structure(extract,
class = c('list', 'rxtracto3D'),
base_url = urlbase,
datasetid = attributes(dataInfo)$datasetid
)
}