Skip to content

Import Part of a GeoTiff using readGDAL{rgdal}

8 messages · Rodrigo Aluizio, Michael Sumner, Robert J. Hijmans

Original
Michael Sumner
# 
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.aluizio at gmail.com> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rodrigo Aluizio
    

    

      
      #
    

  


  

      
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:mdsumner at 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.aluizio at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

        
--
Michael Sumner
Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia
e-mail: mdsumner at gmail.com

  
  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
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.aluizio at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rodrigo Aluizio
    

    

      
      #
    

  


  

      
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:mdsumner at 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.aluizio at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
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.aluizio at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Robert J. Hijmans
    

    

      
      #
    

  


  

      
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 <mdsumner at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Robert J. Hijmans
    

    

      
      #
    

  


  

      
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.hijmans at gmail.com> wrote: