raster CRS changes on tis own
On Tue, 13 Nov 2012, Agustin Lobo wrote:
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.
That is precisely why I asked you to provide another example of the problem. If the problem cannot be reproduced with other data on your system, then it is limited to your initial data set, or the specific code used when the unexpected behaviour occurred. If you can reproduce it for similar data sets, then it is your system. I did note that you are running rgdal 0.7-5, mine is 0.7-22 (the current CRAN version). In my case the crucial startup messages are: Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08 Path to GDAL shared files: /usr/local/share/gdal Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480] Path to PROJ.4 shared files: (autodetected) but yours are unknown. I install GDAL and PROJ.4 from source; like you, I am on: Platform: x86_64-unknown-linux-gnu (64-bit) (you are: x86_64-pc-linux-gnu (64-bit)) My intuition is that the problem is a stale GDAL on your system, but from reading the GDAL tickets and NEWS logs, I can't see an update to the ENVI driver that might explain where the odd extra tags are coming from, so it might be something else, assuming it is your system and not an idiosyncracy of your handling that data set. Roger
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
Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no