2

I'm trying to obtain a GeoTIFF map of the WWF biomes, using their shapefile download: https://www.worldwildlife.org/publications/terrestrial-ecoregions-of-the-world
The purpose of this is to eventually use a raster object to extract biome types at species' occurrences to find out what biomes those species live in. There are two 'layers': biomes which are the larger areas, and ecoregions, which are subdivisions of those biomes. I'm not interested in the ecoregions.

The zip file download includes the following files:

 [1] "layerfiles_readme.txt" "terr_biomes.lyr"       "wwf_terr_ecos.dbf"     "wwf_terr_ecos.htm"    
 [5] "wwf_terr_ecos.lyr"     "wwf_terr_ecos.prj"     "wwf_terr_ecos.sbn"     "wwf_terr_ecos.sbx"    
 [9] "wwf_terr_ecos.shp"     "wwf_terr_ecos.shp.xml" "wwf_terr_ecos.shx"

Something doesn't work because the raster object and GeoTIFF file I created is a 1 x 1 pixel image. This is the code and output from my attempt at obtaining a GeoTIFF file out of the shapefile:

> shapefile1 <- readOGR(getwd(), "wwf_terr_ecos")
OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\horse\OneDrive\Desktop\official\official", layer: "wwf_terr_ecos"
with 14458 features
It has 21 fields
> 
> extent(shapefile1) 
class      : Extent 
xmin       : -180 
xmax       : 180 
ymin       : -89.89197 
ymax       : 83.62313 
> 
> summary(shapefile1) 
Object of class SpatialPolygonsDataFrame
Coordinates:
         min       max
x -179.99999 179.99999
y  -89.89197  83.62313
Is projected: FALSE 
proj4string : [+proj=longlat +datum=WGS84 +no_defs]
Data attributes:
    OBJECTID          AREA           PERIMETER         ECO_NAME            REALM          
 Min.   :    1   Min.   :      0   Min.   :  0.000   Length:14458       Length:14458      
 1st Qu.: 3630   1st Qu.:      7   1st Qu.:  0.140   Class :character   Class :character  
 Median : 7256   Median :     28   Median :  0.321   Mode  :character   Mode  :character  
 Mean   : 7279   Mean   :  10191   Mean   :  3.682                                        
 3rd Qu.:10910   3rd Qu.:    252   3rd Qu.:  1.078                                        
 Max.   :14925   Max.   :9111849   Max.   :815.023                                        
     BIOME           ECO_NUM         ECO_ID         ECO_SYM          GBL_STAT     G200_REGIO       
 Min.   : 1.000   Min.   : 0.0   Min.   :-9999   Min.   : 12.00   Min.   :0.00   Length:14458      
 1st Qu.: 1.000   1st Qu.: 3.0   1st Qu.:40146   1st Qu.: 61.00   1st Qu.:1.00   Class :character  
 Median : 5.000   Median : 9.0   Median :51115   Median : 73.00   Median :2.00   Mode  :character  
 Mean   : 7.218   Mean   :13.7   Mean   :51611   Mean   : 97.08   Mean   :1.83                     
 3rd Qu.:11.000   3rd Qu.:19.0   3rd Qu.:61402   3rd Qu.:111.00   3rd Qu.:3.00                     
 Max.   :99.000   Max.   :82.0   Max.   :81333   Max.   :887.00   Max.   :3.00                     
    G200_NUM          G200_BIOME       G200_STAT        Shape_Leng         Shape_Area      
 Min.   :-9999.00   Min.   : 0.000   Min.   :0.0000   Min.   :  0.0123   Min.   :   0.000  
 1st Qu.:    0.00   1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:  0.1403   1st Qu.:   0.001  
 Median :    0.00   Median : 0.000   Median :0.0000   Median :  0.3212   Median :   0.003  
 Mean   :  -36.36   Mean   : 2.628   Mean   :0.5381   Mean   :  3.6824   Mean   :   1.482  
 3rd Qu.:   72.00   3rd Qu.: 3.000   3rd Qu.:1.0000   3rd Qu.:  1.0792   3rd Qu.:   0.028  
 Max.   :  147.00   Max.   :14.000   Max.   :3.0000   Max.   :814.8107   Max.   :4936.890  
    area_km2         eco_code            PER_area          PER_area_1   PER_area_2
 Min.   :      0   Length:14458       Min.   : 0.00000   Min.   :0    Min.   :0   
 1st Qu.:  17859   Class :character   1st Qu.: 0.00000   1st Qu.:0    1st Qu.:0   
 Median :  77060   Mode  :character   Median : 0.00000   Median :0    Median :0   
 Mean   : 206948                      Mean   : 0.03675   Mean   :0    Mean   :0   
 3rd Qu.: 172182                      3rd Qu.: 0.00000   3rd Qu.:0    3rd Qu.:0   
 Max.   :4629416                      Max.   :17.62724   Max.   :0    Max.   :0   
> 
> r <- raster(extent(shapefile1))
> projection(r) <- proj4string(shapefile1)
Warning message:
In proj4string(shapefile1) :
  CRS object has comment, which is lost in output
>     
> res(r)=1000
> 
> r2 <- rasterize(shapefile1, r, "BIOME")
> 
> writeRaster(r2, filename="biomes.tif")
Vince
  • 20,017
  • 15
  • 45
  • 64
cflax
  • 21
  • 2
  • Your shapefile uses a longlat projection, with units of degrees. When you set res(r)=1000, you are telling R that you want your raster to have raster cells that are 1000 units on a side. If follows your instructions, producing a single raster cell that's as close as possible to 1000 degrees on a side (while still fitting into the raster's extent). Presumably, you actually want a raster with much smaller cells, like 1 degree or 0.1 degrees on a side, so adjust your call to res(r) <- accordingly. – Josh O'Brien Oct 12 '20 at 15:29
  • Why do you have to go through the step of "rasterizing" the biomes? You can perform a simple point-in-polygon or, if all of the data are polygons, an intersection function. This will be much more computationally efficient as well as save the introduction of uncertainty. There are multiple functions that facilitate this analysis and many of them have been covered on this forum eg., https://gis.stackexchange.com/questions/310372/joinn-shapefile-attribute-id-to-point-in-r/310377#310377 – Jeffrey Evans Oct 14 '20 at 15:05

1 Answers1

1

The easiest way of rasterizing that shapefile's biome field to a global raster longlat layer of a 1 x 1 degree resolution is probably the following:

# Load packages
packs <- list("tidyverse", "raster", "sf", "fasterize")
lapply(packs, require, character.only = T)

Load shapefile


wwf_shp <- st_read("C:/Users/horse/OneDrive/Desktop/official/official/wwf_terr_ecos.shp", stringsAsFactors = F)


Rasterize biome field and write to disk


raster(resolution = 1, crs = "+proj=longlat +datum=WGS84", vals = 0) %>% 
   fasterize(wwf_shp, ., field = "BIOME") %>% 
   writeRaster(., "C:/Users/horse/OneDrive/Desktop/official/official/biomes.tif")

You might also want to specify a function in fasterize() via the fun argument, which is similar to the fun argument in rasterize().

Using fasterize::fasterize() instead of raster::rasterize() has the advantages of being a lot faster - especially on large data sets - and being compatible with the sf package, the sp package's successor. Alternatively, you can use the raster package's successor terra. It comes with its own vector data format, SpatVector, and is also much faster than raster.

user
  • 298
  • 1
  • 3
  • 8