-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.Rmd
142 lines (104 loc) · 4.72 KB
/
README.Rmd
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
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%"
)
```
# L3bin
<!-- badges: start -->
[](https://github.com/hypertidy/L3bin/actions/workflows/R-CMD-check.yaml)
<!-- badges: end -->
The goal of L3bin is to provide the [Integerized Sinusoidal Binning Scheme](https://oceancolor.gsfc.nasa.gov/resources/docs/format/l3bins/) used by NASA's Ocean Biology processing Group for Level 3 data.
## Installation
You can install the development version of L3bin from [GitHub](https://github.com/) with:
``` r
# install.packages("devtools")
devtools::install_github("hypertidy/L3bin")
```
## Example
This is a basic example which shows how to create a very low resolution L3 bin scheme, plot it and perform some queries.
```{r example}
library(L3bin)
## basic example code
nr <- 6
bins <- L3bin(nr)
ex <- extent_from_bin(1:bins$totbins, nr)
plot(range(ex[,1:2]), range(ex[,3:4]), type = "n", asp = 1)
rect(ex[,1], ex[,3], ex[,2], ex[, 4])
points(lonlat_from_bin(1:bins$totbins, NUMROWS = nr))
```
Now we increas the resolution and use the available froms to crop, obtain bin boundaries, and provide the centre point of the bin from its index.
```{r twelve}
nr <- 12
bins <- L3bin(nr)
lbins <- crop_bins(bins, extent = c(-20, 30, -60, -10))
ex <- extent_from_bin(lbins, nr)
plot(range(ex[,1:2]), range(ex[,3:4]), type = "n", asp = 1)
rect(ex[,1], ex[,3], ex[,2], ex[, 4])
points(lonlat_from_bin(lbins, NUMROWS = nr))
ex <- extent_from_bin(1:bins$totbins, nr)
plot(range(ex[,1:2]), range(ex[,3:4]), type = "n", asp = 1)
rect(ex[,1], ex[,3], ex[,2], ex[, 4])
points(lonlat_from_bin(1:bins$totbins, NUMROWS = nr))
```
It should be clear that while we can easily up the detail and obtain very large numbers of bins, it's much smarter to keep things lazy and not materialize the actual bin numbers, boundary extents, or centre points. It very quickly can
become unwieldly. For example, the base resolution of SeaWiFS L3 bins in 1998 was a NUMROWS of 2160, which means a lot of unnecessary storage. Much better to avoid that when the functions provide the capability we need without generating the grid explicitly.
```{r rows2160}
(nbins_2160 <- L3bin(2160)$totbins)
## this would be nearly 25Mb just to list the bin index
nbins_2160 * 4 / 1e6
## 190Mb to store the bin boundaries
nbins_2160 * 4 * 8/ 1e6
## another 50Mb just for the centre coordinates
nbins_2160 * 8 / 1e6
```
A grid with 2160 rows represents approximately 9.28km^2 in bin area, and we may have much higher resolution that that for modern tasks. L3bin provides a simple core for this binning scheme.
### Warning
The bin number is currently stored as an integer, so beyond ~41000 ) rows this will overflow and we'd need to use a larger type to store it. Let us know if you need this, or encounter problems. 41000 rows is near a 500m bin size, which is the resolution of GEBCO 2021 for example.
Get lonlat from bin and vice versa.
```{r bin}
plot(range(ex[,1:2]), range(ex[,3:4]))
rect(ex[,1], ex[,3], ex[,2], ex[, 4])
text(lonlat_from_bin(1:bins$totbins, nr), lab = 1:bins$totbins)
abline(v = 147, h = -42, col = "firebrick")
bin_from_lonlat(147, -42, nr)
```
Using the crop from above, reproject after densifying (safely) with sf
```{r reproj}
nr <- 96
bins <- L3bin(nr)
ex <- extent_from_bin(1:bins$totbins, nr)
lbins <- crop_bins(bins, extent = c(100, 180, -90, -20))
ex <- extent_from_bin(lbins, nr)
## generally sf a bad idea because 4 numbers xmin,xmax,ymin,ymax become 10
## and we could use mesh3d or similar from silicate (SC0), but here's an example
## just to show a projected set of the bins
to_sf <- function(ex, crs = "OGC:CRS84") {
idx <- c(1, 2, 2, 1, 1,
3, 3, 4, 4, 3)
proto <- sf::st_polygon(list(cbind(0, c(1:4, 1))))
sf::st_sfc(lapply(split(t(ex), rep(1:nrow(ex), each = 4L)),
function(.x) {
proto[[1L]][] <- .x[idx]
proto
}), crs = crs)
}
sf::sf_use_s2(FALSE)
segit <- function(x, dx = 111111) {
crs <- sf::st_crs(x)
sf::st_set_crs(sf::st_segmentize(sf::st_set_crs(x, NA), dx), crs)
}
prj <- "+proj=laea +lat_0=-80"
## dx is 1 degree in straightforward sense
plot(sf::st_transform(segit(to_sf(ex), 11111), prj))
mm <- reproj::reproj(do.call(cbind, maps::map(plot = FALSE)[1:2]), prj, source = "+proj=longlat")[,1:2]
lines(mm)
```
## Code of Conduct
Please note that the L3bin project is released with a [Contributor Code of Conduct](https://contributor-covenant.org/version/2/1/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms.