[R-sig-Geo] Making GIS and R play nicely - some tips?

[Timothy H. Keitt]

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