3

I am new to R-Gis.

Here is what I need to do, I need to make a raster inside the larger polygon (image above) each cell has to be 10x10 m, mask it to cut cells outside the the larger polygon and then assign a different value for each cell in the smaller polygons. Furthermore, I'd like to have a list (dataframe) with the coordinate of the center of each cell and the value associated to it.

enter image description here

And here is my code (not much):

library(tidyverse)
library(sp)
library(raster)
coords <-structure(list(id = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
                                        3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L), .Label = c("1.1", "2.1", "3.1", 
                                                                                    "3.2"), class = "factor"), lon = c(-47.720315, -47.7135780555556, 
                                                                                                                       -47.7122261111111, -47.7175261111111, -47.7187488888889, -47.715595, 
                                                                                                                       -47.7135130555556, -47.7166461111111, -47.7181161111111, -47.7169469444444, 
                                                                                                                       -47.71553, -47.716925, -47.7161630555556, -47.71552, -47.7139961111111, 
                                                                                                                       -47.7148438888889), lat = c(-22.5507461111111, -22.5456230555556, 
                                                                                                                                                   -22.549705, -22.5537280555556, -22.550825, -22.5479419444444, 
                                                                                                                                                   -22.5494280555556, -22.5521630555556, -22.5507561111111, -22.549725, 
                                                                                                                                                   -22.5507461111111, -22.5515380555556, -22.54924, -22.5486061111111, 
                                                                                                                                                   -22.5494769444444, -22.55026)), class = "data.frame", .Names = c("id", 
                                                                                                                                                                                                                    "lon", "lat"), row.names = c(NA, -16L))

coord.poly <- coords %>%
  nest(-id) %>%
  mutate(Poly = purrr::map(data, ~dplyr::select(., lon, lat) %>% Polygon()),
         polys = map2(Poly, id, ~Polygons(list(.x),.y))
  ) %>%
  {SpatialPolygons(.$polys)}

proj4string(coord.poly) <- CRS("+proj=utm +zone=23 +south +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") 

r <- raster(coord.poly, ncol = 100, nrow = 100) #don't know how to set cell size in meters :/
r[] <- 1:length(r)

plot(mask(r, coord.poly), axes = TRUE)

EDIT: forgot to include the libraries

  • 1
    Welcome to GIS SE. Please edi your post to include a single, focused question. – Aaron Mar 21 '18 at 17:41
2

First, you cannot arbitrary assign a desired projection, you need to actually transform the data into a different projection.

Here we assign the correct projection and then transform the data into a UTM projection. It will also be easier to work with a SpatialPolygonsDataFrame so, we will coerce into this class as well.

utm = "+proj=utm +zone=23 +south +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

proj4string(coord.poly) <-"+proj=longlat +ellps=WGS84 +datum=WGS84" 

polys <- spTransform(coord.poly, utm)

polys <- SpatialPolygonsDataFrame(polys, 
             data.frame(ID=1:length(polys)),
                        match.ID = FALSE)

Now, we can create an empty raster to use in rasterizing your polygons. The result will be a unique raster value for each polygon.

r <- raster(extent(polys), crs=proj4string(polys),
            resolution=10)
r <- rasterize(polys, r, field=polys@data[,"ID"])

# if you just want the encompassing polygon you could use a 
#   bracket index to subset it:
#   r <- rasterize(polys[1,], r, field=polys@data[,"ID"][1])

plot(r)
  plot(polys, add=TRUE)

It is not clear exactly what you are after in manipulating your resulting raster but, here is an example where we assign random values to the raster values corresponding to polygon 4

r[r == 4] <- runif(length(r[r == 4]))   

plot(r)
  plot(polys, add=TRUE)

To get at the coordinates of the raster cells you can just coerce to a SpatialPointsDataFrame. The coordinates slot corresponds to each row in the @data slot containing a data.frame of the attributes.

r.pts <- as(r, "SpatialPointsDataFrame")
  head( coordinates(r.pts) )

It is redundant but you can add the coordinates to the data.frame using

r.pts@data <- data.frame(r.pts@data, coordinates(r.pts))
  head(r.pts@data)
  • Thank you so much @Jeffrey. Since I'm not familiar with GIS terms I think I was not clear about change the raster value, Indeed I was looking for somehow set an "attribute", like an additional information attached to the cell. I fact it would be a random value for each one depending on the polygon it is in – Aníbal Deboni Neto Mar 21 '18 at 18:26
  • In my example that is exactly what is going on. If you want to assign a single random value to a polygon ID you could do something like: r[r == 4] <- rep(runif(1),length(r[r == 4])) I this answered your question please accept the answer so it can be closed. – Jeffrey Evans Mar 21 '18 at 18:44
  • I got that, but it set the same value to the whole polygon, not one for each cell. I tryed thi r[r == 4] <- rep(runif(1, min=45, max=60),length(r[r == 4])) but didn't work – Aníbal Deboni Neto Mar 21 '18 at 18:49
  • The example in the answer r[r == 4] <- runif(length(r[r == 4])) assigns unique random values to each cell corresponding the the given polygon ID (eg., 4). The runif function is for generating random numbers drawn from a Gaussian distribution. You can set desired min and max values as well. – Jeffrey Evans Mar 21 '18 at 18:54

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.