forked from r-spatial/stars
-
Notifications
You must be signed in to change notification settings - Fork 0
/
curvilinear.Rout.save
115 lines (110 loc) · 4.17 KB
/
curvilinear.Rout.save
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
R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> suppressPackageStartupMessages(library(stars))
>
> s5p = system.file("sentinel5p/S5P_NRTI_L2__NO2____20180717T120113_20180717T120613_03932_01_010002_20180717T125231.nc", package = "starsdata")
> if (s5p != "") {
+
+ lat_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/latitude")
+ lon_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/longitude")
+ nit_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/nitrogendioxide_summed_total_column")
+ lat = read_stars(lat_ds)
+ lon = read_stars(lon_ds)
+ nit = read_stars(nit_ds)
+ nit[[1]][nit[[1]] > 9e+36] = NA
+
+ ll = setNames(c(lon, lat), c("x", "y"))
+ nit.c0 = st_as_stars(nit, curvilinear = ll)
+
+ # more direct method:
+ nit.c = read_stars(s5p, sub = "//PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/nitrogendioxide_summed_total_column",
+ curvilinear = c("//PRODUCT/longitude", "//PRODUCT/latitude"), driver = NULL)
+ if (inherits(nit.c[[1]], "units"))
+ nit.c[[1]] = units::drop_units(nit.c[[1]])
+ nit.c[[1]][nit.c[[1]] > 9e+36] = NA
+ all.equal(nit.c0, nit.c)
+ st_crs(nit.c) = 4326
+ print(nit.c)
+
+ if (capabilities()["png"]) {
+ png("nit1.png", 800, 800)
+ plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = TRUE, pch = 16)
+ maps::map('world', add = TRUE, col = 'red')
+ dev.off()
+
+ png("nit2.png", 800, 800)
+ plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = FALSE, border = NA)
+ maps::map('world', add = TRUE, col = 'red')
+ dev.off()
+
+ nit.c = stars:::st_downsample(nit.c, 8)
+ print(nit.c)
+
+ png("nit3.png", 800, 800)
+ plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = TRUE, pch = 16)
+ maps::map('world', add = TRUE, col = 'red')
+ dev.off()
+
+ png("nit4.png", 800, 800)
+ plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = FALSE, border = NA)
+ maps::map('world', add = TRUE, col = 'red')
+ dev.off()
+ }
+
+ }
/PRODUCT/longitude,
/PRODUCT/latitude,
/PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/nitrogendioxide_summed_total_column,
stars object with 3 dimensions and 1 attribute
attribute(s):
nitrogendioxide_summed_total_column
Min. :0e+00
1st Qu.:1e-04
Median :1e-04
Mean :1e-04
3rd Qu.:1e-04
Max. :5e-04
NA's :330
dimension(s):
from to offset delta refsys point values
x 1 450 NA NA WGS 84 NA [450x278] -5.81066 [°],...,30.9468 [°]
y 1 278 NA NA WGS 84 NA [450x278] 28.3605 [°],...,51.4686 [°]
time 1 1 NA NA NA NA NULL
x/y
x [x]
y [y]
time
curvilinear grid
stars object with 3 dimensions and 1 attribute
attribute(s):
nitrogendioxide_summed_total_column
Min. :0.00005
1st Qu.:0.00008
Median :0.00008
Mean :0.00009
3rd Qu.:0.00009
Max. :0.00023
NA's :37
dimension(s):
from to offset delta refsys point values x/y
x 1 57 NA NA WGS 84 NA [57x35] -5.81066 [°],...,30.8108 [°] [x]
y 1 35 NA NA WGS 84 NA [57x35] 28.6622 [°],...,51.4679 [°] [y]
time 1 1 NA NA NA NA NULL
curvilinear grid
null device
1
There were 15 warnings (use warnings() to see them)
>
> proc.time()
user system elapsed
7.340 0.172 7.507