Lib_Validation_biodivMapR.R 15.5 KB
Newer Older
jbferet's avatar
jbferet committed
1
# ==============================================================================
Florian de Boissieu's avatar
Florian de Boissieu committed
2
3
# biodivMapR
# Lib_ValidationbiodivMapR.R
jbferet's avatar
jbferet committed
4
5
6
7
8
9
10
# ==============================================================================
# PROGRAMMERS:
# Jean-Baptiste FERET <jb.feret@irstea.fr>
# Copyright 2018/07 Jean-Baptiste FERET
# ==============================================================================
# This Library contains functions to manipulate & process raster images
# Mainly applicable to ENVI HDR data wth BIL interleave
11
# ==============================================================================
jbferet's avatar
jbferet committed
12
13

#' get projection of a raster or a vector
De Boissieu Florian's avatar
De Boissieu Florian committed
14
15
#' @param file path for a raster or vector (shapefile)
#' @param type 'raster' or 'vector'
16
#' @return projection
De Boissieu Florian's avatar
De Boissieu Florian committed
17
#' @importFrom raster raster shapefile projection
18
#' @importFrom rgdal readOGR
19
#' @import tools
20
#' @export
De Boissieu Florian's avatar
De Boissieu Florian committed
21
projection.file <- function(file, type = 'raster'){
jbferet's avatar
jbferet committed
22
  if (Type == 'raster'){
De Boissieu Florian's avatar
De Boissieu Florian committed
23
    obj <- raster(file)
jbferet's avatar
jbferet committed
24
  } else if (Type == 'vector'){
De Boissieu Florian's avatar
De Boissieu Florian committed
25
26
    Shp.Path      = dirname(file)
    Shp.Name      = file_path_sans_ext(basename(file))
jbferet's avatar
jbferet committed
27
28
    obj  = readOGR(dsn = Shp.Path,layer = Shp.Name,verbose = FALSE)
  }
De Boissieu Florian's avatar
De Boissieu Florian committed
29
30
  projstr <- projection(obj)
  return(projstr)
jbferet's avatar
jbferet committed
31
32
}

33
34
#' Get list of shapefiles in a directory
#'
De Boissieu Florian's avatar
De Boissieu Florian committed
35
#' @param x character or list. Directory containing shapefiles
jbferet's avatar
jbferet committed
36
#' @return list of shapefiles names
37
#' @export
De Boissieu Florian's avatar
De Boissieu Florian committed
38
39
list.shp <- function(x){
  if(typeof(x)=='list'){
jbferet's avatar
jbferet committed
40
41
    List.Shp = c()
    ii = 0
De Boissieu Florian's avatar
De Boissieu Florian committed
42
    for (shp in x){
jbferet's avatar
jbferet committed
43
      ii = ii+1
De Boissieu Florian's avatar
De Boissieu Florian committed
44
      List.Shp[ii]  = dir(shp, pattern = '.shp$', full.names = TRUE, ignore.case = FALSE,include.dirs = FALSE)
jbferet's avatar
jbferet committed
45
46
    }
  } else if(typeof(Path.Shp)=='character'){
De Boissieu Florian's avatar
De Boissieu Florian committed
47
    List.Shp  = dir(Path.Shp, pattern = '.shp', full.names = TRUE, ignore.case = FALSE,include.dirs = FALSE)
jbferet's avatar
jbferet committed
48
49
50
51
  }
  return(List.Shp)
}

52
53
54
55
56
57
# reprojects a vector file and saves it
# @param Initial.File path for a shapefile to be reprojected
# @param Projection projection to be applied to Initial.File
# @param Reprojected.File path for the reprojected shapefile
# @return
#' @importFrom rgdal readOGR
58
#' @import tools
jbferet's avatar
jbferet committed
59
Reproject.Vector = function(Initial.File,Projection,Reprojected.File){
60

jbferet's avatar
jbferet committed
61
62
63
64
  Shp.Path      = dirname(Initial.File)
  Shp.Name      = file_path_sans_ext(basename(Initial.File))
  Vect.Init     = readOGR(Shp.Path,Shp.Name,verbose = FALSE)
  Proj.init     = projection(Vect.Init)
65

jbferet's avatar
jbferet committed
66
67
68
69
70
71
72
73
74
75
  if (!Proj.init==Projection){
    Shp.Path      = dirname(Reprojected.File)
    Shp.Name      = file_path_sans_ext(basename(Reprojected.File))
    Vect.reproj = spTransform(Vect.Init, Projection)
    writeOGR(obj = Vect.reproj, dsn = Shp.Path,layer = Shp.Name, driver="ESRI Shapefile",overwrite_layer = TRUE) #also you were missing the driver argument
  }
  return()
}

