R/dev/gemSST.R

Defines functions readGEM image.GEM

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)
}
mdsumner/mdsutils documentation built on May 22, 2019, 4:45 p.m.