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