0

I am working with some Landsat 8 images from the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center. I am processing the images in R, and after masking the clouds, I want to project from the following CRS "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" (as they are defined in the original data) to NAD 1983 UTM 12N.

I apply de following function (because I have 70 images):

test <- raster(xmn=220785, xmx=455115, ymn=4029285, ymx=4267815, res=30)
crs(test) <- "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

vals <- 1:ncell(test)
test <- setValues(test, vals)

pred.list <- list(test)
names(pred.list) <- c("test")

NAD83proj <- "+init=epsg:26912 +proj=utm +zone=12 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

projtoNAD1983 <- function(x) {
    l8_NAD1983 <- list()
    for (i in 1:length(x)) { 
        print(i)
        l8_NAD1983[[i]] <- projectRaster(from="+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0",crs=NAD83proj, filename= paste0(names(pred.list[i]), "_nad83.tif"), type="GTiff", overwrite=T)
    }
    return(l8_NAD1983)
    gc()
}

test_nad1983 <- projtoNAD1983(pred.list) 

But the output GTiff has a deprecated crs: "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

I do not know what I do wrong. My version of rgdal is 1,1-1, GDAL 1.11.2 (2015/02/10).

Should I just assign the projection epsg:26912 manually?

  • Landsat 8 has either 15 or 30 m resolution. The differences between WGS84 and NAD83 are much less than that so don't bother reprojecting it. Just assign/define EPSG:26912 if needed. – mkennedy Dec 11 '15 at 18:10
1

Please always simplify your example to allow focusing on the essence of your question (no need for a function, or a loop here). It had an error (the from argument is wrong). I also made test smaller.

test <- raster(xmn=220785, xmx=224595, ymn=4263885, ymx=4267815, res=30)
crs(test) <- "+proj=utm +zone=12 +datum=WGS84 +units=m"
values(test) <- 1:ncell(test)
NAD83proj <- "+init=epsg:26912 +proj=utm +zone=12 +datum=NAD83 +units=m"
x <- projectRaster(from=test, crs=NAD83proj)

So you get

crs(x)
# CRS arguments:
# +proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

Which, I think, is the same thing as what you want, expressed differently. Somehow it states the ellipsoid but not the datum. Perhaps because of using an epsg code. I do not know what that code stands for --- I know you can look that up -- but why make things opaque and not simply ask for what you want:

NAD83proj2 <-"+proj=utm +zone=12 +datum=NAD83 +units=m"
y <- projectRaster(from=test, crs=NAD83proj2)
crs(y)

# CRS arguments:
#  +proj=utm +zone=12 +datum=NAD83 +units=m +ellps=GRS80 +towgs84=0,0,0 

Whether it is a good idea to do all this (see mkennedy's comment) is another matter. Note tha tin both cases you have: +towgs84=0,0,0, i.e., no difference.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.