1

I plotted a RasterBrick, but I want to know the coordinates of each corner. I see that xmin = 3263575, xmax = 3300872, etc., but how are those coordinates?

class      : RasterStack 
dimensions : 161, 161, 25921, 12  (nrow, ncol, ncell, nlayers)
resolution : 231.6564, 231.6564  (x, y)
extent     : 3263575, 3300872, -63010.57, -25713.9  (xmin, xmax, ymin, ymax)
crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs 
names      : X2010.01.01, X2010.01.17, X2010.02.02, X2010.02.18, X2010.03.06, X2010.03.22, X2010.04.07, X2010.04.23, X2010.05.09, X2010.05.25, X2010.06.10, X2010.06.26 
min values :        -0.3,        -0.3,        -0.3,        -0.3,        -0.3,        -0.3,        -0.3,        -0.3,        -0.3,        -0.3,        -0.3,        -0.3 
max values :      0.9825,      0.9545,      0.9787,      0.9893,      0.8823,      0.9836,      0.9428,      0.9922,      0.9841,      0.9547,      0.9576,      0.9652 
1
  • 2
    Why do you need all 4 corner coordinates? The extent gives you the diagional corners. Jul 27 '20 at 23:43
2

The extent of that raster object is expressed in units according to the crs. So, to get "coordinates" (presumably you're looking for longitude and latitude) you need to create a spatial object with that extent information, and convert it to the desired crs. Based entirely off this answer, you can convert your extent of interest using the following:

library(sp)
library(raster)

ext <- extent(3263575, 3300872, -63010.57, -25713.9)
ext
# You could directly reference your own raster by uncommenting the next line
# ext <- extent(myRasterBrick)

ext <- as(ext, "SpatialPolygons")
sp::proj4string(ext) <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
# Similarly, use the proj info from your raster to assign crs to extent object
# sp::proj4string(ext) <- proj4string(myRasterBrick)

e.geo <- sp::spTransform(ext, CRS("+proj=longlat +datum=WGS84 +no_defs 
                             +ellps=WGS84 +towgs84=0,0,0"))
e.geo
3
  • 1
    But, then to get the corner coordinates of the SpatialPolygons object you need to access the node coordinates in the @Polygon slot which can be tricky. If you just use sp::coordinates it will return the centroid. Jul 27 '20 at 23:41
  • @JeffreyEvans Once I get the centroid (3282224, -44362.23), how do I convert that to coordinates? Jul 31 '20 at 20:47
  • 1
    @LauraMills What do you mean, those are the [x,y] coordinates? What is it you are expecting? Based on your post your coordinates would be in a sinusoidal orbital projection, Ideally you want spatial data in projected space so units are in distance and not degrees. If you are wanting latitude and longitude in decimal degrees then you have to reproject your data into this coordinate system using projectExtent, projectRaster, st_transform or spTransform depending on data. I think that in order to understand what is going on here you need to read up on projections. Jul 31 '20 at 21:07
2

Here is how you return the four corner coordinates of the raster extent. If you use sp::coordinates on the extent object then the polygon centroid coordinate is returned.

Here is an example of your raster, note that I do not have to assign a projection because everything will remain in the same projection space without it defined.

library(raster)
r <- raster(extent(3263575, 3300872, -63010.57, -25713.9),
            resolution=231.6564)
    r[] <- runif(ncell(r))

Now, create the extent polygon and then pull the node coordinates for the polygon.

e <- as(extent(r), "SpatialPolygons")  
  ( corner.coords <- e@polygons[[1]]@Polygons[[1]]@coords[1:4,] )

Now, the easy way. The sf::st_coordinates function returns polygon node coordinates and not centroids so we can do this in one fell swoop. Note that we are coercing the SpatialPolygons to an sf object in the same call. The index subsetting is removing two reference columns and the last coordinate pair that is the closing coordinate (same as the first) for the polygon.

sf::st_coordinates(as(as(extent(r), "SpatialPolygons"), "sf"))[1:4,][,1:2]  

Lets plot the results and create a points feature class while we are at it.

p=SpatialPoints(sf::st_coordinates(as(as(extent(r), "SpatialPolygons"), "sf"))[1:4,][,1:2] )

plot(p,pch=19,cex=2)
  plot(r,add=TRUE)
2
  • Does this approach generate corner coordinates for all transformations? If we made a "perfect square" in, say, an Albers projection, and then transformed the extent information of that square into lat/lon, would the corner points be true representatives of the corners in the Albers projection? i.e. would the xmin, xmax, ymin, and ymax be retained? Jul 28 '20 at 18:58
  • 1
    It would be the corner coordinates of whatever projection the raster is in. The corners are based on the raster extent which are defined as the lower left (origin), upper right coordinates and array size based on [x,y] cell size. If you want corner coordinates in a different projection space, and do not want to transform your raster, you can incorporate raster::projectExtent into the above code. Jul 28 '20 at 19:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.