Okay, maybe all you want to do is work out where your raster projects to?
Here's a raster in epsg:3035 - you have something like this?
> re
class : RasterLayer
dimensions : 10, 10, 100 (nrow, ncol, ncell)
resolution : 7785.786, 11293.11 (x, y)
extent : 3577016, 3654874, 3714447, 3827378 (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10
+x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
data source : in memory
names : layer
values : 0.001494624, 0.9948669 (min, max)
Now, to see where that is in lat-long, simply do:
> extent(projectExtent(re,CRS("+init=epsg:4326")))
class : Extent
xmin : -2.290917
xmax : -0.7665775
ymin : 56.07447
ymax : 57.08626
(But as I keep banging on, your grid in epsg:3035 won't line up with
lat-long lines in that box)
The other thing you can do is use spTransform to work out where the
four corners of your raster transform to.
This corners function gets the four corner coordinates of a raster:
corners <- function(r){
e = extent(r)
x=c(xmin(e),xmin(e),xmax(e),xmax(e))
y=c(ymin(e),ymax(e),ymin(e),ymax(e))
SpatialPoints(data.frame(x=x,y=y),proj4string=CRS(projection(r)))
}
and returns a SpatialPoints object in the original coordinate system:
> corners(re)
SpatialPoints:
x y
[1,] 3577016 3714447
[2,] 3577016 3827378
[3,] 3654874 3714447
[4,] 3654874 3827378
Coordinate Reference System (CRS) arguments: +init=epsg:3035 +proj=laea
+lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m
+no_defs
Let's see where these transform to:
> spTransform(corners(re),CRS("+init=epsg:4326"))
SpatialPoints:
x y
[1,] -1.972778 55.97462
[2,] -2.290917 56.97278
[3,] -0.738697 56.08479
[4,] -1.025412 57.08626
Coordinate Reference System (CRS) arguments: +init=epsg:4326
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
Useful?
Barry
How to translate coordinates from cylindrical to a Lambert Azimuthal Equal Area Projection
2 messages · Barry Rowlingson, Ivailo
On Tue, Apr 23, 2013 at 7:41 PM, Barry Rowlingson
<b.rowlingson at lancaster.ac.uk> wrote:
Okay, maybe all you want to do is work out where your raster projects to?
... Huge thanks, Barry! I haven't thought of the problem in this way, but your approach is definitely more simple and practical than the one I vaguely imagined ;-)
(But as I keep banging on, your grid in epsg:3035 won't line up with lat-long lines in that box) The other thing you can do is use spTransform to work out where the four corners of your raster transform to.
Yes, having read further to find a way solving my problem I realized that both grids won't match exactly! Your example demonstrated that quite succinctly; now I'll be careful while selecting the corners of the possible extent as planned (in GE)...
This corners function gets the four corner coordinates of a raster:
corners <- function(r){
e = extent(r)
x=c(xmin(e),xmin(e),xmax(e),xmax(e))
y=c(ymin(e),ymax(e),ymin(e),ymax(e))
SpatialPoints(data.frame(x=x,y=y),proj4string=CRS(projection(r)))
}
and returns a SpatialPoints object in the original coordinate system:
> corners(re)
SpatialPoints:
x y
[1,] 3577016 3714447
[2,] 3577016 3827378
[3,] 3654874 3714447
[4,] 3654874 3827378
Coordinate Reference System (CRS) arguments: +init=epsg:3035 +proj=laea
+lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m
+no_defs
Let's see where these transform to:
> spTransform(corners(re),CRS("+init=epsg:4326"))
SpatialPoints:
x y
[1,] -1.972778 55.97462
[2,] -2.290917 56.97278
[3,] -0.738697 56.08479
[4,] -1.025412 57.08626
Coordinate Reference System (CRS) arguments: +init=epsg:4326
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
Useful?
Super useful, indeed -- thanks again for your great input, Barry! Cheers, Ivailo -- UBUNTU: a person is a person through other persons.