R/iucn_spdf_to_raster.R

In harithmorgadinho/wege:

#' @title spdf_to_raster
#' @description A function rasterize iucn shapefiles.
#' @param input  df or raster
#' @param crs projection of data provided ('longlat'/'cea'/'auto')
#' @param zoom_out zoom level for the map
#' @param r_extent an object to use as extent
#' @examples
#' crop_map_world(df)
#' crop_map_world(df,crs='longlat')
#' @export
#'@importFrom sp spTransform
#'@importFrom raster raster
#'@importFrom raster rasterize

ggg_check <- function() {
  if (!requireNamespace('ggg', quietly = TRUE)) {
    msg <- paste0('Package ', 'ggg', ' required. Run ',
                  'devtools::install_github(\'harithmorgadinho/ggg\')',
                  'or similar.')
    stop(msg, call. = FALSE)
  }
}


iucn_spdf_to_raster=function(spdf,r_extent='auto',resolution=NULL,taxa=NULL,output_crs=NULL,filter_taxa='yes',buffer=0){
  ggg_check()
  if(filter_taxa=='yes'){
    cat('loading list of species')
    flpth <- system.file('terrestrial_animals_status.RDS', package = 'wege')
    load(flpth)
  }
  cat('guessing projection of spdf')
  crs_input=ggg::crs_guesser(spdf)


  if(is.null(output_crs)){
    output_crs=crs_input
  }

  if(output_crs=='longlat'){
    if (crs_input=='cea') {
      spdf=spdf
    }
    if(crs_input=='longlat'){
      spdf=spTransform(spdf,CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
    }
  }


  if(output_crs=='cea'){
    if (crs_input=='cea') {
      spdf=spdf
    }
    if(crs_input=='longlat'){
      cat('transforming spdf to cea projection')
      spdf=spTransform(spdf,CRS("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
    }
  }


  crs_input=crs_guesser(spdf)

  if(is.null(r_extent)){

    r_extent=extent(spdf)

  }

  if (is.null(resolution)){
    cat('resolution set to 1','\n')
    resolution=1
  }

  if (crs_input=='cea') {
    resolution=resolution*100000
    r=raster(extent(r_extent), resolution=resolution,crs=CRS("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))

  }

  if(crs_input=='longlat'){

    r=raster(extent(r_extent), resolution=resolution,crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

  }
  if (filter_taxa=='yes'){

    if (taxa=='amphibians'){

      list_rank=amphibian_status

    }

    if (taxa=='mammals'){

      list_rank=mammal_status

    }
    if (taxa=='birds'){

      list_rank=bird_status

    }

    if (taxa=='reptiles'){

      list_rank=reptile_status

    }

    sp_list=list_rank$species
    organismgroup=sp_list
    cat('subsetting iucn shp to selected taxa','\n')
    Ipol2_subset=spdf[spdf$BINOMIAL %in% organismgroup, ]


  }
  if (filter_taxa=='no'){
    Ipol2_subset=spdf
    organismgroup=unique(spdf$BINOMIAL)
    Ipol2_subset=spdf
  }
  ###overlay gbif records
  #gbif_data_amphibia
  IUCN_raster_final=r
  values(IUCN_raster_final)=0


  organismgroup = unique(Ipol2_subset$BINOMIAL)


  buffer=resolution*buffer


  r_grid=rasterToPolygons(r)

  for (i in seq_along(organismgroup)){
    organism <- organismgroup[[i]]
    cat(organism,'[', i, '/', length(organismgroup), ']\n',
        sep = '')



    IUCN_spec <- Ipol2_subset[(Ipol2_subset@data[,"BINOMIAL"]==organism),]
    IUCN_spec=gBuffer(IUCN_spec,width=0,byid=TRUE)
    #plot(IUCN_spec)
    #Rasterize IUCN polygon
    vals <- over(r_grid, IUCN_spec)
    if(nrow(vals)==1 && is.na(vals[,1])){
      cat(organism,'not present','\n')
    }
    else{


      res <- cbind(coordinates(r_grid), vals)
      res=res[!is.na(res$ID),]
      df=res[,c(5,2,1)]
      df_coord=df[,3:2]
      spdf <- SpatialPointsDataFrame(df_coord, df,proj4string =
                                       CRS('+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'))
      IUCN_raster=rasterize(spdf, r, field=1,
                            getCover=T,
                            buffer=buffer)



      #IUCN_raster <- rasterize(IUCN_spec, r, field=1, getCover=T,buffer=buffer)

      #IUCN_raster[IUCN_raster >0] = 1
      #IUCN_raster[IUCN_raster ==0] = NA
      #plot(IUCN_raster)
      #rasterize points

      IUCN_raster[is.na(IUCN_raster)]=0
      IUCN_raster_final=overlay(IUCN_raster_final,IUCN_raster,fun=sum)
      #plot(IUCN_raster_final)

    }
  }
  return(IUCN_raster_final)
}