3

I have a large raster file (5GB) containing only 1's and NA's. I would like to convert this into a multipolygon of the areas with 1's, with adjacent cells dissolved into one polygon.

I have imported the file to R using

r = raster::raster(my_filename)
r
class      : RasterLayer 
dimensions : 17452, 45000, 785340000  (nrow, ncol, ncell)
resolution : 0.008, 0.008  (x, y)
extent     : -180, 180, -55.9875, 83.6285  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : C://...binary_X01_januarysnow.asc 
names      : binary_X01_januarysnow

and I have tried several methods to create the polygon:

  • rasterToPolygons from raster with dissolve==TRUE option (R crashes)
  • isoband from the isoband package (R crashes),

Both of the approaches have worked as expected when I've tried them on a subset of my raster covering appr. the area of Spain, so I assume the problem is only with the size of the data and not my code.

  • Then I have tried to read my raster with read_stars, and use stars::st_as_sf(st, as_points = FALSE, merge = TRUE, connect8 = TRUE). This returned an empty polygon, possibly because the file was read as a stars proxy object, but I'm not sure, I couldn't find any information about that online. Then I have force-read the raster as stars and not as stars proxy by using read_stars(my_filename, proxy=FALSE) and have tried to use the st_as_sf command as above but got the message "Error: cannot allocate vector of size 2.9 Gb"

I know that in the worst case I can probably just decrease the raster resolution and therefore size and will be able to create the polygons I want (but with less precise resolution), but I was wondering if anyone has another suggestion I could try? Both the 1's and NA's are located in large continous areas, so it would be enough to have high resolution on the edges, if that helps.

PS This is my very first question on StackOverflow so I apologize if my problem is not clearly described. I don't know how to provide a reproducible example of a large dataset.

Berbara Boulos
  • 33
  • 1
  • 4
  • 1
    If I understand you correctly, your raster is 70,000 x 70,000 pixels in size. For an area the size of Spain, that means you have about 1 pixel per 10 m. Do you really need this level of detail? Even reducing this by a factor of 10 in each dimension wouldn't make any difference on the final plotted object on a high resolution monitor. Also, if there are multiple polygons, can't you split the raster into regions around each polygon? – Allan Cameron Oct 05 '20 at 14:23
  • I have edited the question to show more info about the raster. The map is global, I have just tried out the rasterToPolygons and the isoband on the Spain area to see if my approach works at all. I'm doing it for a colleague and they have chosen to work with this resolution, yes, if there is no other way we will have to decrease it. I also thought about splitting up the global map into a couple of regions and make polygons separately within each. This is quite fiddly but could be an option. – Berbara Boulos Oct 05 '20 at 15:23

1 Answers1

5

What you are looking for is as.polygons() from the terra package, the raster package's successor. terra handles large data sets better than raster does.

Chr
  • 1,017
  • 1
  • 8
  • 29