# Extracts pixels coordinates from raster corresponding to an area defined by a vector
76
77
78
79
# @param Path.Raster path for the raster file. !! BIL expected
# @param Path.Vector path for the vector file. !! SHP expected
# @return ColRow list of coordinates of pixels corresponding to each polygon in shp
#' @importFrom rgdal readOGR
80
#' @import tools
jbferet's avatar
jbferet committed
81
82
83
84
85
86
87
88
89
Extract.Pixels.Coordinates = function(Path.Raster,Path.Vector){
  # read vector file
  Shp.Path    = dirname(Path.Vector)
  Shp.Name    = file_path_sans_ext(basename(Path.Vector))
  Shp.Crop    = readOGR(Shp.Path,Shp.Name)
  # read raster info
  Raster      = raster(Path.Raster, band = 1)
  # extract pixel coordinates from raster based on vector
  XY  = raster::cellFromPolygon (Raster,Shp.Crop)
90
  # for each polygon in the
jbferet's avatar
jbferet committed
91
92
93
94
95
96
97
  ColRow = list()
  for (i in 1:length(XY)){
    ColRow[[i]] = ind2sub(Raster,XY[[i]])
  }
  return(ColRow)
}

98
99
100
101
# Extracts pixels coordinates from raster corresponding to an area defined by a vector
# @param Path.Raster path for the raster file. !! BIL expected
# @param OGR.Vector  OGR for the vector file obtained from readOGR
# @return ColRow list of coordinates of pixels corresponding to each polygon in shp
jbferet's avatar
jbferet committed
102
103
104
105
106
107
108
109
110
Extract.Pixels.Coordinates.From.OGR = function(Path.Raster,OGR.Vector){
  # read raster info
  Raster      = raster(Path.Raster, band = 1)
  # for each polygon or point in the shapefile
  ColRow = list()
  # extract pixel coordinates from raster based on vector
  if (OGR.Vector@class[1]=='SpatialPointsDataFrame'){
    XY  = raster::cellFromXY (Raster,OGR.Vector)
    ColRow[[1]] = ind2sub(Raster,XY)
111

jbferet's avatar
jbferet committed
112
113
114
115
116
117
118
119
120
  } else if  (OGR.Vector@class[1]=='SpatialPolygonsDataFrame'){
    XY  = raster::cellFromPolygon (Raster,OGR.Vector)
    for (i in 1:length(XY)){
      ColRow[[i]] = ind2sub(Raster,XY[[i]])
    }
  }
  return(ColRow)
}

121
122
123
124
# This function reads information directly from a zipfile
# @param Zip.Path path for the zipped binary image file in .BIL format and
# @param HDR Header corresponding
# @return Img.Data vecotr containing all image data
jbferet's avatar
jbferet committed
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
Read.Image.Zipfile = function(Zip.Path){
  # get HDR corresponding to zipfile
  ImPathHDR       = Get.HDR.Name(Zip.Path)
  HDR             = read.ENVI.header(ImPathHDR)
  nb.Elements     = HDR$samples * HDR$lines * HDR$bands
  Image.Format    = ENVI.Type2Bytes(HDR)
  fid             = unz(description = Zip.Path, filename = zip_list(Zip.Path)$filename, open = "rb", encoding = getOption("encoding"))
  if (Image.Format$Type=='INT'){
    Img.Data      = readBin(fid, integer(), n = nb.Elements, size = Image.Format$Bytes, endian = "little")
  } else if (Image.Format$Type=='FLOAT'){
    Img.Data      = readBin(fid, numeric(), n = nb.Elements, size = Image.Format$Bytes, endian = "little")
  }
  Img.Data        = aperm(array(Img.Data,dim=c(HDR$samples,HDR$bands,HDR$lines)),c(3,1,2))
  Img.Data         = c(array(Img.Data[,,HDR$bands],dim=c(HDR$lines,HDR$samples,length(HDR$bands))))
  close(fid)
  return(Img.Data)
}

143
144
145
# Computes alpha diversity metrics from distribution
# @param Distrib distribution of clusters
# @return Richness, Fisher, Shannon, Simpson
jbferet's avatar
jbferet committed
146
Get.Alpha.Metrics = function(Distrib){
147
  RichnessPlot  = vegan::specnumber(Distrib, MARGIN = 1)        # species richness
jbferet's avatar
jbferet committed
148
149
150
151
152
153
  # if (length(Distrib)>1){
  #   fisherPlot    = fisher.alpha(Distrib, MARGIN = 1)      # fisher's alpha
  # } else {
  #   fisherPlot    = 0
  # }
  FisherPlot    = 0
154
155
156
  ShannonPlot   = vegan::diversity(Distrib, index = "shannon", MARGIN = 1, base = exp(1)) # shannon's alpha
  SimpsonPlot   = vegan::diversity(Distrib, index = "simpson", MARGIN = 1, base = exp(1))
  return(list("Richness" = RichnessPlot,"fisher"=FisherPlot,"Shannon"=ShannonPlot,"Simpson"=SimpsonPlot))
jbferet's avatar
jbferet committed
157
158
159
160
161
162
}

#' gets alpha diversity indicators from plot
#' @param Raster SpectralSpecies file computed from DiverstyMapping method
#' @param Plots list of shapefiles included in the raster
#' @return alpha and beta diversity metrics
163
#' @importFrom rgdal readOGR
164
#' @import tools
165
#' @export
De Boissieu Florian's avatar
De Boissieu Florian committed
166
diversity_from_plots = function(Raster, Plots,NbClusters = 50,Name.Plot = FALSE){
jbferet's avatar
jbferet committed
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
  # get hdr from raster
  HDR           = read.ENVI.header(paste(Raster,'.hdr',sep=''))
  nbRepetitions = HDR$bands
  # get the number of plots
  nbPlots             = length(Plots)
  # define variables where alpha diversity will be stored
  Richness.AllRep = Shannon.AllRep = Fisher.AllRep = Simpson.AllRep = list()
  # initialize alpha diversity for each plot
  Richness       = Shannon    = Fisher     = Simpson    = data.frame()
  Name.Vector       = list()
  # prepare directory to write shapefiles with correct projection
  Dir.Vector            = dirname(Plots[[1]])
  Dir.Vector.reproject  = paste(Dir.Vector,'/Reproject',sep='')
  ###                 Compute alpha diversity             ###
  # for each plot
  Pixel.Inventory.All   = list()
  PolygonID             = 0
184

jbferet's avatar
jbferet committed
185
186
187
188
189
190
  for (ip in 1:nbPlots){
    # prepare for possible reprojection
    File.Vector           = Plots[[ip]]
    Name.Vector[[ip]]     = file_path_sans_ext(basename(File.Vector))
    print(paste('Processing file ',Name.Vector[[ip]],sep=''))
    File.Vector.reproject = paste(Dir.Vector.reproject,'/',Name.Vector[[ip]],'.shp','sep'='')
191

jbferet's avatar
jbferet committed
192
193
194
    if (file.exists(paste(file_path_sans_ext(File.Vector),'.shp',sep=''))){
      Plot                = readOGR(Dir.Vector,Name.Vector[[ip]],verbose = FALSE)
      # check if vector and rasters are in teh same referential
195
      Projection.Plot     = projection.file(File.Vector,'vector')
jbferet's avatar
jbferet committed
196
197
198
199
200
201
202
203
204
205
206
      # if not, convert vector file
      if (!Projection.Raster==Projection.Plot){
        if (!file.exists(File.Vector.reproject)){
          dir.create(Dir.Vector.reproject, showWarnings = FALSE,recursive = TRUE)
          Reproject.Vector(File.Vector,Projection.Raster,File.Vector.reproject)
        }
        Plot              = readOGR(Dir.Vector.reproject,Name.Vector[[ip]],verbose = FALSE)
      }
    } else if (file.exists(paste(VectorFile,'kml','sep'='.'))){
      print('Please convert vector file to shpfile')
    }
207

jbferet's avatar
jbferet committed
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
    # extract data corresponding to the raster
    XY          = Extract.Pixels.Coordinates.From.OGR(Raster,Plot)
    # if the plot is included in the raster
    if (length(XY)==1 & length(XY[[1]]$Column)==0){
      if (length(Name.Plot)==nbPlots){
        Name.Plot[ip] = NA
      }
    }
    if (length(XY)>1 | length(XY[[1]]$Column)>0){
      ID.Poly     = list()
      idPix       = 1
      if (length(XY)==1){
        coordPix    = cbind(XY[[1]]$Row,XY[[1]]$Column)
        ID.Poly[[1]]= seq(idPix,idPix+length(XY[[1]]$Row)-1,by = 1)
      } else {
        coordPix    = list()
        ii          = 0
        for (pp in 1:length(XY)){
          ii        = ii+1
          if (length(XY[[pp]]$Row)>0){
            coordPix[[pp]]  = cbind(XY[[pp]]$Row,XY[[pp]]$Column)
            ID.Poly[[pp]]   = seq(idPix,idPix+length(XY[[pp]]$Row)-1,by = 1)
            idPix           = idPix+length(XY[[pp]]$Row)
          }
        }
        coordPix    = do.call(rbind,coordPix)
      }
      ExtractIm   = Extract.Pixels(coordPix,Raster,HDR)
      # compute alpha diversity for each repetition
      Pixel.Inventory = list()
      Richness.tmp = Shannon.tmp = Fisher.tmp = Simpson.tmp = vector(length = nbRepetitions)
      # eliminate spectral species contributing to less than pcelim percent of the total valid pixels
      pcelim  = 0.02
      for (ii in 1:length(ID.Poly)){
        if (length(ID.Poly[[ii]])>0){
          for (i in 1:nbRepetitions){
            Distritab     = table(ExtractIm[ID.Poly[[ii]],i])
            Pixel.Inventory[[i]]  = as.data.frame(Distritab)
            SumPix        = sum(Pixel.Inventory[[i]]$Freq)
            ThreshElim    = pcelim*SumPix
            ElimZeros = which(Pixel.Inventory[[i]]$Freq<ThreshElim)
            if (length(ElimZeros)>=1){
              Pixel.Inventory[[i]]   = Pixel.Inventory[[i]][-ElimZeros,]
            }
            if (length(which(Pixel.Inventory[[i]]$Var1==0))==1){
              Pixel.Inventory[[i]]   = Pixel.Inventory[[i]][-which(Pixel.Inventory[[i]]$Var1==0),]
            }
            Alpha = Get.Alpha.Metrics(Pixel.Inventory[[i]]$Freq)
            # Alpha diversity
            Richness.tmp[i]   = as.numeric(Alpha$Richness)
            Fisher.tmp[i]     = Alpha$fisher
            Shannon.tmp[i]    = Alpha$Shannon
            Simpson.tmp[i]    = Alpha$Simpson
          }
          PolygonID = PolygonID+1
          Richness.AllRep[[PolygonID]] = Richness.tmp
          Shannon.AllRep[[PolygonID]]  = Shannon.tmp
          Fisher.AllRep[[PolygonID]]   = Fisher.tmp
          Simpson.AllRep[[PolygonID]]  = Simpson.tmp
          Richness   = rbind(Richness, mean(Richness.tmp), row.names = NULL, col.names = NULL)
          Fisher     = rbind(Fisher, mean(Fisher.tmp), row.names = NULL, col.names = NULL)
          Shannon    = rbind(Shannon, mean(Shannon.tmp), row.names = NULL, col.names = NULL)
          Simpson    = rbind(Simpson, mean(Simpson.tmp), row.names = NULL, col.names = NULL)
          Pixel.Inventory.All[[PolygonID]] = Pixel.Inventory
        }
      }
    }
  }
  Richness.AllRep = do.call(rbind,Richness.AllRep)
  Shannon.AllRep  = do.call(rbind,Shannon.AllRep)
  Fisher.AllRep   = do.call(rbind,Fisher.AllRep)
  Simpson.AllRep  = do.call(rbind,Simpson.AllRep)
280

jbferet's avatar
jbferet committed
281
282
283
284
285
286
287
288
289
290
  ###                 Compute beta diversity             ###
  # for each pair of plot, compute beta diversity indices
  BC = list()
  for(i in 1:nbRepetitions){
    MergeDiversity = matrix(0,nrow = NbClusters,ncol = PolygonID)
    for(j in 1:PolygonID){
      SelSpectralSpecies  =  as.numeric(as.vector(Pixel.Inventory.All[[j]][[i]]$Var1))
      SelFrequency        = Pixel.Inventory.All[[j]][[i]]$Freq
      MergeDiversity[SelSpectralSpecies,j] = SelFrequency
    }
291
    BC[[i]] = vegan::vegdist(t(MergeDiversity),method="bray")
jbferet's avatar
jbferet committed
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
  }
  BC_mean   = 0*BC[[1]]
  for(i in 1:nbRepetitions){
    BC_mean = BC_mean+BC[[i]]
  }
  if (length(Name.Plot)>1){
    elim = which(is.na(Name.Plot))
    if (length(elim)>0){
      Name.Plot =   Name.Plot[-elim]
    }
  }
  BC_mean = as.matrix(BC_mean/nbRepetitions)
  return(list("Richness" = Richness,"Fisher"=Fisher,"Shannon"=Shannon,"Simpson"=Simpson,'BCdiss' = BC_mean,"fisher.All"=Fisher.AllRep,"Shannon.All"=Shannon.AllRep,"Simpson.All"=Simpson.AllRep,'BCdiss.All' = BC,'Name.Plot'=Name.Plot))
}

