Hi Robert,

Thank you for the detailed information.

I downloaded the 10 min SRTM30 DEM but was unable to run the aggregate
command..  The error message was :
Error in .local(x, ...) : object 'startrow' not found

I tried to use readGDAL() but then ran into the same issue with NAs after
the call to overlay().  It seems it does not quite cover the entire globe..


On Mon, Sep 14, 2009 at 2:01 AM, Robert J. Hijmans <r.hijmans@gmail.com>wrote:

> I consider GTOPO30 to be obsolete because of SRTM30, at least for tropical
> mountain areas where GTOPO30 is pretty bad in some areas (although not
> perhaps not that much different when aggregated to a 1 degree
> resolution)...
>
> There are a few different SRTM30 products available. On the worldclim
> website, http://www.worldclim.org/current you can download them at
> different resolutions, including at 10 minutes; so that gets close to what
> you need.
>
> if you download the 10 min alt file, you could further aggregate it with
> the raster package. You can use the esri format or the generic format, but
> in case of the generic format you need to add this line to the .hdr file
> first:
> PIXELTYPE SIGNEDINT
>
> then:
>
> install.packages("raster", repos="http://R-Forge.R-project.org")
> # there is a problem with raster and R-forge right now, but this should
> work again soon
> *
> *
> library(raster)
> # create an object
> r <- raster('alt.bil')
> # aggregate to 1 degree
> ra <- aggregate(r, 6)
> # extract x and y values for the cells
> xy <- xyFromCell(ra, 1:ncell(ra))
> #  extract elevation values
> v  <-  values(ra)
>
> xyv <- cbind(xy, v)
>
>
>
> On Sun, Sep 13, 2009 at 3:13 PM, Ashton Shortridge <ashton@msu.edu> wrote:
>
>> Hello Abe,
>>
>> 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:
>> http://eros.usgs.gov/products/elevation/gtopo30/gtopo30.html
>>
>> 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,
>> etc. See:
>> http://dss.ucar.edu/catalogs/ranges/range750.html
>> 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.
>>
>> Good luck!
>>
>> Ashton
>>
>>
>> --
>> Ashton Shortridge
>> Associate Professor                     ashton@msu.edu
>> Dept of Geography                       http://www.msu.edu/~ashton<http://www.msu.edu/%7Eashton>
>> 235 Geography Building          ph (517) 432-3561
>> Michigan State University               fx (517) 432-1671
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo@stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

	[[alternative HTML version deleted]]

