[R-sig-Geo] raster CRS changes on tis own

[Roger Bivand]

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:
>>
>>> http://dl.dropbox.com/u/3180464/mosaico_envi.dat
>>> http://dl.dropbox.com/u/3180464/mosaico_envi.hdr
>>>
>>
>> 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