Sorry, this should have been: Alternatively, you could create three RasterLayer objects etopo1 <- raster("Etopo1.tif", band=1) etopo2 <- raster("Etopo1.tif", band=2) etopo3 <- raster("Etopo1.tif", band=3)
On Fri, Jul 23, 2010 at 9:53 AM, Robert J. Hijmans <r.hijm...@gmail.com> wrote: > If you want a multi-layer object you can create a RasterBrick in stead > of a RasterLayer. I.e. 'brick()' in stead of 'raster()' > > library(raster) > etopo1 <- brick("Etopo1.tif") # change here > e <- extent(-90, -32, -60, 15) > r <- crop(etopo1, e) > > > Alternatively, you could create three RasterLayer objects > etopo1 <- brick("Etopo1.tif", band=1) > etopo2 <- brick("Etopo1.tif", band=2) > etopo3 <- brick("Etopo1.tif", band=3) > > > Robert > > On Fri, Jul 23, 2010 at 7:54 AM, Michael Sumner <mdsum...@gmail.com> wrote: >> I know this is rather simplistic, and has bugs that I just don't have >> the impetus to track down at the moment, but it should work for your >> problem. I didn't realize that raster doesn't support (AFAICS) >> multiband rasters, such as 3-band images. >> >> Please use with caution - there are bugs! and very bad practices . . . >> >> readpartGDAL <- function(x, xlim = NULL, ylim = NULL, ...) { >> require(rgdal) >> info <- GDALinfo(x) >> offs <- info[c("ll.x", "ll.y")] >> scl <- info[c("res.x", "res.y")] >> dimn <- info[c("columns", "rows")] >> ## brain dead, but easy >> xs <- seq(offs[1], by = scl[1], length = dimn[1]) + scl[1]/2 >> ys <- seq(offs[2], by = scl[2], length = dimn[2]) + scl[2]/2 >> >> if (!is.null(xlim)) { >> if (!is.numeric(xlim)) stop("xlim must be numeric") >> if (!length(xlim) == 2) stop("xlim must be of length 2") >> if (!diff(xlim) > 0) stop("xlim[1] must be less than xlim[2]") >> xind <- which(xs >= xlim[1] & xs <= xlim[2]) >> } >> >> if (!is.null(ylim)) { >> if (!is.numeric(ylim)) stop("ylim must be numeric") >> if (!length(ylim) == 2) stop("ylim must be of length 2") >> if (!diff(ylim) > 0) stop("ylim[1] must be less than ylim[2]") >> yind <- which(ys >= ylim[1] & ys <= ylim[2]) >> } >> ## probably need a sign check for info["ysign"] >> ## reverse for y/x order in readGDAL >> rgdal.offset <- rev(c(min(xind), dimn[2] - max(yind))) >> rgdal.dim <- rev(c(length(xind), length(yind))) >> readGDAL(x, offset = rgdal.offset, region.dim = rgdal.dim, ...) >> } >> >> f <- "Etopo1.tif" >> >> d <- readpartGDAL(f, xlim = c(-90, -32), ylim = c(-60, 15)) >> >> ## assuming you have 3-bands describing RGB >> ## (see ?SGDF2PCT >> cols <- SGDF2PCT(d) >> d$idx <- cols$idx >> image(d, "idx", col = cols$ct) >> >> >> >> On Sat, Jul 24, 2010 at 12:09 AM, Rodrigo Aluizio <r.alui...@gmail.com> >> wrote: >>> Well Etopo1 have the colored image as well >>> (ftp://ftp.ngdc.noaa.gov/mgg/global/relief/ETOPO1/image/color_etopo1_ice_full.tif.gz). >>> I'm not able to import it using readGDAL (because of the vector size), but >>> normally when I do so with georeferenced images, the object created have, >>> in the data slot, the three band (RGB) color codes as columns. Using >>> raster{raster} and then as(x,'SpatialGridDataFrame') these data seems to >>> get truncated and unrecognizable to image{sp}. >>> >>> Regards. >>> >>> Rodrigo >>> >>> -----Mensagem original----- >>> De: Michael Sumner [mailto:mdsum...@gmail.com] >>> Enviada em: sexta-feira, 23 de julho de 2010 10:53 >>> Para: Rodrigo Aluizio >>> Cc: R Help >>> Assunto: Re: [R-sig-Geo] RES: Import Part of a GeoTiff using readGDAL{rgdal} >>> >>> I don't know what you mean by "band data structure" and "original colors". >>> I just have the generic (ice/cell registered) Etopo1 which is elevation >>> values as 16-bit integers. Perhaps you could describe the actual file you >>> have more, where you got it, and how you have used it before. >>> >>> Cheers, Mike. >>> >>> On Fri, Jul 23, 2010 at 8:43 PM, Rodrigo Aluizio <r.alui...@gmail.com> >>> wrote: >>>> Excellent suggestion, but I'm facing another little problem. When >>>> importing through raster{raster}, I lost the band data structure (but the >>>> data seems to still be there) which turns impossible to reproduce the >>>> original colors. Is there a way to bypass this? >>>> >>>> Thanks. >>>> >>>> Rodrigo. >>>> >>>> -----Mensagem original----- >>>> De: Michael Sumner [mailto:mdsum...@gmail.com] Enviada em: >>>> quinta-feira, 22 de julho de 2010 12:50 >>>> Para: Rodrigo Aluizio >>>> Cc: R Help >>>> Assunto: Re: [R-sig-Geo] Import Part of a GeoTiff using >>>> readGDAL{rgdal} >>>> >>>> With readGDAL you must figure out the impled coordinates, referenced from >>>> the top left - it's confusing, but with some experimenting you will get it. >>>> >>>> But, far easier is to use the raster package, which will let you deal with >>>> the entire Etopo1 if you want by accessing the file as needed, and there >>>> are simpler methods to crop the big raster - and coerce to >>>> SpatialGridDataFrame if need be. >>>> >>>> E.g. >>>> ## all this works pretty fast as raster takes care of the memory >>>> library(raster) >>>> etopo1 <- raster("Etopo1.tif") >>>> e <- extent(-90, -32, -60, 15) >>>> r <- crop(etopo1, e) >>>> plot(r) >>>> >>>> ## if we want to enforce the full memory use as an sp object x <- >>>> as(r, "SpatialGridDataFrame") >>>> image(x) >>>> summary(x) >>>> Object of class SpatialGridDataFrame >>>> Coordinates: >>>> min max >>>> s1 -90 -32 >>>> s2 -60 15 >>>> Is projected: FALSE >>>> proj4string : >>>> [+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0] Number >>>> of points: 2 Grid attributes: >>>> cellcentre.offset cellsize cells.dim >>>> s1 -89.99167 0.01666667 3480 >>>> s2 -59.99167 0.01666667 4500 Data attributes: >>>> Min. 1st Qu. Median Mean 3rd Qu. Max. >>>> -8086 -4188 -2583 -2025 148 6494 >>>> >>>> >>>> >>>> sessionInfo() >>>> R version 2.11.1 (2010-05-31) >>>> x86_64-pc-mingw32 >>>> >>>> locale: >>>> [1] LC_COLLATE=English_Australia.1252 LC_CTYPE=English_Australia.1252 >>>> [3] LC_MONETARY=English_Australia.1252 LC_NUMERIC=C [5] >>>> LC_TIME=English_Australia.1252 >>>> >>>> attached base packages: >>>> [1] stats graphics grDevices utils datasets methods base >>>> >>>> other attached packages: >>>> [1] rgdal_0.6-27 raster_1.2-6 sp_0.9-65 >>>> >>>> loaded via a namespace (and not attached): >>>> [1] grid_2.11.1 lattice_0.18-8 tools_2.11.1 >>>> >>>> On Fri, Jul 23, 2010 at 1:14 AM, Rodrigo Aluizio <r.alui...@gmail.com> >>>> wrote: >>>>> Hi List. >>>>> >>>>> I’m trying to import only a part of interest from a raster image map >>>>> (let say, ETOPO 1 Global Relief GeoTiff). I realize that, PROBABLY, >>>>> the readGDAL function can deal with my issue through the arguments: >>>>> >>>>> >>>>> >>>>> offset - Number of rows and columns from the origin (usually the >>>>> upper left >>>>> corner) to begin reading from; presently ordered (y,x) - this may >>>>> change >>>>> >>>>> region.dim - The number of rows and columns to read from the dataset; >>>>> presently ordered (y,x) - this may change >>>>> >>>>> >>>>> >>>>> But, I don’t know exactly how to specify that I want the rows (and >>>>> columns?) which show only South America, I do know the SA coordinate >>>>> limits, but not the corresponding rows. >>>>> >>>>> The real problem is that I don’t know how to use the function >>>>> resources properly. I’m trying to limit the region during import to >>>>> avoid memory limit issues, once importing the ETOPO 1 Geotiff tries to >>>>> generate a huge vector. >>>>> >>>>> Really Sorry About It, Thanks for Your Patience and Attention. >>>>> >>>>> >>>>> >>>>> ------------------------------------------------------------- >>>>> >>>>> MSc. <mailto:r.alui...@gmail.com> Rodrigo Aluizio >>>>> >>>>> Centro de Estudos do Mar/UFPR >>>>> Laboratório de Micropaleontologia >>>>> Avenida Beira Mar s/n - CEP 83255-000 Pontal do Paraná - PR - Brasil >>>>> >>>>> >>>>> >>>>> >>>>> [[alternative HTML version deleted]] >>>>> >>>>> >>>>> _______________________________________________ >>>>> R-sig-Geo mailing list >>>>> R-sig-Geo@stat.math.ethz.ch >>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>> >>>>> >>>> >>>> >>>> >>>> -- >>>> Michael Sumner >>>> Institute for Marine and Antarctic Studies, University of Tasmania Hobart, >>>> Australia >>>> e-mail: mdsum...@gmail.com >>>> >>>> _______________________________________________ >>>> R-sig-Geo mailing list >>>> R-sig-Geo@stat.math.ethz.ch >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>> >>> >>> >>> >>> -- >>> Michael Sumner >>> Institute for Marine and Antarctic Studies, University of Tasmania >>> Hobart, Australia >>> e-mail: mdsum...@gmail.com >>> >>> >> >> >> >> -- >> Michael Sumner >> Institute for Marine and Antarctic Studies, University of Tasmania >> Hobart, Australia >> e-mail: mdsum...@gmail.com >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@stat.math.ethz.ch >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo