ghislainv / forestatrisk

📦🐍 Python package to model and forecast the risk of deforestation
https://ecology.ghislainv.fr/forestatrisk
GNU General Public License v3.0
116 stars 27 forks source link

Forest cover input #66

Closed ricardozwarg closed 2 years ago

ricardozwarg commented 2 years ago

Hi @ghislainv

I'm quite confused about the forest cover map that is requested as input for the function far.sample. In the documentation I found the following information associate to it: param input_forest_raster: Name of the forest raster file (1=forest, 0=deforested) in the var_dir directory..

Considering it, I'm assuming for instance that 1 should be equal to the remaining forest, and 0 to the amout of deforestation observed during the time period considered. However, I am not sure how to generate this kind of input, since in my head I will need to assign null values for all other pixels that are not remaing forest (class 1) or even deforestation during the period of interest (class 0).

So, the first question is, how to generate this kind of input? using GEE or Python.

Nevertheless, I tried to running the model with two explanation vaiables (SRTM and slope) to see what will happen using as input a raster file with 1 been the remaining forest, and 0 all the rest, because I was not capable to eliminate the pixels that are not deforestation during the period analysed. As a result I got an error, that I will be sharing below. My guess is that this observed error is associate with my input.

I also sharing my progress with a link for google earth engine, and also information about the data I'm using, as well the error message that I got trying to run the package with jupyter notebook.

Thank you in advance @ghislainv

Best from Brazil,

Ricardo

Link to the script on GEE: https://code.earthengine.google.com/d505cf980983d44be7efe4b7d8517dc6 Link to the deforestation asset: https://code.earthengine.google.com/?asset=users/Jurupari/PDigital2000_2021_AMZ Link to the study area asset: https://code.earthengine.google.com/?asset=users/Jurupari/RRL

error observed.txt

PDigital2000_2021_AMZ.txt

ghislainv commented 2 years ago

Dear @ricardozwarg,

As input, you need a forest cover change raster with values 0 (deforestation between date 1 and date 2) and 1 (remaining forest at date 2). All other pixels should be set to any other value (eg., nodata value 255 for a Byte type raster).

I would advise you to have a close look at the Get Started tutorial of forestatrisk. In the data folder, you can inspect raster fcc23.tif:

gdalinfo -stats data/fcc123.tif

You can obtain such raster with any GIS software (GDAL, QGis, GRASS) or as follow using GEE. Pay attention to the following point: in GEE, the mask is set by default for values different from zero but this can be changed (see below).

// Importing study area
var rrl = ee.FeatureCollection('users/Jurupari/RRL');

// Importing deforestation data from PRODES (Official Monitoring Project in the Legal Amazon)
var prodes = ee.Image('users/Jurupari/PDigital2000_2021_AMZ').clip(rrl);

// Getting loss for the period of interest (2009 - 2019);
var historicalLoss = prodes.gte(8).and(prodes.lte(28)).selfMask();
Map.addLayer(historicalLoss, {palette: 'red'}, 'Forest loss (2009-2019)');

// Getting forest in 2009
var forest2009 = prodes.gte(8).and(prodes.lte(29)).add(prodes.eq(31)).add(prodes.eq(33)).add(prodes.eq(1)).rename('forest2009').selfMask();
Map.addLayer(forest2009, {palette: 'green'}, 'Forest in 2009');

// Getting forest in 2019
var forest2019 = forest2009.where(historicalLoss, 0).rename('forest2019').selfMask();
Map.addLayer(forest2019, {palette: ('green')}, 'Forest in 2019');

// Forest cover change
var fcc = ee.Image.constant(255).clip(rrl);
var fcc = fcc.where(forest2009.eq(1), 0).where(forest2019.eq(1), 1);
var mask = fcc.neq(255);
var fcc = fcc.updateMask(mask);
Map.addLayer(fcc, {min:0, max:1, palette: ['red', 'green']}, 'fcc');

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: fcc.unmask(255),
  description: 'fcc',
  scale: 30,
  region: rrl
});

Optionally, set 255 as nodata value for the raster fcc.tif with gdal_translate locally:

gdal_translate -a_nodata 255 fcc.tif fcc_nodata.tif -co COMPRESS=DEFLATE -co PREDICTOR=2

Best regards,

Ghislain

ricardozwarg commented 2 years ago

Thank you @ghislainv

I have corrected the forest input map using the code above and after it the error in the forestatrisk package have been solved.

I have another question about the package, but I will be posting it by opening a new issue, okay?

All the best,

Ricardo

ghislainv commented 2 years ago

Hi @ricardozwarg, OK, great if you manage to solve your problem with the input file. Yes, please, open a new issue for a new problem. Regards, Ghislain