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
