[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