Grass module for Effective Mesh Size, CBC

effemeshsize_cbc.R

# Effective Mesh Size (CBC)

library(vroom)
library(dplyr)

# pour surface en ha
divisurf <- 10000

nom_colonnes <- c("id", "id_patch", "surf_m2", "nb_pixels")
tbl_stats_grass <- vroom("FR/stats_completes.txt", delim = "\t", col_names= nom_colonnes, col_types = "iidi")

# surface en a, ou ha, ou km2 (selon divisurf)
tbl_stats_grass <- mutate(tbl_stats_grass, a = surf_m2 / divisurf)
tbl_stats_grass

# calculer tableau avec surf complete des patchs
# remarque : enlever id_patch=NA (pixels hors patches) avant de calculer sum 
tbl_patch_cmpl <- tbl_stats_grass %>% 
  filter(!is.na(id_patch)) %>%
  group_by(id_patch) %>% 
  summarise(a_cmpl = sum(a), nb_cells = n())
tbl_patch_cmpl

# calculer tableau avec surf totale des mailles (patch et non patch)
# remarque : enlever id=NA (pixels hors de la grille) 
tbl_msize  <- tbl_stats_grass %>%
  filter(!is.na(id)) %>%
  group_by(id) %>% 
  summarise(aT = sum(a))
tbl_msize

# jointure aI-acpmlI + calcul produit
tbl_stats_2 <- tbl_stats_grass %>% 
  filter(!is.na(id) & !is.na(id_patch)) %>%
  left_join(tbl_patch_cmpl, by="id_patch")
tbl_stats_2 <-  mutate(tbl_stats_2, a2=a*a_cmpl)
tbl_stats_2

# calcule sommpe de produits (aires interieur / complet)
tbl_stats_grid <- tbl_stats_2 %>% group_by(id) %>% 
  summarise(sum_a2 = sum(a2, na.rm = TRUE), 
            patches=n(), sum_aI = sum(a))
tbl_stats_grid

# calculer eff mesh size
tbl_msize_2 <- left_join(tbl_msize, tbl_stats_grid, by="id")
tbl_msize_2 <- mutate(
  tbl_msize_2, sum_a2=coalesce(sum_a2, 0), 
  patches=coalesce(patches, 0), sum_aI=coalesce(sum_aI, 0))
tbl_msize_2 <- mutate(tbl_msize_2, cbc_msize=sum_a2/aT, divisurf)
tbl_msize_2

vroom_write(tbl_msize_2, "effmeshsize_cbc.txt")

library(sf)
sf_grid <- st_read("FR/grid-L93.shp")
sf_grid_2 <- left_join(sf_grid, tbl_msize_2)
st_write(sf_grid_2, "FR/grid_cbc_msize.shp", delete_dsn=TRUE)

effmeshsize_cbc.md

Workflow Calcul de l'indice "Effective Mesh Size CBC"

Fichiers en entrée :

`OCS_2018_CESBIO.tif` = Couche occupation du sol CESBIO
`grid-WGS84.shp` = grille WGS84 
`emprise_fr.shp` = emprise zone analyse (pour masquer raster)

Import données

• import raster OCS_2018_CESBIO.tif (SCR=EPSG:2154) => OCS_2018_CESBIO
• import vecteur grid-WGS84.shp (SCR=EPSG:4326) => grid_L93 (SCR=EPSG:2154)
• import vecteur emprise_fr.shp => emprise_5km

Traitement grille vecteur

• ajoute colonne id_int (INT)
• MAJ colonne id_int=id

Définir zone analyse

• zone = emprise_5km alignée sur OCS_2018_CESBIO

Convertir grille vecteur en raster

• conversion grid_L93 => grid_rast
• conversion emprise_5km => rast_emprise_5km (masque pour le calcul à venir) avec val=1

Filtrer OCS zones naturelles (Calculatrice raster)

• si OCS dans (13,18,19) et rast_emprise_5km = Vrai alors OCS_cat131819_null = Vrai sinon OCS_cat131819_null = null
• exporter OCS_cat131819_null => result131819_clip.tif (COMPRESS=DEFLATE,PREDICTOR=2)

Identifier patches

• identifier les patchs dans OCS_cat131819_null => OCS_cat131819_clump

Stats patches X mailles

• superposer grid_rast X OCS_cat131819_clump, pour compter le nombre de pixels de chaque patch dans chaque maille de la grille => stats_completes.txt

Calculer Effective Mesh Size

• lire stats_completes.txt dans un tableau
• calculer a = surface patch X maille dans le tableau tbl_stats_grass
• calculer a_cmpl = surface totale des patchs (grouper par id_patch puis somme) dans le tableau tbl_patch_cmpl
• calculer aT = surface totale des mailles (grouper par id puis somme) dans le tableau tbl_msize
• joindre tbl_stats_grass > tbl_patch_cmpl
• calculer a2 = produit a * a_cmpl
• calculer (somme des a2) / aT
• joindre le résultat à la grille vectorielle pour visualiser l'indice de fragmentation