307
308
309
310
311
312
313
314
315
316
# build a vector file from raster footprint
# borrowed from https://johnbaumgartner.wordpress.com/2012/07/26/getting-rasters-into-shape-from-r/
# @param x path for a raster or raster object
# @param outshape path for a vector to be written
# @param gdalformat
# @param pypath
# @param readpoly
# @param quiet
# @return NULL
#' @importFrom rgdal readOGR
317
#' @importFrom raster writeRaster
jbferet's avatar
jbferet committed
318
319
gdal_polygonizeR = function(x, outshape=NULL, gdalformat = 'ESRI Shapefile',
                            pypath=NULL, readpoly=TRUE, quiet=TRUE) {
320

jbferet's avatar
jbferet committed
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
  if (is.null(pypath)) {
    pypath <- Sys.which('gdal_polygonize.py')
  }
  if (!file.exists(pypath)) stop("Can't find gdal_polygonize.py on your system.")
  owd <- getwd()
  on.exit(setwd(owd))
  setwd(dirname(pypath))
  if (!is.null(outshape)) {
    outshape  = sub('\\.shp$', '', outshape)
    f.exists  = file.exists(paste(outshape, c('shp', 'shx', 'dbf'), sep='.'))
    if (any(f.exists)){
      print('Footprint of raster file already defined')
    }
    # stop(sprintf('File already exists: %s',
    #              toString(paste(outshape, c('shp', 'shx', 'dbf'),
    #                             sep='.')[f.exists])), call.=FALSE)
  } else outshape <- tempfile()
  if (is(x, 'Raster')) {
    writeRaster(x, {f <- tempfile(fileext='.tif')})
    rastpath <- normalizePath(f)
  } else if (is.character(x)) {
    rastpath <- normalizePath(x)
  } else stop('x must be a file path (character string), or a Raster object.')
  outshape.Shp = paste(outshape,'.shp',sep='')
  if (!file.exists(outshape.Shp)){
    system2('python', args=(sprintf('"%1$s" "%2$s" -f "%3$s" "%4$s.shp"',
                                    pypath, rastpath, gdalformat, outshape)))
  }
  if (readpoly==TRUE) {
    shp <- readOGR(dirname(outshape), layer = basename(outshape), verbose=!quiet)
    return(shp)
  }
  return(NULL)
}

