-
Notifications
You must be signed in to change notification settings - Fork 0
/
.Rhistory
512 lines (512 loc) · 20.7 KB
/
.Rhistory
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
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd)
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd)
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd)
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd)
plot(result_raster)
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd)
plot(result_raster)
#' @title Compute Solar Radiation from NetCDF Data using Angstrom-Prescott Model
#'
#' @description
#' The \code{radNC} function computes daily solar radiation for each pixel in a SpatRaster object based on sunshine duration data from a NetCDF file. It uses the Angstrom-Prescott model implemented in the \code{sol_rad} function.
#'
#' @param ssd SpatRaster object. A multi-layer raster containing sunshine duration data (in hours) for each pixel, with one layer per day.
#' @param show_progress Logical. Whether to display a progress bar during the computation. Defaults to \code{TRUE}.
#'
#' @return A SpatRaster object where each layer contains the computed solar radiation (in MJ/m²/day) for each corresponding day.
#'
#' @details
#' This function extracts the sunshine duration from each pixel in the \code{ssd} SpatRaster, and computes solar radiation for each pixel and each layer (day) using the Angstrom-Prescott model. The model requires sunshine duration, latitude, and the corresponding date as inputs. If any required attributes are missing, a warning or error is issued.
#'
#' @examples
#' \dontrun{
#' library(terra)
#' filename <- system.file("extdata", "ssd.nc", package = "sola")
#' ssd <- rast(filename)
#' result_raster <- radNC(ssd)
#' plot(result_raster)
#' }
#'
#' @export
radNC <- function(ssd, show_progress = TRUE) {
# Check if there is a time attribute
if (is.null(terra::time(ssd))) {
stop("The input SpatRaster object does not have a valid time attribute. Please ensure the input is a valid NetCDF file with time layers.")
}
# Extract the longitude and latitude of all pixels
lons <- terra::xFromCell(ssd, 1:ncell(ssd))
lats <- terra::yFromCell(ssd, 1:ncell(ssd))
# Check for valid longitude and latitude attributes
if (is.null(lons) || is.null(lats)) {
stop("The input SpatRaster object does not have valid longitude and latitude information.")
}
# Extract all time attributes
dates <- terra::time(ssd)
# Check if the time attribute is empty
if (any(is.na(dates))) {
warning("Some layers in the SpatRaster object do not have valid time information.")
}
# Define the calculation function
calc_solar_radiation <- function(ssd_layer, date) {
ssd_values <- terra::values(ssd_layer)
solar_radiation <- sol_rad(lat = lats, date = rep(date, length(lats)), ssd = ssd_values)
return(solar_radiation)
}
# Initialize a raster to store the calculation results
solar_radiation_raster <- terra::rast(ssd, nlyr = terra::nlyr(ssd))
# Start time measurement
start_time <- Sys.time()
# Use a standard for loop, with a progress bar
if (show_progress) {
pb <- txtProgressBar(min = 0, max = terra::nlyr(ssd), style = 3)
}
for (i in 1:terra::nlyr(ssd)) {
# Calculate solar radiation for each layer
solar_radiation_raster[[i]] <- terra::setValues(solar_radiation_raster[[i]], calc_solar_radiation(ssd[[i]], dates[i]))
# Update the progress bar
if (show_progress) {
setTxtProgressBar(pb, i)
# Estimate remaining time
elapsed_time <- Sys.time() - start_time
estimated_time <- (elapsed_time / i) * (terra::nlyr(ssd) - i)
cat(sprintf("\rEstimated remaining time: %s", format(estimated_time, digits = 2)))
}
}
# Close the progress bar
if (show_progress) {
close(pb)
}
# Calculate total elapsed time
total_time <- Sys.time() - start_time
# Display total time if progress is shown
if (show_progress) {
cat(sprintf("\nEstimated remaining time: %s\n", format(estimated_time, digits = 2)))
cat(sprintf("Total computation time: %s\n", format(total_time, digits = 2)))
}
# Return the raster with the calculation results
return(solar_radiation_raster)
}
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd)
plot(result_raster)
?declination
declination(as.Date("2024-03-15")) # Returns in radians (default)
declination(75) # Julian day in degrees
declination(as.Date("2024-03-15")) # Returns in radians (default)
declination(as.Date("2024-03-15")) # Returns in radians (default)
declination(as.Date("2024-03-15")) # Returns in radians (default)
declination(as.Date("2024.03.15"), TRUE) # Returns in degrees
declination(as.Date("2024-03-15")) # Returns in radians (default)
declination(as.Date("2024.03.15"), TRUE) # Returns in degrees
declination(as.Date("2024-03-15")) # Returns in radians (default)
declination("2024-03-15", TRUE) # Returns in degrees
declination(75) # Julian day in degrees
library(FAO56)
declination(as.Date("2024-03-15"))
SolDec(as.Date("2024-03-15"))
declination(as.Date("2024-03-15"))
SolDec(as.Date("2024-03-15"))
declination(as.Date("2024-05-15"))
SolDec(as.Date("2024-05-15"))
declination(as.Date("2024-05-15"))
declination(as.Date("2024-05-15"))
SolDec(as.Date("2024-05-15"))
SolDec(as.Date("2024-05-15"))
declination(as.Date("2024-05-15"))
declination("2024-05-03") # Returns in radians (default)
declination(as.Date("2024-05-03"), TRUE) # Returns in degrees
declination(75) # Julian day in degrees
# Handling a large dataset of dates (e.g., 100,000 dates)
large_dates <- as.Date("2024-01-01") + 0:99999
result_large <- declination(large_dates)
declination("2024-05-03") # Returns in radians (default)
declination(as.Date("2024-05-03"), TRUE) # Returns in degrees
declination(75) # Input Julian day
?calc_Ra
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = -15.0, date = "2024-06-21")
ExRad(d_r = 0.985, omega_s = 1.527, phi = -0.35, delta = 0.12)
?calc_Rs
# Estimate solar radiation for a single location and date
calc_Rs(lat = 35.0, date = as.Date("2023-03-15"), ssd = 8)
# Estimate solar radiation for a single location and date
calc_Rs(lat = 35.0, date = as.Date("2023-03-15"), ssd = 8)
?daylength
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-15, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-15, 45), date = c("2024-06-21", "2023-12-21"))
?calc_Ra
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = -15.0, date = "2024-06-21")
?calc_Rs
# Estimate solar radiation for a single location and date
calc_Rs(lat = 35.0, date = as.Date("2023-03-15"), ssd = 8)
# Estimate solar radiation for a single location and date
calc_Rs(lat = 35.0, date = as.Date("2023-03-15"), ssd = 8)
# Estimate solar radiation for a single location and date
calc_Rs(lat = 35.0, date = as.Date("2023-03-15"), ssd = 8)
# Example call to the function
calc_Rs(lat = 35.0, date = as.Date("2023-03-15"), ssd = 8)
SolDec("2020/08/25")
calc_Ra("2020/08/25")
?ap
calc_Rs(lat = 35.0, date = as.Date("2023-03-15"), Ra=12)
# 确保sirad包已经安装并加载
if (!require(sirad)) install.packages("sirad")
library(sirad)
# 定义共用的参数
lat <- 35.0
date <- as.Date("2023-03-15")
ssd <- 8
# 使用sirad包中的ap函数
ap_result <- ap(days = date, lat = lat, lon = NA, extraT = NULL, A = NA, B = NA, SSD = ssd)
ap(days = date, lat = lat, lon = NA, extraT = NULL, A = 0.25, B = 0.5, SSD = ssd)
calc_Rs(lat = lat, date = date, ssd = ssd, N = NULL, Ra = NULL, A = 0.25, B = 0.50)
calc_Rs(lat = lat, date = date, ssd = ssd, N = NULL, Ra = NULL, A = 0.25, B = 0.50)
?radNC
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd)
plot(result_raster)
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd)
plot(result_raster)
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd,na_neg = FLASE)
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd,na_neg = FALSE)
plot(result_raster)
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd,na_neg = FALSE)
# 保存结果为NetCDF文件
writeCDF(rad, "D:/CN05/roi/radiation_2020_2022.nc", overwrite = TRUE)
?sunset_ha
sunset_ha(as.Date("2023-03-15"), lat = 35) # Returns in radians (default)
sunset_ha("2023-03-15", lat = 35, in_degrees = TRUE) # Returns in degrees
sunset_ha(as.Date("2023-03-15"), lat = 90)
sunset_ha(as.Date("2023-03-15"), lat = 70)
sunset_ha(as.Date("2023-03-15"), lat = 80)
sunset_ha(as.Date("2023-03-15"), lat = 85)
?calc_Ra
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = -15.0, date = "2024-06-21")
# Example with a large dataset
set.seed(123)
test_lat <- runif(10000, -89, 89) # Generate random latitudes within recommended range
test_date <- as.Date("2023-01-01") + sample(0:364, 10000, replace = TRUE) # Generate random dates in 2023
result_ext_rad <- calc_Ra(test_lat, test_date)
print(head(result_ext_rad))
test_lat
calc_Ra(lat = 78.4031766, date = "2024-06-21")
calc_Ra(lat = 35.0, date = "2023-03-15")
calc_Ra(lat = 60, date = "2024-06-21")
calc_Ra(lat = 65, date = "2024-06-21")
calc_Ra(lat = 70, date = "2024-06-21")
calc_Ra(lat = 69, date = "2024-06-21")
calc_Ra(lat = 65, date = "2024-06-21")
sunset_ha(lat = 65, date = "2024-06-21")
sunset_ha(lat = 69, date = "2024-06-21")
sunset_ha(lat = 70, date = "2024-06-21")
library(solrad)
sunset(lat = 70, date = "2024-06-21")
Sunset(lat = 70, date = "2024-06-21")
Sunset(Lat = 70, doy("2024-06-21"))
1+NA
?daylength
daylength(lat = 35, date = as.Date("2023-03-15"))
?daylength
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-15, 45), date = c("2024-06-21", "2023-12-21"))
x<-daylength(lat = c(-15, 45), date = c("2024-06-21", "2023-12-21"))
x
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-60, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-65, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-70, 45), date = c("2024-06-21", "2023-12-21"))
?sunset_ha
sunset_ha(as.Date("2023-03-15"), lat = 35) # Returns in radians (default)
sunset_ha("2023-03-15", lat = 35, in_degrees = TRUE) # Returns in degrees
sunset_ha(as.Date("2023-03-15"), lat = 65)
sunset_ha(as.Date("2023-03-15"), lat = 70)
sunset_ha(as.Date("2023-06-15"), lat = 70)
sunset_ha(as.Date("2023-04-15"), lat = 70)
sunset_ha(as.Date("2023-05-15"), lat = 70)
sunset_ha(as.Date("2023-06-15"), lat = 70)
sunset_ha(as.Date("2023-06-05"), lat = 70)
sunset_ha(as.Date("2023-06-03"), lat = 70)
sunset_ha(as.Date("2023-06-03"), lat = 70)
sunset_ha(as.Date("2023-06-04"), lat = 70)
sunset_ha(as.Date("2023-06-05"), lat = 70)
sunset_ha(as.Date("2023-06-05"), lat = 70)
sunset_ha(as.Date("2023-06-04"), lat = 70)
sunset_ha(as.Date("2023-05-04"), lat = 70)
?calc_Ra
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = -15.0, date = "2024-06-21")
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = 65.0, date = "2024-06-21")
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = 70.0, date = "2024-06-21")
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = 80.0, date = "2024-06-21")
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = 90.0, date = "2024-06-21")
calc_Ra(lat = 35.0, date = as.Date("2023-03-15"))
calc_Ra(lat = 99.0, date = "2024-06-21")
?calc_Rs
# Estimate for another date and location
calc_Rs(lat = -15.0, date = "2024-06-21", ssd = 10)
# Estimate for another date and location
calc_Rs(lat = -70.0, date = "2024-06-21", ssd = 10)
# Estimate for another date and location
calc_Rs(lat = 65.0, date = "2024-06-21", ssd = 10)
# Estimate for another date and location
calc_Rs(lat = 70.0, date = "2024-06-21", ssd = 10)
# Estimate for another date and location
calc_Rs(lat = 71.0, date = "2024-06-21", ssd = 10)
# Estimate for another date and location
calc_Rs(lat = 72.0, date = "2024-06-21", ssd = 10)
# Estimate for another date and location
calc_Rs(lat = 73.0, date = "2024-06-21", ssd = 10)
?daylength
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-15, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-15, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(65, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(85, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(86, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(80, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(80, 45), date = c("2024-03-21", "2023-12-21"))
?calc_Rs
# Estimate solar radiation for a single location and date
calc_Rs(lat = 35.0, date = as.Date("2023-03-15"), ssd = 8)
# Estimate for another date and location
calc_Rs(lat = -15.0, date = "2024-06-21", ssd = 10)
# Estimate solar radiation for multiple locations and dates
calc_Rs(lat = c(35, 20), date = c(as.Date("2023-03-15"), as.Date("2022-02-15")), ssd = c(8, 5))
# Performance test with large datasets
set.seed(123)
test_lat <- runif(10000, -90, 90)
test_date <- as.Date("2023-01-01") + sample(1:365, 10000, replace = TRUE)
test_ssd <- runif(10000, 0, 24)
system.time({
calc_Rs(lat = test_lat, date = test_date, ssd = test_ssd)
})
?radNC
## Not run:
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd)
plot(result_raster)
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd,na_neg = FALSE)
?sunset_ha
sunset_ha(as.Date("2023-03-15"), lat = 35) # Returns in radians (default)
sunset_ha("2023-03-15", lat = 35, in_degrees = TRUE) # Returns in degrees
sunset_ha(as.Date("2023-03-15"), lat = 35) # Returns in radians (default)
sunset_ha("2023-03-15", lat = 35, in_degrees = TRUE) # Returns in degrees
?radNC
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd)
plot(result_raster)
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd)
?declination
declination(as.Date("2024-05-03")) # Returns in radians (default)
declination("2024-05-03", TRUE) # Returns in degrees
declination(75) # Input Julian day
# Handling a large dataset of dates (e.g., 100,000 dates)
large_dates <- as.Date("2024-01-01") + 0:99999
result_large <- declination(large_dates)
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd)
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd)
?daylength
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-15, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(75, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(-15, 0), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 35, date = as.Date("2023-03-15"))
daylength(lat = c(75, 0), date = c("2024-06-21", "2023-12-21"))
daylength(lat = c(75, 0), date = c("2024-06-21", "2023-12-21"))
# Check declination values for summer and winter solstices
declination_summer <- declination(as.Date("2024-06-21"), in_degrees = TRUE)
declination_winter <- declination(as.Date("2023-12-21"), in_degrees = TRUE)
declination_summer
declination_winter
# Check sunset hour angle values for 75 degrees latitude and summer and winter solstices
sunset_ha_summer <- sunset_ha(as.Date("2024-06-21"), lat = 75, in_degrees = TRUE)
sunset_ha_winter <- sunset_ha(as.Date("2023-12-21"), lat = 75, in_degrees = TRUE)
sunset_ha_summer
sunset_ha_winter
# Convert latitude from degrees to radians
lat_rad <- 75 * (pi / 180)
# Check for polar day or night
is_polar_day <- lat_rad + declination_summer * (pi / 180) > pi / 2
is_polar_night <- lat_rad + declination_winter * (pi / 180) < -pi / 2
is_polar_day
is_polar_night
# Convert latitude from degrees to radians
lat_rad <- 75 * (pi / 180)
# Check for polar day or night
is_polar_day <- lat_rad + declination_summer * (pi / 180) > pi / 2
is_polar_night <- lat_rad + declination_winter * (pi / 180) < -pi / 2
# Print the results
cat("Is it polar day on summer solstice at 75 degrees latitude?", is_polar_day, "\n")
cat("Is it polar night on winter solstice at 75 degrees latitude?", is_polar_night, "\n")
# Convert latitude from degrees to radians
lat_rad <- 75 * (pi / 180)
# Check for polar day or night
is_polar_day <- lat_rad + declination_summer * (pi / 180) > pi / 2
is_polar_night <- lat_rad + declination_winter * (pi / 180) < -pi / 2
# Print the results
cat("Is it polar day on summer solstice at 75 degrees latitude?", is_polar_day, "\n")
cat("Is it polar night on winter solstice at 75 degrees latitude?", is_polar_night, "\n")
# Calculate day length for winter solstice at 75 degrees latitude
daylength_winter <- daylength(lat = 75, date = as.Date("2023-12-21"))
# Print the result
daylength_winter
daylength(lat = c(75, 0), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 75, date = "2024-06-21")
daylength(lat = 0, date = "2023-12-21")
daylength(lat = c(75, 0), date = c("2024-06-21", "2023-12-21"))
daylength(lat = 75, date = "2024-06-21")
daylength(lat = 0, date = "2023-12-21")
?sunset_ha
sunset_ha(as.Date("2023-03-15"), lat = 35) # Returns in radians (default)
sunset_ha("2023-03-15", lat = 35, in_degrees = TRUE) # Returns in degrees
?sunset_ha
sunset_ha(as.Date("2023-03-15"), latitudes = c(35, 45)) # Returns in radians (default)
sunset_ha(as.Date("2023-03-15"), latitudes = c(35, 45)) # Returns in radians (default)
sunset_ha("2023-03-15", latitudes = c(35, 45), in_degrees = TRUE) # Returns in degrees
sunset_ha(as.Date("2023-03-15"), latitudes = c(35, 65)) # Returns in radians (default)
sunset_ha("2023-03-15", latitudes = c(35, 45), in_degrees = TRUE) # Returns in degrees
sunset_ha(as.Date("2023-03-15"), latitudes = c(35, 65)) # Returns in radians (default)
sunset_ha("2023-03-15", latitudes = c(35, 65), in_degrees = TRUE) # Returns in degrees
sunset_ha(as.Date("2023-03-15"), latitudes = c(35, 65)) # Returns in radians (default)
sunset_ha(c("2023-03-15","2023-04-20", latitudes = c(35, 65), in_degrees = TRUE) # Returns in degrees
sunset_ha(as.Date("2023-03-15"), latitudes = c(35, 65)) # Returns in radians (default)
sunset_ha(c("2023-03-15","2023-04-20"), latitudes = c(35, 65), in_degrees = TRUE) # Returns in degrees
sunset_ha(as.Date("2023-03-15"), latitudes = c(35, 65)) # Returns in radians (default)
sunset_ha(c("2023-03-15","2023-04-20"), latitudes = c(0, 65), in_degrees = TRUE) # Returns in degrees
daylength(lat = c(-15, 45), date = c("2024-06-21", "2023-12-21"))
daylength(lat = -15, date = "2024-06-21")
daylength(lat = 45, date = "2023-12-21")
?calc_Ra
# Example with a large dataset
set.seed(123)
test_lat <- runif(10000, -89, 89) # Generate random latitudes within recommended range
test_date <- as.Date("2023-01-01") + sample(0:364, 10000, replace = TRUE) # Generate random dates in 2023
result_ext_rad <- calc_Ra(test_lat, test_date)
print(head(result_ext_rad))
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd)
?calc_Rs
# Performance test with large datasets
set.seed(123)
test_lat <- runif(10000, -90, 90)
test_date <- as.Date("2023-01-01") + sample(1:365, 10000, replace = TRUE)
test_ssd <- runif(10000, 0, 24)
system.time({
calc_Rs(lat = test_lat, date = test_date, ssd = test_ssd)
})
?radNC
## Not run:
library(terra)
filename <- system.file("extdata", "ssd.nc", package = "sola")
ssd <- rast(filename)
result_raster <- radNC(ssd)
plot(result_raster)
library(sola)
library(terra)
library(pbapply)
# 读取NetCDF文件
ssd <- rast("D:/CN05/roi/ssd_roi.nc")
# 选择子集
ssd <- ssd[[time(ssd) >= "2020-01-01"]]
rad <- radNC(ssd,na_neg = FALSE)
writeCDF(rad, "D:/CN05/roi/radiation_2020_2022.nc", overwrite = TRUE)