Script shell

    r.in.gdal -o input=D:\Travail\Support\Leandro\OCS_2018_CESBIO.tif output=OCS_2018_CESBIO
    v.in.ogr -o input=D:\Travail\Support\Leandro\FR\emprise_fr.shp output=emprise_5km

    # importer grille vecteur en lambert 93 => nelle couche vecteur grid_L93
    v.import input='C:\Support\Leandro\grid-WGS84.shp' layer=grid-WGS84 output=grid_L93
    # ajoute champ id_int INTEGER
    v.db.addcolumn map=grid_L93@PERMANENT columns='id_int INT'
    v.db.update map=grid_L93@PERMANENT column=id_int query_column=id

    g.region -s -p vector=emprise_5km@PERMANENT align=OCS_2018_CESBIO@PERMANENT
    v.to.rast input=grid_L93@PERMANENT type=area output=grid_rast use=attr attribute_column=id_int
    v.to.rast input=emprise_5km@PERMANENT type=area output=rast_emprise_5km use=val

    r.mapcalc --overwrite expression=OCS_cat131819 = (OCS_cat131819 = OCS_2018_CESBIO@PERMANENT == 13 || OCS_2018_CESBIO@PERMANENT == 18 || OCS_2018_CESBIO@PERMANENT == 19) && rast_emprise_5km@PERMANENT == 1
    r.out.gdal input=OCS_cat131819@PERMANENT output=D:\Travail\Support\Leandro\result131819_clip.tif format=GTiff type=Byte createopt=COMPRESS=DEFLATE,PREDICTOR=2 nodata=255
    r.mapcalc --overwrite expression=OCS_cat131819_null = OCS_cat131819@PERMANENT == 1 ? 1 : null()
    r.clump -d input=OCS_cat131819_null@PERMANENT output=OCS_cat131819_clump
    r.stats -a -c input=grid_rast@PERMANENT,OCS_cat131819_clump@PERMANENT output=C:\Support\Leandro\FR\stats_completes.txt separator=tab null_value=NA

module_OCS.py

#!/usr/bin/env python3
#
#########################################################################
#
# MODULE:       Effective Mesh Size
#
# AUTHOR(S):    Bernard C
#
# PURPOSE:      Module description
#
# DATE:         Mon Sep  6 18:26:24 2021
#
#########################################################################

#%module
#% description: Effective Mesh Size, CBC
#% keyword: landscape analysis
#% keyword: fragmentation
#%end

import sys
import os
import atexit

import grass.script as gscript

RAST_REMOVE = []

def cleanup():
    gscript.run_command('g.remove', flags='f', type='raster',
                          name=RAST_REMOVE)

def main():
    options, flags = gscript.parser()

        # donnees en entree
    v_grid = 'grid_L93@PERMANENT'
    v_emprise = 'emprise_5km@PERMANENT'
    r_ocs = 'OCS_2018_CESBIO@PERMANENT' 


    # definir region
    gscript.run_command('g.region', vector=v_emprise, align=r_ocs, flags='ps')
    # rasterisation
    gscript.run_command('r.mask', vector=v_emprise)
    gscript.run_command(
        'v.to.rast', input=v_grid, type='area', output='grid_rast', 
        use='attr', attribute_column='id_int')
    #gscript.run_command('v.to.rast', input=v_buffer_routes, type='area', output='rast_emprise_5km', use='val')
    # mapcalc
    exp_calc = f'OCS_cat131819 = if(({r_ocs} == 13 || {r_ocs} == 18 || {r_ocs} == 19), 1, null())'
    gscript.mapcalc(exp_calc)
    gscript.run_command(
        'r.clump', input='OCS_cat131819', output='OCS_cat131819_clump', flags='d')
    result = gscript.read_command(
        'r.stats', input=['grid_rast', 'OCS_cat131819_clump'], separator='pipe', null_value='NA', flags='ac')

    gscript.run_command('g.remove', flags='f', type='raster',
                        name=RAST_REMOVE)

    return 0

if __name__ == "__main__":
    atexit.register(cleanup)
    sys.exit(main())

script_import.py

#!/usr/bin/env python3

# IMPORT DES DONNEES DANS LA BASE 
# - COUCHE RASTER OCS_2018_CESBIO
# - COUCHE VECTOR grid_L93
# - COUCHE VECTOR emprise_5km

import grass.script as gscript

def main():

    # donnees en entree
    shp_grid = 'D:/misc/Leandro/shp/grid-WGS84.shp'
    shp_emprise = 'D:/misc/Leandro/shp/emprise_fr.shp'
    tif_ocs = 'D:/misc/Leandro/tif/OCS_2018_CESBIO.tif' 

    # couches GRASS en sortie
    v_grid = 'grid_L93'
    v_emprise = 'emprise_5km'
    r_ocs = 'OCS_2018_CESBIO' 

    # ecraser les données existantes ?
    ecraser = True    

    # COUCHE RASTER OCS_2018_CESBIO # flags='o' pour passer outre la reprojection
    gscript.run_command('r.in.gdal', input=tif_ocs, output=r_ocs, flags='o', overwrite=ecraser)
    # COUCHE VECTOR grid_L93 # reprojection WGS84 -> L93
    gscript.run_command('v.import', input=shp_grid, output=v_grid, overwrite=ecraser)
    # COUCHE VECTOR emprise_5km # flags='o' pour passer outre la reprojection
    gscript.run_command('v.import', input=shp_emprise, output=v_emprise, flags='o', overwrite=ecraser)

if __name__ == '__main__':
    main()

script_OCS.py

#!/usr/bin/env python3
#
# script: EffMeshSize

import os
import grass.script as gscript

TMP_MAPS = []

def cleanup():

    # effacer couches temporaires
    gscript.run_command('g.remove', flags='f', type='raster', name=TMP_MAPS)

def main():

    # donnees en entree : vecteur grille + raster OCS + colonne id grille
    v_grid = 'grid_L93@PERMANENT'
    r_ocs = 'OCS_2018_CESBIO@PERMANENT' 
    nom_colonne_id = 'id_int'

    # codes zone naturelle
    zones_nat = [13,18,19]

    # optionnel : couche vecteur masque (delimite zone a traiter)
    v_mask = 'emprise_5km@PERMANENT'

    # donnees en sortie (temporaire)
    tmp_rast_grid = 'grid_rast_' 
    tmp_rast_mask = 'emprise_rast_'
    tmp_rast_grid_mask = 'grid_rast_mask_'
    tmp_rast_zonesnat = 'ocs_znat_'
    tmp_rast_patches = 'ocs_znat_patches_'
    # tmp_rast_grid = 'grid_rast_' + str(os.getpid())
    # tmp_rast_mask = 'emprise_rast_' + str(os.getpid())
    # tmp_rast_grid_mask = 'grid_rast_mask_' + str(os.getpid())
    # tmp_rast_zonesnat = 'ocs_znat_' + str(os.getpid())
    # tmp_rast_patches = 'ocs_znat_patches_' + str(os.getpid())
    TMP_MAPS = [tmp_rast_grid, tmp_rast_mask, tmp_rast_zonesnat]

    # definir region temporaire
    gscript.run_command('g.region', vector=v_mask, align=r_ocs, flags='ps')

    # rasterisation grille
    gscript.run_command(
        'v.to.rast', input=v_grid, type='area', output=tmp_rast_grid, 
        use='attr', attribute_column=nom_colonne_id)

    # rasterisation emprise
    if (v_mask != None):
        gscript.run_command(
            'v.to.rast', input=v_mask, type='area', output=tmp_rast_mask, use='val', value=1)
        # application du masque sur le raster grille 
        formule_calc = tmp_rast_grid_mask + ' = if(' + tmp_rast_mask + '==1, ' + tmp_rast_grid + ', null())'
        gscript.mapcalc(formule_calc)
        # critere additionnel masque pour masquer OCS
        crit_mask = ' && ' + tmp_rast_mask + '==1'
    else:
        gscript.run_command('g.copy', raster=[tmp_rast_grid, tmp_rast_grid_mask])
        crit_mask = ''

    # creer raster zones naturelles avec mapcalc
    # exemple de formule : 
    # "ocs_znat = if((OCS_2018_CESBIO == 13 || OCS_2018_CESBIO == 18 || OCS_2018_CESBIO == 19) && emprise_rast==1, 1, null())"
    liste_crit = [r_ocs + ' == ' + str(x) for x in zones_nat]
    formule_calc = tmp_rast_zonesnat + ' = if((' + (' || '.join(liste_crit)) + ')' + crit_mask + ', 1, null())'
    gscript.mapcalc(formule_calc)
    gscript.run_command('r.clump', input=tmp_rast_zonesnat, output=tmp_rast_patches, flags='d')
    cleanup()

if __name__ == '__main__':
    gscript.use_temp_region()
    main()

script_OCS.R

library(terra)
terraOptions(datatype="INT2U")

# lire donnees OCS, grille, zone etude
rast_ocs <- rast("tif/OCS_2018_CESBIO.tif")
vect_grd <- vect("shp/grid-WGS84.shp")
vect_emp <- vect("shp/emprise_dep11.shp")

# harmoniser sys coord
ref_crs <- crs(rast_ocs)
vect_grd <- project(vect_grd, ref_crs)
vect_emp <- project(vect_emp, ref_crs)

# crop + mask
ext_analyse <- ext(vect_emp)
rast_ocs <- crop(rast_ocs, ext_analyse)
rast_ocs <- mask(rast_ocs, vect_emp)

# raster zones naturelles
rast_znat <- rast_ocs %in% c(13, 18, 19)

# zone nat par maille
#df_agreg <- extract(rast_znat, vect_grd, fun=sum)
rast_grd <- rasterize(vect_grd, rast_ocs, field="id")
df_agreg2 <- zonal(x=rast_znat, z=rast_grd, fun=sum)


# detecte patches
rast_znat_patches <- patches(rast_znat, directions=8, zeroAsNA=TRUE)

# stats par maille