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

[Agustin Lobo]

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:
>
>> 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
>