GEM
function to read xlim/ylim/tlim based on lat/lon/time
- needs to convert from Mercator grid to rectilinear lat/lon
- need to flip over dateline as necessary
mkLookup function
- modified to detect irregularity and optimize appropriately (does it matter?)
- modified to optionally return the values, not a test (because the thing may not be a bit array)
by.segment is required - TRUE/FALSE
bit.array = TRUE/FALSE
Then we can
- read big set of SST (assume not too big)
- push through masking into bit array
- or make big array of data, with attributes for x, y, z
readGEM <- function(xlim = NULL, ylim = NULL, tlim = NULL,
longlat = TRUE,
filepath = "E:/backup/WD6400/GEM/GEMdata/temperature/temp_depth25.vrt")
{
## TODO
## generate Mercator versions
## detect non longlat case and subst based on that
gemx <- seq(-179.8332, 179.8332, length = 1080)
gemy <- c(-70.0568463846126, -69.942996316589, -69.8285222315121, -69.71342115295,
-69.5976900979515, -69.4813260771779, -69.3643260950375, -69.2466871498237,
-69.128406233856, -69.0094803336247, -68.8899064299388, -68.7696814980772,
-68.6488025079433, -68.5272664242239, -68.4050702065501, -68.2822108096636,
-68.1586851835854, -68.0344902737883, -67.9096230213743, -67.7840803632539,
-67.6578592323311, -67.5309565576911, -67.4033692647926, -67.2750942756636,
-67.1461285091017, -67.0164688808783, -66.886112303947, -66.755055688656,
-66.6232959429652, -66.4908299726674, -66.3576546816134, -66.2237669719426,
-66.0891637443168, -65.9538418981595, -65.8177983318995, -65.6810299432187,
-65.5435336293052, -65.4053062871107, -65.2663448136131, -65.1266461060833,
-64.9862070623577, -64.8450245811149, -64.7030955621578, -64.560416906701,
-64.4169855176622, -64.2727982999604, -64.127852160818, -63.9821440100684,
-63.8356707604697, -63.6884293280223, -63.5404166322937, -63.3916295967472,
-63.2420651490769, -63.0917202215482, -62.9405917513438, -62.7886766809152,
-62.6359719583404, -62.4824745376867, -62.3281813793799, -62.1730894505791,
-62.0171957255573, -61.860497186088, -61.7029908218378, -61.5446736307656,
-61.3855426195264, -61.2255948038831, -61.0648272091229, -60.9032368704808,
-60.7408208335688, -60.5775761548122, -60.413499901891, -60.2485891541888,
-60.0828410032473, -59.9162525532276, -59.7488209213781, -59.5805432385084,
-59.4114166494702, -59.2414383136441, -59.0706054054341, -58.8989151147674,
-58.7263646476019, -58.5529512264392, -58.3786720908458, -58.2035244979795,
-58.0275057231234, -57.8506130602261, -57.6728438224491, -57.4941953427208,
-57.3146649742968, -57.1342500913275, -56.9529480894323, -56.7707563862805,
-56.5876724221792, -56.4036936606673, -56.218817589117, -56.0330417193416,
-55.8463635882107, -55.6587807582706, -55.4702908183734, -55.2808913843112,
-55.0905800994578, -54.8993546354164, -54.707212692675, -54.5141520012669,
-54.320170321439, -54.1252654443258, -53.9294351926303, -53.7326774213113,
-53.5349900182766, -53.336370905083, -53.1368180376425, -52.9363294069343,
-52.7349030397235, -52.5325369992853, -52.3292293861357, -52.1249783387674,
-51.9197820343925, -51.71363868969, -51.5065465615596, -51.2985039478803,
-51.0895091882757, -50.8795606648832, -50.6686568031295, -50.4567960725105,
-50.2439769873767, -50.0301981077236, -49.815458039986, -49.5997554378375,
-49.383089002995, -49.1654574860262, -48.9468596871623, -48.7272944571145,
-48.5067606978942, -48.2852573636369, -48.0627834614295, -47.8393380521414,
-47.614920251258, -47.389529229718, -47.1631642147703, -46.9358244907462,
-46.7075094000082, -46.4782183437277, -46.2479507827503, -46.0167062384466,
-45.7844842935631, -45.5512845930755, -45.3171068450421, -45.0819508214588,
-44.8458163591133, -44.6087033604407, -44.3706117943781, -44.1315416972195,
-43.8914931734692, -43.6504663966952, -43.4084616103809, -43.1654791287748,
-42.921519337739, -42.6765826955953, -42.4306697339686, -42.1837810586277,
-41.9359173503226, -41.6870793656187, -41.4372679377270, -41.1864839773296,
-40.9347284734013, -40.6820024940259, -40.4283071872069, -40.1736437816736,
-39.9180135876793, -39.6614179977949, -39.4038584876942, -39.1453366169329,
-38.8858540297192, -38.6254124556771, -38.3640137106002, -38.1016596971973,
-37.8383524058289, -37.5740939152328, -37.3088863932405, -37.042732097483,
-36.7756333760845, -36.5075926683457, -36.2386125054145, -35.9686955109443,
-35.6978444017396, -35.4260619883883, -35.1533511758802, -34.8797149642104
)
gemt <- seq(ISOdatetime(1993, 10, 12, 0, 0, 0), length = 728, by = "7 days")
if (is.null(xlim)) xlim <- range(gemx)
if (is.null(ylim)) ylim <- range(gemy)
if (is.null(tlim)) tlim <- range(gemt[1:3])
if (!(length(xlim) == 2 & length(ylim) == 2 & length(tlim) == 2)) {stop("xlim, ylim, tlim must be of length 2, or unspecified")}
indX <- which(gemx >= xlim[1] & gemx <= xlim[2])
indY <- which(gemy >= ylim[1] & gemy <= ylim[2])
indT <- which(gemt >= tlim[1] & gemt <= tlim[2])
if (length(indX) < 1) {stop("x range out of bounds")}
if (length(indY) < 1) {stop("y range out of bounds")}
if (length(indT) < 1) {stop("t range out of bounds")}
if (!require(rgdal) ) {stop("rgdal package not available")}
print(range(indX))
print(range(indY))
print(range(indT))
print(c(length(gemy) - max(indY), min(indX) - 1))
print(c(length(indY), length(indX)))
d <- readGDAL(filepath, offset = c(length(gemy) - max(indY), min(indX) - 1), region.dim = c(length(indY), length(indX)), band = indT)
if (length(indT) > 1) {
a <- array(dim = c(length(indX), length(indY), length(indT)))
for (idim in 1:length(indT))
{
a[,,idim] <- as.image.SpatialGridDataFrame(d[idim])$z
}
} else
{
a <- as.image.SpatialGridDataFrame(d[1])$z
}
rm(d)
a <- list(x = gemx[indX], y = gemy[indY], z = a)
attr(a, "times") <- gemt[indT]
class(a) <- c(class(a), "GEM")
a
}
image.GEM <- function(x, bands = NULL, ...) {
#if (is.null(bands)) {print(paste("GEM contains ", dim(x$z), " bands, plotting the first one"))}
brks <- quantile(x$z, prob = seq(0, 1, by = 0.01), na.rm = TRUE)
cols <- heat.colors(length(brks) - 1)
for (i in 1:dim(x$z)[3]) {
image(x = x$x, y = x$y, z = x$z[,,i], col = cols, breaks = brks, ...)
}
invisible(x)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.