Import Part of a GeoTiff using readGDAL{rgdal}
Locked
Import Part of a GeoTiff using readGDAL{rgdal}
 
|
Hi List.
Im 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 dont 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 dont know how to use the function resources properly. Im 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:[hidden email]> 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 [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Re: 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 <[hidden email]> 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:[hidden email]> 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 > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > -- Michael Sumner Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia e-mail: [hidden email] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
RES: Import Part of a GeoTiff using readGDAL{rgdal}
|
|
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:[hidden email]] 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 <[hidden email]> 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:[hidden email]> 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 > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > -- Michael Sumner Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia e-mail: [hidden email] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Re: 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 <[hidden email]> 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:[hidden email]] > 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 <[hidden email]> 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:[hidden email]> 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 >> [hidden email] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> >> > > > > -- > Michael Sumner > Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia > e-mail: [hidden email] > > _______________________________________________ > R-sig-Geo mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Michael Sumner Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia e-mail: [hidden email] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
RES: RES: Import Part of a GeoTiff using readGDAL{rgdal}
|
|
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:[hidden email]] 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 <[hidden email]> 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:[hidden email]] 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 <[hidden email]> 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:[hidden email]> 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 >> [hidden email] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> >> > > > > -- > Michael Sumner > Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia > e-mail: [hidden email] > > _______________________________________________ > R-sig-Geo mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Michael Sumner Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia e-mail: [hidden email] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Re: RES: Import Part of a GeoTiff using readGDAL{rgdal}
|
|
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 <[hidden email]> 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:[hidden email]] > 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 <[hidden email]> 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:[hidden email]] 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 <[hidden email]> 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:[hidden email]> 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 >>> [hidden email] >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >>> >> >> >> >> -- >> Michael Sumner >> Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia >> e-mail: [hidden email] >> >> _______________________________________________ >> R-sig-Geo mailing list >> [hidden email] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > > > -- > Michael Sumner > Institute for Marine and Antarctic Studies, University of Tasmania > Hobart, Australia > e-mail: [hidden email] > > -- Michael Sumner Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia e-mail: [hidden email] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Re: RES: Import Part of a GeoTiff using readGDAL{rgdal}
|
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 <[hidden email]> 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 <[hidden email]> 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:[hidden email]] >> 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 <[hidden email]> 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:[hidden email]] 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 <[hidden email]> 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:[hidden email]> 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 >>>> [hidden email] >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>> >>>> >>> >>> >>> >>> -- >>> Michael Sumner >>> Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia >>> e-mail: [hidden email] >>> >>> _______________________________________________ >>> R-sig-Geo mailing list >>> [hidden email] >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> >> >> >> -- >> Michael Sumner >> Institute for Marine and Antarctic Studies, University of Tasmania >> Hobart, Australia >> e-mail: [hidden email] >> >> > > > > -- > Michael Sumner > Institute for Marine and Antarctic Studies, University of Tasmania > Hobart, Australia > e-mail: [hidden email] > > _______________________________________________ > R-sig-Geo mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Re: RES: Import Part of a GeoTiff using readGDAL{rgdal}
|
|
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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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:[hidden email]] >>> 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 <[hidden email]> 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:[hidden email]] 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 <[hidden email]> 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:[hidden email]> 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 >>>>> [hidden email] >>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>> >>>>> >>>> >>>> >>>> >>>> -- >>>> Michael Sumner >>>> Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia >>>> e-mail: [hidden email] >>>> >>>> _______________________________________________ >>>> R-sig-Geo mailing list >>>> [hidden email] >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>> >>> >>> >>> >>> -- >>> Michael Sumner >>> Institute for Marine and Antarctic Studies, University of Tasmania >>> Hobart, Australia >>> e-mail: [hidden email] >>> >>> >> >> >> >> -- >> Michael Sumner >> Institute for Marine and Antarctic Studies, University of Tasmania >> Hobart, Australia >> e-mail: [hidden email] >> >> _______________________________________________ >> R-sig-Geo mailing list >> [hidden email] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
|
In reply to this post by Michael Sumner-2
The readpartGDAL function by Michael Sumner is simple and helpful for pruning larger raster/DEM maps.
Please bear in mind the function code is consistent with a simple and important assumption behind readGDAL that the data origin from which offset row and column counting begins to define new dim regions is "usually the upper left corner." Care is advised since some data sources (e.g. GOTOPO30) do not conform to the usual assumption. The function code needs a simple check to confirm the origin location and adjust as needed. When the origin is the lower left, for example, a rev() statement is needed before the Yaxis sequence named "ys". The section commented "brain dead and easy" will easily cause the function to die absent a control check. Many thanks to Michael for his original post. |
«
Return to R-sig-geo
|
1 view|%1 views
Loading...
| Powered by Nabble | Edit this page |