Search code examples
r-rasterrgdalrasterize

Rasterizing a list of shps inR


I've been trying to rasterize a list of shp imported in R. There is an error when I call the list because it is not recognize as dataframe but as a character. I've tryed to use as.dataframe(list[i]) without success.

Error: Unable to find an inherited method for function ‘rasterize’ for signature ‘"character", "RasterLayer"’

## import Shp
shp_ELC <- readOGR("EconomicLandConsesions/ELC.shp")
shp_CF <- readOGR("Communityforestry/Community_forestryPolygon.shp")
shp_NPA <- readOGR("Natural-Protected-Areas/Natural Protected Areas.shp")
shp_Ecoregion <- readOGR("TerrestrialEcoRegionsWWF/KH_TerrEcoregions_Dslv.shp")

## create Raster template
utm     <- proj4string(shp_CF)
my_proj <- spTransform(shp_Ecoregion,utm)
my_ext  <- extent(my_proj)
my_res  <- 30
tempR    <- raster(resolution = my_res, ext = my_ext, crs = my_proj)
## list SHP
listShp <- ls(pattern = glob2rx("*.df"))

## looping list
for(i in 1:length(listShp)){
  rst_CF  <- rasterize(data.frame(listShp[i]),tempR)
}

Solution

  • Here are some notes

    readOGR would not work with the arguments you provided. The shapefile function would. I changed the variable names because after reading the data into R, they are no longer shapefiles. Presumably the are SpatialPolygonsDataFrame objects

    library(raster)
    p_ELC <- shapefile("EconomicLandConsesions/ELC.shp")
    p_CF <- shapefile("Communityforestry/Community_forestryPolygon.shp")
    p_NPA <- shapefile("Natural-Protected-Areas/Natural Protected Areas.shp")
    p_Ecoregion <- shapefile("TerrestrialEcoRegionsWWF/KH_TerrEcoregions_Dslv.shp")
    
    # you can create a list right away, but with only a few objects that might make things harder for you
    # f <- c("EconomicLandConsesions/ELC.shp", "Communityforestry/Community_forestryPolygon.shp", "Natural-Protected-Areas/Natural Protected Areas.shp", "TerrestrialEcoRegionsWWF/KH_TerrEcoregions_Dslv.shp")
    # pList <- lapply(f, shapefile)
    

    Apparently some of the polygons have a UTM CRS, while others do not. Make sure that they are all transformed to the same CRS.

    utmcrs  <- crs(p_CF)
    p_Ecoregion <- spTransform(p_Ecoregion, utmcrs)
    
    # create Raster template
    r <- raster(resolution=30, ext = extent(p_Ecoregion), crs = utmcrs)
    
    # make a list of the SpatialPolygonDataFrame objects
    listP <- list(p_ELC, p_CF, p_NPA, p_Ecoregion)
    

    I am not sure what to do now. In your example you overwrite the output three times. I assume you want four rasters

    x <- list()
    for(i in 1:length(listP)){
      x[[i]] <- rasterize(listP[[i]]), r)
    }
    
    s <- stack(x)
    plot(s)