0

I have the following grid of polygons (cf. at the end of the message). I would like to export it as a raster file that I could open with ArcMap.

I tried to actually add it to a grid that I already have (raster):

> grid <- readOGR("W:/FlorentDrafts/Grid24Cells","Grid1km_BalticSea_SAMPLE_24Cells")
OGR data source with driver: ESRI Shapefile 
Source: "W:/FlorentDrafts/Grid24Cells", layer: "Grid1km_BalticSea_SAMPLE_24Cells"
with 24 features
It has 1 fields
Integer64 fields read as strings:  Id 
> grid
class       : SpatialPolygonsDataFrame 
features    : 24 
extent      : 4530000, 4536000, 3517000, 3521000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
variables   : 1
names       : Id 
min values  : -1 
max values  :  0 

But I can not succeed using rasterize. I used the following lines:

coordinates(grid)
test_grid <- rasterFromXYZ(coordinates(grid))
rr <- rasterize(as(polys, "Spatial"), test_grid, getCover = TRUE, progress = "text").

It is possible to export such a grid of polygons as a ESRI raster file?

Here is the grid:

polys <-
structure(list(Id = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23), total = c(7, 3, 
8, 5, 2, 1, 2, 6, 4, 9, 8, 11, 2, 10, 1, 3, 4, 3, 2, 5, NA, 1, 
NA, NA), geometry = structure(list(structure(list(structure(c(4530000, 
4530000, 4531000, 4531000, 4530000, 3517000, 3518000, 3518000, 
3517000, 3517000), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
"sfg")), structure(list(structure(c(4531000, 4531000, 4532000, 
4532000, 4531000, 3517000, 3518000, 3518000, 3517000, 3517000
), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4532000, 4532000, 4533000, 4533000, 4532000, 
    3517000, 3518000, 3518000, 3517000, 3517000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4533000, 4533000, 4534000, 4534000, 4533000, 
    3517000, 3518000, 3518000, 3517000, 3517000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4534000, 4534000, 4535000, 4535000, 4534000, 
    3517000, 3518000, 3518000, 3517000, 3517000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4535000, 4535000, 4536000, 4536000, 4535000, 
    3517000, 3518000, 3518000, 3517000, 3517000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4530000, 4530000, 4531000, 4531000, 4530000, 
    3518000, 3519000, 3519000, 3518000, 3518000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4531000, 4531000, 4532000, 4532000, 4531000, 
    3518000, 3519000, 3519000, 3518000, 3518000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4532000, 4532000, 4533000, 4533000, 4532000, 
    3518000, 3519000, 3519000, 3518000, 3518000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4533000, 4533000, 4534000, 4534000, 4533000, 
    3518000, 3519000, 3519000, 3518000, 3518000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4534000, 4534000, 4535000, 4535000, 4534000, 
    3518000, 3519000, 3519000, 3518000, 3518000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4535000, 4535000, 4536000, 4536000, 4535000, 
    3518000, 3519000, 3519000, 3518000, 3518000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4530000, 4530000, 4531000, 4531000, 4530000, 
    3519000, 3520000, 3520000, 3519000, 3519000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4531000, 4531000, 4532000, 4532000, 4531000, 
    3519000, 3520000, 3520000, 3519000, 3519000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4532000, 4532000, 4533000, 4533000, 4532000, 
    3519000, 3520000, 3520000, 3519000, 3519000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4533000, 4533000, 4534000, 4534000, 4533000, 
    3519000, 3520000, 3520000, 3519000, 3519000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4534000, 4534000, 4535000, 4535000, 4534000, 
    3519000, 3520000, 3520000, 3519000, 3519000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4535000, 4535000, 4536000, 4536000, 4535000, 
    3519000, 3520000, 3520000, 3519000, 3519000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4530000, 4530000, 4531000, 4531000, 4530000, 
    3520000, 3521000, 3521000, 3520000, 3520000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4531000, 4531000, 4532000, 4532000, 4531000, 
    3520000, 3521000, 3521000, 3520000, 3520000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4532000, 4532000, 4533000, 4533000, 4532000, 
    3520000, 3521000, 3521000, 3520000, 3520000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4533000, 4533000, 4534000, 4534000, 4533000, 
    3520000, 3521000, 3521000, 3520000, 3520000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4534000, 4534000, 4535000, 4535000, 4534000, 
    3520000, 3521000, 3521000, 3520000, 3520000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
    structure(c(4535000, 4535000, 4536000, 4536000, 4535000, 
    3520000, 3521000, 3521000, 3520000, 3520000), .Dim = c(5L, 
    2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
    epsg = NA_integer_, proj4string = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"), .Names = c("epsg", 
"proj4string"), class = "crs"), class = c("sfc_POLYGON", "sfc"
), precision = 0, bbox = structure(c(4530000, 3517000, 4536000, 
3521000), .Names = c("xmin", "ymin", "xmax", "ymax"), class = "bbox"))), .Names = c("Id", 
"total", "geometry"), row.names = c(NA, -24L), class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_), class = "factor", .Label = c("constant", "aggregate", 
"identity"), .Names = c("Id", "total")))
0

Use the approach that I gave you in an earlier question and take account that rasterize doesn't work with sf objects:

polys <- as(polys,'Spatial')

r <- raster(nrows = 4, ncols = 6, xmn = 4530000, xmx = 4536000, ymn = 3517000, ymx = 3521000,
       crs = '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs')

rr <- rasterize(x = polys, y = r, field = 'total')

plot(rr)

enter image description here

  • rasterize does work with sf objects, even if not documented... – Robert Hijmans Dec 14 '17 at 3:47
1

Since you already have a regular raster, albeit represented as polygons, rasterize is a bit inefficient (not that it matters with the small example data set). Instead, you can do the following.

Create a raster that matches your data:

r <- raster(nrow=4, ncol=6, ext=extent(polys), crs="+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m")

Extract the centroids of the polygons

# with sf: 
xy <- do.call(rbind, st_geometry(st_centroid(polys)))
# with sp
# xy <- coordinates(polys)

Match coordinates to cells:

i <- cellFromXY(r, xy)

Assign cell values

r[i] <- polys$Id

# or 

r[i] <- polys$total

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.