356
357
358
359
360
# create a binary mask file based on a matrix and the corresponding header
# @param Mask matrix of a binary mask
# @param HDR header corresponding to an image to be masked
# @param Path.Mask path for the mask
# @return
jbferet's avatar
jbferet committed
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
Create.Shademask=function(Mask,HDR,Path.Mask){
  ipix        = HDR$lines
  jpix        = HDR$samples
  Mask        = array(Mask,c(ipix,jpix,1))
  Mask        = aperm(Mask,c(2,3,1))
  fidOUT      = file(description = Path.Mask, open = "wb", blocking = TRUE,
                     encoding = getOption("encoding"), raw = FALSE)
  writeBin(c(as.integer(Mask)), fidOUT, size = 1,endian = .Platform$endian, useBytes = FALSE)
  close(fidOUT)
  # write updated shademask
  HDR.Update             = HDR
  HDR.Update$bands       = 1
  HDR.Update$`data type` = 1
  HDR.Update$`band names`= {'Mask'}
  HDR.Update$wavelength  = NULL
  HDR.Update$fwhm        = NULL
  HDR.Update$resolution  = NULL
  HDR.Update$bandwidth   = NULL
  HDR.Update$purpose     = NULL
  HDRFpath = paste(Path.Mask,'.hdr',sep='')
  write.ENVI.header(HDR.Update, HDRFpath)
  return()
}