raster CRS changes on tis own
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
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