rspatial / raster

R raster package https://rspatial.github.io/raster/reference/raster-package.html
GNU General Public License v3.0
161 stars 53 forks source link

Mask function fail inside a for loop #328

Closed rondon-d closed 8 months ago

rondon-d commented 8 months ago

hello all, I am having a strange issue with a mask function in raster in R. I want to crop the data of my raster file to a Simple feature that has a hexagon as a geometry type. I need to crop the data for 300 polygons and count how many characteristics are found in each one.

My R code looks as follow:

library(sp)
library(sf)
library(s2)
library(dplyr) 
library(raster)
library(tidyverse)
library(stars)

load(file="Cov_Data.Rdata")
# print some data
str(DataArea)
Hex.G<-Hex.G$x

  Shape_1<-raster("Clc2018_FI20m.tif")

  Hex.G <- st_as_sf(Hex.G)
  Hex.G <- st_transform(Hex.G, st_crs(Shape_1)) 

  Data_1<-NULL
  Data_1<- list()

  Aux_2<-Aux_3<-NULL
  i<-1

  for(i in 1:300){
    print(i)
    Aux_2<-mask(Shape_1,  
                Hex.G[i,] 
                )   
    Aux_3<-crop(Aux_2, 
                extent(Hex.G[i,]) )  
    Data_1[[i]]<-table(values(Aux_3))  
    Aux_2<-Aux_3<-NULL
  }

Data looks as follow:

Shape_1 class : RasterLayer dimensions : 61978, 35989, 2230526242 (nrow, ncol, ncell) resolution : 20, 20 (x, y) extent : 20000, 739780, 6597200, 7836760 (xmin, xmax, ymin, ymax) crs : +proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs source : Clc2018FI20m.tif names : Count values : 1, 49 (min, max) attributes : ID Count_ LEVEL4 LEVEL3 LEVEL2 LEVEL1 Level4Suo Level4Eng from: 1 427379 1111 111 11 1 Kerrostaloalueet Continuous urban fabric to : 49 130493524 5231 523 52 5 Meret Sea and ocean

Hex.G Simple feature collection with 313 features and 0 fields Geometry type: POLYGON Dimension: XY Bounding box: xmin: 63747.59 ymin: 6613938 xmax: 763747.6 ymax: 7803279 Projected CRS: +proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs First 10 features: x 1 POLYGON ((83747.59 6683220,... 2 POLYGON ((103747.6 6648579,... 3 POLYGON ((123747.6 6683220,... 4 POLYGON ((143747.6 6648579,... 5 POLYGON ((163747.6 6613938,... 6 POLYGON ((163747.6 6683220,... 7 POLYGON ((183747.6 6648579,... 8 POLYGON ((183747.6 6717861,... 9 POLYGON ((183747.6 6787143,... 10 POLYGON ((183747.6 6856425,...

The code works properly if i extract the information individually (not in a for loop). Inside a for usually works for 2 or 3 loops and then it provides the following error:

Error in result[, i] <- readBin(raster@file@con, what = dtype, n = ncols, : replacement has length zero In addition: There were 50 or more warnings (use warnings() to see the first 50)

Thank in advance!

rondon-d commented 8 months ago

Hello, I just my code to Terra and seems to be working with out issues. :)

rhijmans commented 8 months ago

You can probably reduce your script to:

library(terra)
x <- rast("Clc2018_FI20m.tif")
hex <- vect(st_as_sf(Hex.G)) |> project(x)
d <- extract(x, hex, fun="table")