Making GIS and R play nicely - some tips?
Hi Andy, I've had pretty good luck reading large images line-by-line in rgdal. Its kind of slow, but not a bad as one might imagine. ArcGIS went down in flames trying to do any sort of processing on these files (10K x 10K). A better solution I think though will be to write space-time points into a postgis database (given lots of disk space) and use queries to subset the data as needed. I've found the database approach increadibly useful for working with BBS + climate data (~1e7 records). Ciao, THK
On Wed, 2005-03-02 at 11:26 -0500, abunn wrote:
Great. Thanks Roger and other others that replied.
But isn't the real challenge going to be shoehorning 250 by 1000 by 1000 by 8 into R at once - or are the data layers needed just one after the other? If I could grasp the workflow better, the advice above could become more focussed, I think.
I'm not positive on the workflow yet! That's one of the problems. I'm
finding it tough to wrap my head around the data since it is in both time
and space.
The data are not as huge as I thought at first. Each layer has 567 rows and
409 columns (231903 elements) but 137113 are no data. So, that makes the
actual amount of data to crunch much more manageable at 94790 elements per
layer.
The goal is to predict satellite reflectance as a function of climate. The
response data are monthly satellite reflectance over 19 years (So, that's
228 images). The predictor data are interpolated climate data and there
should be a minimum of three predictors and I'd like to use two more if they
are reliable. So, assuming three predictors the data would only run to 500mb
or so.
R > object.size(rep(runif(94790), times = 19*12*3)) / 2^20
[1] 494.6622
And if I can get all five then I'd be under a gig still.
R > object.size(rep(runif(94790), times = 19*12*5)) / 2^20
[1] 824.437
That's cumbersome but possible.
I've managed to read in a test example using erdas img files using rgdal
(thanks Tim Keitt):
R > library(rgdal)
R > junk <- list()
R > try1 <- GDAL.open("test.img")
R > getDriverLongName(getDriver(try1))
[1] "Erdas Imagine Images (.img)"
R > junk[[1]] <- c(getRasterData(try1))
R > GDAL.close(try1)
Closing GDAL dataset handle 00CACBF0... destroyed ... done.
R > try2 <- GDAL.open("test2.img")
R > getDriverLongName(getDriver(try2))
[1] "Erdas Imagine Images (.img)"
R > junk[[2]] <- c(getRasterData(try2))
R > GDAL.close(try2)
Closing GDAL dataset handle 00C7F180... destroyed ... done.
R > str(junk)
List of 2
$ : int [1:231903] 0 0 0 0 0 0 0 0 0 0 ...
$ : int [1:231903] 0 0 0 0 0 0 0 0 0 0 ...
So, to read all this in I can pipe a script to Arc to convert the grids to
images and loop through all the data. Does that sound like a reasonable
plan?
Thanks for all the help, Andy
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo