How can I obtain cell values and coordinate data from a (.tif) raster when R rasterToPoints function-CodePudding

I am interested in extracting the cell values alongside their corresponding x and y coordinates from a tif file accessible from the WorldPop database [ https://hub.worldpop.org/geodata/summary?id=49920 ].

I have converted this file alongside other tif files available on this website into rasters and then used the rasterToPoints function in R to extract this information. However, although this approach has worked for most of the files, it has failed for this particular file amongst a few others. It's like the R session remains stuck and the code just never runs when I attempt to convert the rasters to spdf data.

library(raster)
Raster <- raster("C:/file path/aus_ppp_2020_UNadj_constrained.tif")
Raster <- rasterToPoints(Raster, spatial = TRUE)

As an alternative, I thought I could extract the coordinates after obtaining the cell values using getValues() or readAll() but due to the size of the raster being too large I run into the following error:

Error: cannot allocate vector of size 17.8 Gb.

sessionInfo()
# R version 4.2.0 (2022-04-22 ucrt)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 22000)

library(memuse)
memuse::Sys.meminfo()
# Totalram:  31.781 GiB 
# Freeram:   26.164 GiB

I then tried to see if I could modify the usable memory using memory.limit() but this code has been retired from R version 4.2 and I cannot find an alternative.

memory.limit() 
# Warning: 'memory.limit()' is no longer supported[1] Inf

I was wondering if anyone knows:

1. If there is a way I could get the rasterToPoints function to work for this raster.

2. If there is a way to subset the raster to smaller rasters, while retaining all data, so that I could use the rasterToPoints function on each subset and then merge the resulting spatial point dataframes.

3. If there is an alternative way to extract the coordinates alongside the cell values for this tif file.

Any help is greatly appreciated.

CodePudding user response：

First step would be to question why you are converting the raster to a data frame in the first place. There is a reason R can handle the full raster but not all the points as a data frame because a raster is a very efficient way of working with data.

To chop the raster up hopefully this method will work:

library(raster)
library(sf)
library(dplyr)

## download raster from source
r <-  raster("https://data.worldpop.org/GIS/Population/Global_2000_2020_Constrained/2020/BSGM/AUS/aus_ppp_2020_UNadj_constrained.tif")

## determine the extent of the raster layer
x_min <- r@extent@xmin
x_max <- r@extent@xmax
y_min <- r@extent@ymin
y_max <- r@extent@ymax

## create a df
d <- data.frame(X = c(x_min, x_max, x_max, x_min), 
                     Y = c(y_max, y_max, y_min, y_min))

## generate an sf of the raster domain
d_sf <- d %>%
  st_as_sf(coords = c("X", "Y"), crs = 4326) %>%
  dplyr::summarise(geometry = st_combine(geometry)) %>% 
  st_cast("POLYGON")

## define the grid resolution (i.e. the size of the grids the domain will be divided into)
res <- 10

## create the grid
d_grid <- st_make_grid(d_sf, square = T, n = c(res, res), crs = 4326) %>% # the grid, covering bounding box
  st_sf()

## pick a cell
random_cell <- d_grid[44,]
  
## crop the raster to this cell
  r_crop <- crop(r, random_cell)

  ## check you have data
plot(r_crop)

## raster to points
pts <- rasterToPoints(r_crop)

## to trim the data to a specific location can use the osmdata package to get the extent of a place and crop to that
library(osmdata)

## go for Melbourne
p <- osmdata::getbb("Melbourne")

r_place_crop <- crop(r, p)
plot(r_place_crop)
place_pts <- rasterToPoints(r_place_crop)

CodePudding user response：

Although it seems like an overall questionable idea to break an efficient storage apart in general - it's not like the relation between values and coordinates is lost in a raster - there is a pretty convenient way to do so using terra.

Import raster data, convert raster to points, convert points to data frame with coordinates listed:

library(terra)
#> terra 1.6.7

r <- system.file("ex/elev.tif", package="terra") |> rast()

df <- as.points(r) |> as.data.frame(geom = "XY") 

head(df)
#>   elevation        x        y
#> 1       529 6.004167 50.17917
#> 2       542 6.012500 50.17917
#> 3       547 6.020833 50.17917
#> 4       535 6.029167 50.17917
#> 5       485 5.970833 50.17083
#> 6       497 5.979167 50.17083

However, with your raster aus_ppp_2020_UNadj_constrained.tif, as.points(r) seems to be failing with Error in x@ptr$as_points(values, na.rm, na.all, opt) : std::bad_alloc, but it seems like you can work around this using as.polygons(r) in a preceding call:

library(terra)
#> terra 1.6.7

r <- rast("aus_ppp_2020_UNadj_constrained.tif")

df <- as.polygons(r) |> as.points() |> as.data.frame(geom = "XY") 

head(df)
#>   aus_ppp_2020_UNadj_constrained        x          y
#> 1                              0 143.3087  -9.862916
#> 2                              0 143.3087  -9.863750
#> 3                              0 143.3096  -9.863750
#> 4                              0 143.3096  -9.862916
#> 5                              0 143.3212 -10.104583
#> 6                              0 143.3212 -10.106250

However, be aware of the size of this monstrosity you created with approx. 5 million rows.