RASTERIZE A SPATVECTOR MAKES R SESSION TO ABORT

[anjelinejeline]

Hi I am using terra 1.7-46
I am trying to rasterize a large point dataset (worldwide coverage) but when I run the function rasterize the R session aborts for a fatal error.
I would like to have a raster with cell value equals to the number of points
Could you please fix this bug?

Thank you !

df<-read.csv("df.csv")[,-1]

df<-st_as_sf(df,coords=c("x", "y"), crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0")

data(World)
World<-World |> 
  filter(continent %in% c("Asia","Africa","Europe","Oceania","South America","North America")) |> 
  st_transform("+proj=moll +lon_0=0 +x_0=0 +y_0=0")

# Create a raster template of the World with 1kmx1km resolution

r <- terra::rast(World)
terra::res(r) <- 1000


Vect<-terra::vect(df)

Raster<-terra::rasterize(Vect, r, fun="count", field="ID")

[anjelinejeline]

If I change the above code with the following it runs but an empty raster is returned

r <- terra::rast(World,ncols=1000,nrows=1000,
                 crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0")

Vect<-terra::vect(df)

Raster<-terra::rasterize(Vect,r, fun="count",field="ID")

terra::plot(Raster)

[Rapsodia86]

Is "World" a point dataset?
If so, use "length" instead of "counts"; see: https://rdrr.io/cran/terra/man/rasterize.html

fun - summarizing function for when there are multiple geometries in one cell. For lines and polygons you can only use "min", "max", "mean", "count" and "sum" For points you can use any function that returns a single number; for example mean, length (to get a count), min or max

[rhijmans]

I cannot help with this because I do not have df nor World. And if I make a similar example that works fine.

library(terra)
r <- project(rast(ymin=-60), "+proj=moll +lon_0=0 +x_0=0 +y_0=0")
# res(r) <- 1000
# easier to see the results with a larger resolution
res(r) <- 500000

s <- seq(-500000, 5000000, 1000)
df <- data.frame(x=s, y=s)
v <- vect(df, c("x", "y"))
x <- terra::rasterize(v, r, fun="count")

[anjelinejeline]

Hi @rhijmans @Rapsodia86 the World shapefile is from the package tmap.
The df is point dataframe in sf format

library("tmap") # World shapefile
data(World)
World <- World[World$continent %in% c("Asia","Africa","Europe","Oceania","South America","North America"), ]
World <- st_transform(World, "+proj=moll +lon_0=0 +x_0=0 +y_0=0")

r <- terra::rast(World, resolution=1000)

The idea is to have a grid of **1kmx1km **covering the entire World and then count the points falling in each grid...
I also tried with rasterizeGeom
but once again the R session aborts by itself...
Can you help me understand what I am doing wrong, is the way I create the raster of the World fine? and if the terra package can work with such a big grid of the world?

Vect<-terra::vect(df)

Raster<-terra::rasterizeGeom(Vect, r,
                                      fun="count",unit="m",
                                      resolution=c(1000,1000))

Thank you
Angela

[rhijmans]

I have already shown above how you can do this (in a simpler way). If that does not work for you, it may be related to "df". But I do not have that, and you do not even show what it looks like. With the information at hand I cannot say what is going on.

[anjelinejeline]

Hi @rhijmans ,
Please find attached df...
It's quite big but I believe that the terra package should be handle it .. so I don't know why the rasterization does not work
df.csv

[HassanMasoomi]

It works fine for me ("df.csv" is the one you provided above)...

library(sf)
library(terra)
library(tmap)
library(data.table)

# sample raster
data(World)
World <- World[World$continent %in% c("Asia","Africa","Europe","Oceania","South America","North America"),]
World <- sf::st_transform(World, "+proj=moll +lon_0=0 +x_0=0 +y_0=0")
sample_raster <- terra::rast(World, resolution = 1000)

# assuming that x & y provided in df are already in the same crs as sample_raster
df <- fread("df.csv")
df_vect <- vect(df, geom = c("x", "y"), 
                crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0")
final_raster <- terra::rasterize(df_vect, sample_raster, fun = "count")

# another method (maybe better if your df is huge?!!); result will be on memory!
# 1: find the cell each row in located in
df[, cell := cellFromXY(sample_raster, .SD[, .(x, y)])]
# 2: summarize df by cell
new_df <- df[, .(count = .N), keyby = cell]
# 3: set each cell values in your raster
final_raster2 <- sample_raster
final_raster2[new_df$cell] <- new_df$count

@rhijmans , regarding the last line in the code above, is there any function to set the values of a bunch of cells for a raster (a function with filename argument)?

[anjelinejeline]

Hi @HassanMasoomi
Thank you !!
I was able to create it using the first method .. the R session did not abort this time around.
If I look at the statistics using global the values seem fine but if I try to plot it the plot is empty .. I believe that it may be due to the fine resolution..
Have you tried to plot it ?
I'd like to visualize the results somehow ..
Any advice on how to do that ?

Thank you

[HassanMasoomi]

I don't know your goal of using the raster with that high resolution here. But, maybe you want to explore tmap and plot you results as points:

library(sf)
library(terra)
library(tmap)
library(data.table)

# sample raster
data(World)
World <- World[World$continent %in% c("Asia","Africa","Europe","Oceania","South America","North America"),]
World <- sf::st_transform(World, "+proj=moll +lon_0=0 +x_0=0 +y_0=0")
sample_raster <- terra::rast(World, resolution = 1000)

# assuming that x & y provided in df are already in the same crs as sample_raster
df <- fread("df.csv")

# 1: find the cell each row in located in
df[, cell := cellFromXY(sample_raster, .SD[, .(x, y)])]
# 2: summarize df by cell
new_df <- df[, .(count = .N), keyby = cell]

# plot
new_df[, x := xFromCell(sample_raster, cell)]
new_df[, y := yFromCell(sample_raster, cell)]
new_df_vect <- vect(new_df, geom = c("x", "y"),
                    crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0")
tmap_mode("view")
tm_shape(sf::st_as_sf(new_df_vect)) + tm_dots(size = "count", col = "darkred", 
                                              alpha = 0.75, id = "count")

[anjelinejeline]

@HassanMasoomi
Thank you but I would like to plot it as a raster grid .. is there no way I can visualise it as a raster ?

[anjelinejeline]

Hi @rhijmans it now works ... I changed the machine .. I guess that the crash was due to my laptop ..
I have just another question ..

If I want to assign to each cell of my sample raster [created with sample_raster=rast(World, resolution = 1000)]
not the count of points of a spatial vector x but rather the sum of an attribute (myattribute) of all points within that specific raster cell ..should I do like this:

# Sum the variable of interest
  ras=rasterize(x,sample_raster,field="myattribute", fun="sum")

I tried and it seems that it works but I would like to double check with you whether this is fine...

Best wishes
Angela