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 <[email protected]> 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 <[email protected]> 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:[email protected]]
>> 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 <[email protected]> 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:[email protected]] 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 <[email protected]>
>>> 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:[email protected]> 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
>>>> [email protected]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>>
>>>
>>>
>>>
>>> --
>>> Michael Sumner
>>> Institute for Marine and Antarctic Studies, University of Tasmania Hobart,
>>> Australia
>>> e-mail: [email protected]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [email protected]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>>
>> --
>> Michael Sumner
>> Institute for Marine and Antarctic Studies, University of Tasmania
>> Hobart, Australia
>> e-mail: [email protected]
>>
>>
>
>
>
> --
> Michael Sumner
> Institute for Marine and Antarctic Studies, University of Tasmania
> Hobart, Australia
> e-mail: [email protected]
>
> _______________________________________________
> R-sig-Geo mailing list
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
