raster CRS changes on tis own
Roger, I'm not dumping anything on anybody: those are links that you download if you want to, and that I posted only after you requested an example. Also, could you post your session info? According to what you say and what I reproduce we get different results for the same file and objects. Agus
On Tue, Nov 13, 2012 at 12:40 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
On Tue, 13 Nov 2012, Agustin Lobo wrote:
I do not think dumping 40MB on people is a helpful example. The idea is that
you rather go through the trouble of re-creating the problem from a small
data set, preferably built-in, because very oftem the process of re-creating
the problem solves it, that is, you find out which step you took that
deserved more thought.
It is also really helpful to have the input code verbatim, rather than
interspersed with output, as this is easier to copy & paste. In addition,
please use spaces after commas in included code, to avoid problems with
representation by different email clients; using T for TRUE is bad practice,
as I could have assigned T <- logical(NA) or "boo!"!
After downloading your data, and:
library(raster)
library(rgdal)
vinedomosaicMCA03 <- brick("mosaico_envi.dat")
projection(vinedomosaicMCA03)
says: "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs"
and:
ds <- GDAL.open("mosaico_envi.dat")
getProjectionRef(ds)
GDAL.close(ds)
says: "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs "
Why then assign a projection? After:
proj <- CRS("+init=epsg:3042")
projection(vinedomosaicMCA03) <- proj
projection(vinedomosaicMCA03)
we have: "+init=epsg:3042 +proj=utm +zone=30 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
which is the same projection. Following subsetting:
vinedomosaicMCA03v2 <- subset(vinedomosaicMCA03, c(3, 1, 2, 4, 5, 6))
projection(vinedomosaicMCA03v2)
we again have the same projection: "+proj=utm +zone=30 +ellps=GRS80 +units=m
+no_defs"; as you know, the dropped terms are simply duplicates, and do not
change the spatial reference at all.
writeRaster(vinedomosaicMCA03v2, file="mosaico_enviv2", overwrite=TRUE,
format="GTiff")
which takes a long time as this is a very large object.
ds <- GDAL.open("mosaico_enviv2.tif")
getProjectionRef(ds)
GDAL.close(ds)
shows: "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs "
which is exactly as read into R from the input file. Using the ENVI format:
writeRaster(vinedomosaicMCA03v2, file="mosaico_enviv2", overwrite=TRUE,
format="ENVI")
ds <- GDAL.open("mosaico_enviv2.envi")
getProjectionRef(ds)
GDAL.close(ds)
shows: "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs "
So what is your problem? Please do check everything carefully, use another
data set, and check that you understand what you are doing before suggesting
that others should devote their time to your muddle?
Roger