1# Author: Robert J. Hijmans
2# Date : June 2008
3# Version 1.0
4# Licence GPL v3
5
6
7.gdFixGeoref <- function(mdata) {
8	gdversion <- getOption('rasterGDALVersion')
9	test <- gdversion < '1.8.0'
10	if (test) {
11		if (! is.null(mdata) ) {
12			for (i in 1:length(mdata)) {
13				if (mdata[i] == "AREA_OR_POINT=Area") {
14					return(FALSE)
15				} else if (mdata[i] == "AREA_OR_POINT=Point") {
16					return(TRUE)
17				}
18			}
19		}
20	}
21	return(FALSE)
22}
23
24
25
26.rasterFromGDAL <- function(filename, band, type, sub=0, RAT=TRUE, silent=TRUE, warn=TRUE, crs="", ...) {
27
28	.requireRgdal()
29
30	if (sub > 0) {
31		gdalinfo <- rgdal::GDALinfo(filename, silent=TRUE, returnRAT=FALSE, returnCategoryNames=FALSE)
32		sub <- round(sub)
33		subdsmdata <- attr(gdalinfo, 'subdsmdata')
34
35		i <- grep(paste("SUBDATASET_", sub, "_NAME", sep=''), subdsmdata)
36		if (length(i) > 0) {
37			x <- subdsmdata[i[1]]
38			filename <- unlist(strsplit(x, '='))[2]
39		} else {
40			stop(paste('subdataset "sub=', sub, '" not available', sep=''))
41		}
42	}
43
44	suppressWarnings(
45		gdalinfo <- try ( rgdal::GDALinfo(filename, silent=silent, returnRAT=RAT, returnCategoryNames=RAT) )
46	)
47
48	if ( inherits(gdalinfo, "try-error")) {
49		gdalinfo <- rgdal::GDALinfo(filename, silent=silent, returnRAT=FALSE, returnCategoryNames=FALSE)
50		warning('Could not read RAT or Category names')
51	}
52
53	nc <- as.integer(gdalinfo[["columns"]])
54	nr <- as.integer(gdalinfo[["rows"]])
55
56	xn <- gdalinfo[["ll.x"]]
57	xn <- round(xn, digits=9)
58
59	xx <- xn + gdalinfo[["res.x"]] * nc
60	xx <- round(xx, digits=9)
61
62	yn <- gdalinfo[["ll.y"]]
63	yn <- round(yn, digits=9)
64	yx <- yn + gdalinfo[["res.y"]] * nr
65	yx <- round(yx, digits=9)
66
67	nbands <- as.integer(gdalinfo[["bands"]])
68
69	if (isTRUE(attr(gdalinfo, "ysign") == 1)) {
70		warning("data seems flipped. Consider using: flip(x, direction='y')")
71	}
72
73	rotated <- FALSE
74	if (gdalinfo['oblique.x'] != 0 | gdalinfo['oblique.y'] != 0) {
75		rotated <- TRUE
76
77		## adapted from rgdal::getGeoTransFunc
78		if (warn) {
79			warning('\n\n This file has a rotation\n Support for such files is limited and results of data processing might be wrong.\n Proceed with caution & consider using the "rectify" function\n')
80		}
81		rotMat <- matrix(gdalinfo[c('res.x', 'oblique.x', 'oblique.y', 'res.y')], 2)
82		ysign <- attr(gdalinfo, 'ysign')
83		rotMat[4] <- rotMat[4] * ysign
84
85		invMat <- solve(rotMat)
86
87		offset <- c(xn, yx)
88		trans <- function(x, inv=FALSE) {
89			if (inv) {
90				x <- t(t(x) - c(offset[1], offset[2]))
91				x <- round( x %*% invMat  + 0.5 )
92				x[x < 1] <- NA
93				x[x[,1] > nc  | x[,2] > nr, ] <- NA
94			} else {
95				x <- (x - 0.5) %*% rotMat
96				x <- t(t(x) + c(offset[1], offset[2]))
97			}
98			return(x)
99		}
100
101		crd <- trans(cbind(c(0, 0, nc, nc), c(0, nr, 0, nr))+0.5)
102		rot <- methods::new(".Rotation")
103
104		gtr <- gdalinfo[c('ll.x', 'res.x', 'oblique.x', NA, 'oblique.y', 'res.y')]
105		gtr[4] <- yx
106		gtr[6] <- gtr[6] * ysign
107
108		rot@geotrans <- gtr
109		rot@transfun <- trans
110
111		xn  <- min(crd[,1])
112		xx  <- max(crd[,1])
113		yn  <- min(crd[,2])
114		yx  <- max(crd[,2])
115
116	}
117
118	mdata <- attr(gdalinfo, 'mdata')
119	fixGeoref <- FALSE
120	try( fixGeoref <- .gdFixGeoref(mdata), silent=TRUE )
121
122	# for ENVI files
123	bnames <- unique(mdata[grep("Band_", mdata)])
124	if (length(bnames) > 0) {
125		bn <- sapply(strsplit(bnames, '='), function(x) x[2])
126		bi <- gsub("Band_", "", sapply(strsplit(bnames, '='), function(x) x[1]))
127		bnames <- try(bn[order(as.integer(bi))], silent=TRUE)
128		if ( inherits(bnames, "try-error") ) {
129			bnames <- NULL
130		}
131	} else {
132		gobj <- rgdal::GDAL.open(filename)
133		bnames <- rep("", nbands)
134		for (i in 1:nbands) {
135			objbnd <- rgdal::getRasterBand(gobj, i)
136			bnames[i] <- rgdal::getDescription(objbnd)
137		}
138		rgdal::GDAL.close(gobj)
139	}
140
141	if (type == 'RasterBrick') {
142
143		r <- brick(ncols=nc, nrows=nr, xmn=xn, ymn=yn, xmx=xx, ymx=yx, crs="")
144		r@file@nbands <- r@data@nlayers <- nbands
145		band <- 1:nbands
146		#RAT <- FALSE
147
148	} else {
149
150		r <- raster(ncols=nc, nrows=nr, xmn=xn, ymn=yn, xmx=xx, ymx=yx, crs="")
151		r@file@nbands <- as.integer(nbands)
152		band <- as.integer(band)
153		if ( band > nbands(r) ) {
154			stop(paste("band too high. Should be between 1 and", nbands))
155			#if (warn) {
156				#stop("band too high. Set to nbands")
157			#}
158			#band <- nbands(r)
159		}
160		if ( band < 1) {
161			stop(paste("band should be 1 or higher"))
162			#if (warn) {
163				#stop("band too low. Set to 1")
164			#}
165			#band <- 1
166		}
167		r@data@band <- as.integer(band)
168		nbands <-1
169	}
170	if (rotated) {
171		r@rotated <- TRUE
172		r@rotation <- rot
173	}
174
175	prj <- attr(gdalinfo, 'projection')
176	if (!is.na(prj)) {
177		prjcom <- attr(prj, 'comment')
178		if ((!is.null(prjcom) && !is.na(prjcom))) {
179			prj <- prjcom
180		}
181	}
182	crs <- .getProj(prj, crs)
183	r@crs <- .CRS(crs, TRUE)
184	#r@crs <- .CRS(crs, FALSE)
185	# F to avoid warnings about other than WGS84 datums or ellipsoids
186
187#  	r@history[[1]] <- mdata
188
189
190	bi <- attr(gdalinfo, 'df')
191	GDType <- as.character(bi[['GDType']])
192	hasNoDataValues <- bi[['hasNoDataValue']]
193	NoDataValue <- bi[['NoDataValue']]
194
195#	if (getOption('rasterNewRGDALVersion')) {
196#		sbi <- attr(gdalinfo, 'sdf')
197#		Bmin <- sbi[['Bmin']]
198#		Bmax <- sbi[['Bmax']]
199#	} else {
200		Bmin <- bi[['Bmin']]
201		Bmax <- bi[['Bmax']]
202#	}
203
204
205	RATlist <- attr(gdalinfo, 'RATlist')
206	CATlist <- attr(gdalinfo, 'CATlist')
207
208	blockrows <- integer(nbands)
209	blockcols <- integer(nbands)
210
211	x <- rgdal::GDAL.open(filename, silent=TRUE)
212	ct <- rgdal::getColorTable( x )
213	if (! is.null(ct)) {
214		r@legend@colortable <- ct
215	}
216	for (i in 1:nbands) {
217		bs <- rgdal::getRasterBlockSize( rgdal::getRasterBand(x, i) )
218		blockrows[i] <- bs[1]
219		blockcols[i] <- bs[2]
220	}
221	rgdal::GDAL.close(x)
222
223	r@file@blockrows <- blockrows
224	r@file@blockcols <- blockcols
225
226
227	if (fixGeoref) {
228		message('Fixing "AREA_OR_POINT=Point" georeference')
229		rs <- res(r)
230		xmin(r) <- xmin(r) - 0.5 * rs[1]
231		xmax(r) <- xmax(r) - 0.5 * rs[1]
232		ymin(r) <- ymin(r) + 0.5 * rs[2]
233		ymax(r) <- ymax(r) + 0.5 * rs[2]
234	}
235
236	if (type == 'RasterBrick') {
237		ub <- unique(bnames)
238		if ((!all(ub == "")) && (length(ub) == nlayers(r))) {
239			names(r) <- bnames
240		} else {
241			names(r) <- rep(gsub(" ", "_", extension(basename(filename), "")), nbands)
242		}
243	} else {
244		lnames <- gsub(" ", "_", extension(basename(filename), ""))
245		if (nbands > 1) {
246			lnames <- paste0(lnames, '_', band)
247		}
248		names(r) <- lnames
249
250	}
251	r@file@name <- filename
252	r@file@driver <- 'gdal'
253
254
255	r@data@fromdisk <- TRUE
256
257	datatype <- "FLT4S"
258	minv <-	rep(Inf, nlayers(r))
259	maxv <-	rep(-Inf, nlayers(r))
260	try ( minv <- as.numeric( Bmin ) , silent=TRUE )
261	try ( maxv <- as.numeric( Bmax ) , silent=TRUE )
262	minv[minv == -4294967295] <- Inf
263	maxv[maxv == 4294967295] <- -Inf
264	try ( datatype <- .getRasterDType ( GDType[1] ), silent=TRUE )
265
266	if ( all(c(is.finite(minv), is.finite(maxv)))) {
267		r@data@haveminmax <- TRUE
268	}
269	r@file@datanotation <- datatype
270
271	r@data@min <- minv[band]
272	r@data@max <- maxv[band]
273
274	rats <- ! sapply(RATlist, is.null)
275	if (any(rats)) {
276		att <- vector(length=nlayers(r), mode='list')
277		for (i in 1:length(RATlist)) {
278			if (! is.null(RATlist[[i]])) {
279				dr <- data.frame(RATlist[[i]], stringsAsFactors=TRUE)
280				wv <- which(colnames(dr)=='VALUE')
281				if (length(wv) > 0) {
282					if (wv != 1) {
283						dr <- data.frame(dr[,wv,drop=FALSE], dr[,-wv,drop=FALSE])
284					}
285					colnames(dr)[1] <- 'ID'
286				} else {
287					if (all((colnames(dr) %in% c('Red', 'Green', 'Blue', 'Opacity', 'Histogram')))) {
288						# this is really a color table
289						rats[i] <- FALSE
290						if (is.null(ct)) {
291							r@legend@colortable <- grDevices::rgb(dr$Red, dr$Green, dr$Blue, dr$Opacity)
292						}
293						next
294					} else {
295						j <- which(colnames(dr) == 'Histogram')
296						if (isTRUE(j>0) & ncol(dr) > 1) {
297							dr <- data.frame(ID=0:(nrow(dr)-1), COUNT=dr[,j], dr[,-j,drop=FALSE])
298						} else {
299							dr <- data.frame(ID=0:(nrow(dr)-1), dr)
300						}
301					}
302				}
303				att[[i]] <- dr
304			}
305		}
306
307		r@data@attributes <- att[band]
308		r@data@isfactor <- rats[band]
309
310	} else {
311		cats <- ! sapply(CATlist, is.null)
312		if (any(cats)) {
313			att <- vector(length=nlayers(r), mode='list')
314			for (i in 1:length(CATlist)) {
315				if (! is.null(CATlist[[i]])) {
316					att[[i]] <- data.frame(ID=(1:length(CATlist[[i]]))-1, category=CATlist[[i]], stringsAsFactors=TRUE)
317				}
318			}
319			r@data@attributes <- att[band]
320			r@data@isfactor <- cats[band]
321		}
322	}
323	return(r)
324}
325
326
327