[R-sig-Geo] raster CRS changes on tis own
Roger Bivand
Roger.Bivand at nhh.no
Tue Nov 13 12:40:28 CET 2012
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
> to be read in with
> vinedomosaicMCA03 <- brick("mosaico_envi.dat")
>
> (both files must be present in the same directory)
>
> Agus
>
> On Tue, Nov 13, 2012 at 11:01 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>> Please provide a completely reproducible example, otherwise nobody can say
>> that this isn't simply a feature of your data. I can't see why helpers
>> should be expected to create an example for you, which may not correctly
>> represent your problem.
>>
>> Roger
>>
>>
>>
>> On Tue, 13 Nov 2012, Agustin Lobo wrote:
>>
>>> Sorry I forgot my session info:
>>>>
>>>> sessionInfo()
>>>
>>> R version 2.15.2 (2012-10-26)
>>> Platform: x86_64-pc-linux-gnu (64-bit)
>>>
>>> locale:
>>> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
>>> LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
>>> LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
>>> [7] LC_PAPER=C LC_NAME=C
>>> LC_ADDRESS=C LC_TELEPHONE=C
>>> LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>>>
>>> attached base packages:
>>> [1] stats graphics grDevices utils datasets methods base
>>>
>>> other attached packages:
>>> [1] rgdal_0.7-5 raster_2.0-12 sp_0.9-99
>>>
>>> loaded via a namespace (and not attached):
>>> [1] grid_2.15.2 lattice_0.19-30 tools_2.15.2
>>>
>>> Agus
>>>
>>> On Tue, Nov 13, 2012 at 10:41 AM, Agustin Lobo <alobolistas at gmail.com>
>>> wrote:
>>>>
>>>> Hi!
>>>> I modify the CRS of a raster object:
>>>>>
>>>>> vinedomosaicMCA03 <- brick(paste(mosdir,"mosaico_envi.dat",sep="/"))
>>>>> proj <- CRS("+init=epsg:3042")
>>>>> proj
>>>>
>>>> CRS arguments:
>>>> +init=epsg:3042 +proj=utm +zone=30 +ellps=GRS80
>>>> +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
>>>>>
>>>>> projection(vinedomosaicMCA03) <- proj
>>>>> vinedomosaicMCA03
>>>>
>>>> class : RasterBrick
>>>> dimensions : 1486, 1162, 1726732, 6 (nrow, ncol, ncell, nlayers)
>>>> resolution : 0.170881, 0.170881 (x, y)
>>>> extent : 578587.4, 578785.9, 4722025, 4722279 (xmin, xmax, ymin,
>>>> ymax)
>>>> coord. ref. : +init=epsg:3042 +proj=utm +zone=30 +ellps=GRS80
>>>> +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
>>>> data source :
>>>> /media/Iomega_HDD/UAVetal/vinedo/Clara_final/mosaico_envi.dat
>>>> names : mosaico_envi.1, mosaico_envi.2, mosaico_envi.3,
>>>> mosaico_envi.4, mosaico_envi.5, mosaico_envi.6
>>>>
>>>> But then, after a reordering, the CRS changes again on its own:
>>>>>
>>>>> vinedomosaicMCA03v2 <- subset(vinedomosaicMCA03,c(3,1,2,4,5,6))
>>>>> vinedomosaicMCA03v2
>>>>
>>>> class : RasterStack
>>>> dimensions : 1486, 1162, 1726732, 6 (nrow, ncol, ncell, nlayers)
>>>> resolution : 0.170881, 0.170881 (x, y)
>>>> extent : 578587.4, 578785.9, 4722025, 4722279 (xmin, xmax, ymin,
>>>> ymax)
>>>> coord. ref. : +proj=utm +zone=30 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0
>>>> +units=m +no_defs
>>>> names : mosaico_envi.3, mosaico_envi.1, mosaico_envi.2,
>>>> mosaico_envi.4, mosaico_envi.5, mosaico_envi.6
>>>> min values : NA, NA, NA, NA, NA, NA
>>>> max values : NA, NA, NA, NA, NA, NA
>>>>
>>>> The correct CRS is conserved if I write to GTiff:
>>>>>
>>>>> projection(vinedomosaicMCA03v2) <- proj
>>>>> vinedomosaicMCA03v2
>>>>
>>>> class : RasterStack
>>>> dimensions : 1486, 1162, 1726732, 6 (nrow, ncol, ncell, nlayers)
>>>> resolution : 0.170881, 0.170881 (x, y)
>>>> extent : 578587.4, 578785.9, 4722025, 4722279 (xmin, xmax, ymin,
>>>> ymax)
>>>> coord. ref. : +init=epsg:3042 +proj=utm +zone=30 +ellps=GRS80
>>>> +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
>>>> names : mosaico_envi.3, mosaico_envi.1, mosaico_envi.2,
>>>> mosaico_envi.4, mosaico_envi.5, mosaico_envi.6
>>>> min values : NA, NA, NA, NA, NA, NA
>>>> max values : NA, NA, NA, NA, NA, NA
>>>>
>>>>>
>>>>> writeRaster(vinedomosaicMCA03v2, file="mosaico_enviv2",
overwrite=T, format="GTiff")
>>>>
>>>> class : RasterBrick
>>>> dimensions : 1486, 1162, 1726732, 6 (nrow, ncol, ncell, nlayers)
>>>> resolution : 0.170881, 0.170881 (x, y)
>>>> extent : 578587.4, 578785.9, 4722025, 4722279 (xmin, xmax, ymin,
>>>> ymax)
>>>> coord. ref. : +proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs
>>>> data source :
>>>> /media/Iomega_HDD/UAVetal/vinedo/RNavarra/mosaico_enviv2.tif
>>>> names : mosaico_enviv2.1, mosaico_enviv2.2, mosaico_enviv2.3,
>>>> mosaico_enviv2.4, mosaico_enviv2.5, mosaico_enviv2.6
>>>> min values : 0, 0, 0, 0, 0, 0
>>>> max values : 175, 224, 179, 199, 196, 204
>>>>
>>>> But it is modified if I write to ENVI format:
>>>>
>>>>>
>>>>> writeRaster(vinedomosaicMCA03v2,file="mosaico_enviv2",overwrite=T,format="ENVI")
>>>>
>>>> class : RasterBrick
>>>> dimensions : 1486, 1162, 1726732, 6 (nrow, ncol, ncell, nlayers)
>>>> resolution : 0.170881, 0.170881 (x, y)
>>>> extent : 578587.4, 578785.9, 4722025, 4722279 (xmin, xmax, ymin,
>>>> ymax)
>>>> coord. ref. : +proj=utm +zone=30 +datum=NAD27 +units=m +no_defs
>>>> +ellps=clrk66 +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat
>>>> data source :
>>>> /media/Iomega_HDD/UAVetal/vinedo/RNavarra/mosaico_enviv2.envi
>>>> names : mosaico_enviv2.1, mosaico_enviv2.2, mosaico_enviv2.3,
>>>> mosaico_enviv2.4, mosaico_enviv2.5, mosaico_enviv2.6
>>>> min values : 0, 0, 0, 0, 0, 0
>>>> max values : 175, 224, 179, 199, 196, 204
>>>>
>>>>
>>>> Is this a bug in raster or rgdal, or am I doing something wrong?
>>>>
>>>> Agus
>>>
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> --
>> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
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
More information about the R-sig-Geo
mailing list