I am trying to rasterize a polygon in R using teh rasterize function but the output is not correct. I am comparing the result with the output produced by the rasterize tool available in QGIS

May data are:

> Fields2
class       : SpatialPolygonsDataFrame 
features    : 1293 
extent      : 219807.1, 230535.3, 4097488, 4104080  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 4
names       :  ID,       area, id_2,      mean 
min values  :   1, 100044.627,    1, 1001.9653 
max values  : 999,   9990.943,    1,  998.5807 

enter image description here

> test
class       : RasterLayer 
dimensions  : 660, 1074, 708840  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 219800, 230540, 4097480, 4104080  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       : fortesting 
values      : 126.7778, 1965.884  (min, max)

I am using the following command poly_to_raster<-rasterize(x=Fields2, y=test, field="mean")

The results looks like this: enter image description here

The white pixels do not have any value and are NA. Confirmed by looking at the pixel info of this raster in QGIS

enter image description here

I thouhgh there is someting wrong wity the shapefile I am using as reference in rasterize, but when I do the same process in QGIS, I get the right result: enter image description here

enter image description here

--- EDIT ---

The mean attribute comes from an extract function

meancalc<-extract(S_images_ts[[1]], Fields, fun=mean)
Fields@data<-data.frame(Fields@data, mean=meancalc) 

For confirmation

> str(Fields2$mean)
 num [1:1293] 275 1350 442 430 1364 ...
  • 1
    I would suspect something is not quite right with the polygons (eg they are encoded clockwise instead of anticlockwise, or something), but its impossible to tell without the data. If you want identical behaviour in R as QGIS, then use the gdal_rasterize function from the gdalUtils package, which calls the same command-line function as listed in the QGIS dialog. – Spacedman Apr 24 at 11:00
  • 2
    I noticed something unusual with the min and max values of the SpatialPolygonsDataFrame. Why is min area > max area, and min mean > max mean? – Micha Apr 24 at 11:01
  • Well spotted - I wonder if these are character values? Try Fields2$mean2 = as.numeric(as.character(Fields2$mean)) and rasterizing with mean2... – Spacedman Apr 24 at 11:06
  • Note that instead of printing an sp-object if you do library(raster) and summary(sp_object) it will tell you the class of the column and you should see immediately if its not numeric. – Spacedman Apr 24 at 11:11
  • I have edited my question to add more information. True, the min/max value is strange to me. mean is numeric although I have tried the option as.numeric to double check and the same result is obtained. Any way to share the shapefile to have a closer look? – GCGM Apr 24 at 12:42

The last feature in your shapefile is a bit.... off - it looks like a union of all the polygon lines:

> plot(Fields2[1293,])

enter image description here

I'm pretty sure that's messing up rasterize. Not sure how QGIS is coping, but I reckon if you can figure out where this has come from and remove it then you'll stop having problems in R.

Hmmm actually it looks like the polygon of the negative space that isn't the first 1292 polygons within the rectangle, but its convoluted form perhaps means rasterize isn't sure what's the inside and what's the outside...

And indeed if you read the shapefile in with sf::st_read and test with st_is_valid you get to see that the last feature is not valid.

  • You are right, Fields2[1293,] comes from a union operation that I did in QGIS. I had those field boundaries and I need a polygon to be created in the empty spaces in between. For that I used the shapefile of the study area and the union operator – GCGM May 7 at 8:37

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.