ACCUEIL / HOME
DOWNLOAD : Tutoriel R_SIG 2019. Script et jeu de données.
Useful refs
https://science.nature.nps.gov/im/datamgmt/statistics/r/advanced/spatial.cfm
Sub select an area from a raster
|
library(raster)
r1 <- raster(nrow=10, ncol=10)
r1[] <- runif(ncell(r1))
plot(r1)
bb <- drawExtent() # This function requires the user to select directly on the plot the area he wants to extract.
# Very handy when you want to take some data around some interesting points.
# This function takes a raster layer and reduced it by taking the required extent.
# The extent can come from the drawExtent call (just above) or directly from another raster or shapefile.
r = crop(r1, bb)
plot(r) # Very useful to go from a France climatic layer to the French Alps for instance
|
MASK
|
# Create a new RasterLayer where all cells that are NA in mask are set to NA, and that has the same values as x in the other cells
# Very useful to remove some values in the sea for instance, or simply to make sure that all the raster layers which might be used
# later for projections for instance or to build a stack of rasters have all the same extent and boundaries.
r <- raster(ncol=10, nrow=10)
m <- raster(ncol=10, nrow=10)
r[] <- runif(ncell(r)) * 10
m[] <- runif(ncell(r))
m[m <0.5] <- NA
mr <- mask(r, m)
|
Extract values of a unique raster or a raster stack on a set of XY points
Equivalent to the 'Sample' function in ArcGIS, 'intersect points' in Hawth Tools in ArcGIS.
Similarly to ArcGIS, these methods return the values of a Raster object for the cells in which a set of points fall.
Values can also be extracted through bilinear interpolation of values of the four cells nearest to each point.
It is also possible to supply a buffer to select values for cells within a certain distance around each point (see ?extract).
|
# If all data have the same extent, projection and bounding box (which is not the case with the data above), this is much faster to first create a stack of raster.
MyData <- stack(ras1, ras2, ras3)
# Coord: whatever dataframe containing at least 2 columns with XY coordinates. Let's assume here X and Y are the first two columns.
MyXY_env <- extract(MyData, Coord[, 1:2])
# Then all the data in the stack are sempled onto the XY coordinates. This is rather fast, especially using stack.
# If the data do not have the same resolution, or extent or bounding box, you cannot make a stack.
# Find below a very simple script using the eval(parse(text)) call.
# Just create a vector with the names of the raster you want.
raster_names <- c("A", "B", "C")
for(i in raster_names)
eval(parse(text=paste("Coord$", i, " <- extract(", i, ", Coord[,1:2])", sep="")))
# This simple function takes the raster loaded previously and sample it onto the Coord[,1:2] and store directly into the one column of Coord named by the name of the raster.
# Note that the reading of the rasters can be done directly into the loop as well.
|
Clip Raster BY polygon
|
# polygon: POLY
# raster: RAS
plot(RAS)
plot(POLY, add=T)
limits <- POLY@polygons[[1]]@Polygons[[1]]@coords # change index if necessary
points(limits) # verification
xyRaster <- xyFromCell(RAS, 1:ncell(RAS))
index <- inout(xyRaster, limits, bound=F)
RAS2 <- RAS
RAS2[!index]=NA
RAS2 <- crop(RAS2, extent(POLY))
plot(RAS2)
|
Extract polygon (SpatialPolygon) values from points (SpatialPoints)
|
pip <- overlay(POINTS, POLYGONS)
values <- unlist(lapply(pip, function(x){POLYGONS@data[x,"FIELD_VALUE"] }))
|
Spatial sample
see spsample(sp)
Projections
see GIS projections
Select polygons by attribute:
|
polygon2 <- polygon[polygon$NOM_ATTRIBUT==Attribute_value_to_select,]
|
Select polygons by location:
Select POLY1 in POLY2
1) Get vertices of polygons that determines the limits:
class: list (from shapefiles: read.shapefile)
limits <- POLY2$shp$shp[[1]]
class: SpatialPolygonsDataFrame
limits <- POLY2@polygons[[1]]@Polygons[[2]]@coords # change index if necessary
plot(limits) #verification
see also get.pts (gpclib) !!!
2) Clipping:
library(splancs)
inPoly <- inout(coordinates(POLY1),limits)
NEWPOLY<- POLY1[inPoly,]
ACCUEIL / HOME
Comments (0)
You don't have permission to comment on this page.