Geographic Information Systems Asked by Laura Mills on September 27, 2021
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
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
Correct answer by JepsonNomad on September 27, 2021
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)
Answered by Jeffrey Evans on September 27, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP