[R-sig-Geo] Raster crosstab output: -1 values?

Lyndon Estes lyndon.estes at gmail.com
Fri Jan 4 04:37:53 CET 2013


Hi Megan,

You might be able to do both the resampling and reprojecting quite
handily with gdal, using a system call out of R.

Here are some lines from one of my scripts (I haven't looked at this
in a while and haven't verified now, so caveat emptor) that might
help:

# lc is a landsat landcover I was working with
TARGET_SRS <- "SR-ORG:6842" # Your projection. See spatialreference.org
TARGET_EXTENT <- paste(extent(lc)@xmin, extent(lc)@ymin,
extent(lc)@xmax, extent(lc)@ymax)
# In target units!
TARGET_RES <- "30" # Whatever you want the output resolution to be
output_name <- "your_grid_name.tif"
gdcall <- "/Library/Frameworks/GDAL.framework/Programs/gdalwarp"  #
path to gdalwarp on my mac
gdwarp <- paste(gdcall, "-t_srs", TARGET_SRS,  "-te", TARGET_EXTENT,
"-tr", TARGET_RES, lc at file@name,
                           "-of GTiff -co 'COMPRESS=LZW'", output_name)
# In the above you will probably want to insert the -r option with
bilinear specified (e.g. "-r bilinear" worked
# into the paste command), otherwise you will have a nearest neighbor
resample done.
setwd("path/where/your/grids/live")
system(gdwarp, intern = F)

Hope it works, but if so it should be pretty quick.  Note, you can
also specify NA values:

http://www.gdal.org/gdalwarp.html

Best, Lyndon



More information about the R-sig-Geo mailing list