ashton at msu.edu
Mon Sep 14 00:13:20 CEST 2009
my thoughts to your questions are listed inline...
On Sunday 13 September 2009 15:05:04 Aben Woha wrote:
> Hello All,
> My questions are regarding the best way to obtain DEMs for the entire
> globe. I am trying to krige a dataset of global scope. The dataset has
> points at all latitude/longitude intersects and so there are 64,800 sample
> points. The problem that I am having is that the DEMs I have been using
> are too large (memory) to be used as R variables. I also have had trouble
> finding a way to download DEMs for the entire globe.
> 1. Given the size of my sample dataset: What is the maximum resolution of
> DEM that I should be using? Do I only need values at prediction points or
> will higher res yield higher accuracy?
As coarse as possible, since you only appear to need elevations at each 1x1
degree interval (note that many of these locations are nearby or identical for
extremely high and low latitudes, due to convergence of longitudes at the
poles). This question is more complicated if what you really want are average
elevations per 1-degree 'chunk', or you want to calculate aspect, slope, or
some other derivative.
I would check out GTOPO30, a 30-arc second global DEM:
The accuracy of this depends on the source data, which is reported at the site
listed above. There are other possibilities, but I don't know if the sources
are as recent or whether you'd get much in the way of accuracy statistics,
for more possibilities.
> 2. Is there a way to programmatically access DEMs using GDAL or R in
> general? ie. internally or through integration with McIDAS, NetCDF, etc.
Yes, check out the rgdal library. You can also use gdal (and maybe the python
bindings) to preprocess the dataset outside of R - perhaps to thin it to just
the lat-lons you need - prior to importing. I have used PyGDAL to access (and
sample from) 10s of gigabytes of tiled (e.g. in separate files) elevation data,
which I subsequently loaded into R. I think others here have done the same
from within GIS as well.
> 3. Do I need to break the DEMs and sample dataset into smaller tiles to
> avoid the memory issues? If so how do I do that in R / GDAL?
GTOPO30 comes in maybe 3 dozen tiles, ranging in size from 5-25 MB zipped. On
modern hardware it's not really that big, and R can probably load it all
without problem if you have enough RAM. My personal preference would probably
be to extract the elevations I wanted in home-brewed Python, write to some
sort of comma-delimited file, and read that into R - and I have no problems on
my 2 year old pc loading 64,000 records of xyz data, or actually many times
that number. There are surely simpler solutions that wouldn't involve much if
any coding, but I have direct experience with that approach.
Associate Professor ashton at msu.edu
Dept of Geography http://www.msu.edu/~ashton
235 Geography Building ph (517) 432-3561
Michigan State University fx (517) 432-1671
More information about the R-sig-Geo