ghislainv / forestatrisk

📦🐍 Python package to model and forecast the risk of deforestation

forestatrisk Python package

Overview

The forestatrisk Python package can be used to model the tropical deforestation spatially, predict the spatial risk of deforestation, and forecast the future forest cover in the tropics. It provides functions to estimate the spatial probability of deforestation as a function of various spatial explanatory variables.

Spatial explanatory variables can be derived from topography (altitude, slope, and aspect), accessibility (distance to roads, towns, and forest edge), deforestation history (distance to previous deforestation), or land conservation status (eg. protected area) for example.

Scientific publication

Vieilledent G. 2021. forestatrisk: a Python package for modelling and forecasting deforestation in the tropics. Journal of Open Source Software. 6(59): 2975. [doi: 10.21105/joss.02975].

Statement of Need

Spatial modelling of the deforestation allows identifying the main factors determining the spatial risk of deforestation and quantifying their relative effects. Forecasting forest cover change is paramount as it allows anticipating the consequences of deforestation (in terms of carbon emissions or biodiversity loss) under various technological, political and socio-economic scenarios, and informs decision makers accordingly. Because both biodiversity and carbon vary greatly in space, it is necessary to provide spatial forecasts of forest cover change to properly quantify biodiversity loss and carbon emissions associated with future deforestation.

The forestatrisk Python package can be used to model the tropical deforestation spatially, predict the spatial risk of deforestation, and forecast the future forest cover in the tropics. The spatial data used to model deforestation come from georeferenced raster files, which can be very large (several gigabytes). The functions available in the forestatrisk package process large rasters by blocks of data, making calculations fast and efficient. This allows deforestation to be modeled over large geographic areas (e.g. at the scale of a country) and at high spatial resolution (eg. ≤ 30 m). The forestatrisk package offers the possibility of using logistic regression with auto-correlated spatial random effects to model the deforestation process. The spatial random effects make possible to structure the residual spatial variability of the deforestation process, not explained by the variables of the model and often very large. In addition to these new features, the forestatrisk Python package is open source (GPLv3 license), cross-platform, scriptable (via Python), user-friendly (functions provided with full documentation and examples), and easily extendable (with additional statistical models for example). The forestatrisk Python package has been used to model deforestation and predict future forest cover by 2100 across the humid tropics (https://forestatrisk.cirad.fr).

Installation

You will need several dependencies to run the forestatrisk Python package. The best way to install the package is to create a Python virtual environment, either through conda (recommended) or virtualenv.

Using virtualenv

The easiest way to install the forestatrisk Python package is via pip in the OSGeo4W Shell for Windows or in a virtual environment for Linux.

For Linux, create and activate a virtual environment before installing geefcc with pip:

cd ~
# Create a directory for virtual environments
mkdir venvs
# Create the virtual environment with venv
python3 -m venv ~/venvs/venv-geefcc
# Activate (start) the virtual environment
source ~/venvs/venv-geefcc/bin/activate

Install Python dependencies and forestatrisk in the OSGeo4W Shell or in the newly created virtual environment:

# Upgrade pip, setuptools, and wheel
python3 -m pip install --upgrade pip setuptools wheel
# Install numpy
python3 -m numpy
# Install gdal Python bindings (the correct version)
python3 -m pip install gdal==$(gdal-config --version)
# Install forestatrisk. This will install all other dependencies
python3 -m pip install forestatrisk

If you want to install the development version of forestatrisk, replace the last line with:

python3 -m pip install https://github.com/ghislainv/forestatrisk/archive/master.zip

To deactivate and delete the virtual environment:

deactivate
rm -R ~/venvs/venv-forestatrisk # Just remove the repository

In case of problem while installing GDAL Python bindings, try the following command:

python3 -m pip install --no-cache-dir --force-reinstall gdal==$(gdal-config --version)

Using conda

You first need to have miniconda3 installed (see here).

Then, create a conda environment (details here) and install the forestatrisk package with the following commands:

conda create --name conda-far -c conda-forge python gdal numpy matplotlib pandas patsy pip statsmodels earthengine-api --yes
conda activate conda-far
pip install pywdpa scikit-learn # Packages not available with conda
pip install forestatrisk # For PyPI version
# pip install https://github.com/ghislainv/forestatrisk/archive/master.zip # For GitHub dev version
# conda install -c conda-forge python-dotenv --yes  # Additional libraries if needed

To deactivate and delete the conda environment:

conda deactivate
conda env remove --name conda-far

Installation testing

You can test that the package has been correctly installed using the command forestatrisk in a terminal:

forestatrisk

This should return a short description of the forestatrisk package and the version number:

# forestatrisk: modelling and forecasting deforestation in the tropics.
# https://ecology.ghislainv.fr/forestatrisk/
# forestatrisk version x.x.

You can also test the package following the Get started tutorial.

Main functionalities

Sample

Function .sample() sample observations points from a forest cover change map. The sample is balanced and stratified between deforested and non-deforested pixels. The function also retrieves information from explanatory variables for each sampled point. Sampling is done by block to allow computation on large study areas (e.g. country or continental scale) with a high spatial resolution (e.g. 30m).

Model

Function .model_binomial_iCAR() can be used to fit the deforestation model. A linear Binomial logistic regression model is used in this case. The model includes an intrinsic Conditional Autoregressive (iCAR) process to account for the spatial autocorrelation of the observations. Parameter inference is done in a hierarchical Bayesian framework. The function calls a Gibbs sampler with a Metropolis algorithm written in pure C code to reduce computation time.

Other models (such as a simple GLM or a Random Forest model) can also be used.

Predict and project

Function .predict() allows predicting the deforestation probability on the whole study area using the deforestation model fitted with .model_*() functions. The prediction is done by block to allow the computation on large study areas (e.g. country or continental scale) with a high spatial resolution (e.g. 30m).

Function .deforest() predicts the future forest cover map based on a raster of probability of deforestation (rescaled from 1 to 65535), which is obtained from function .predict(), and an area (in hectares) to be deforested.

Validate

A set of functions (eg. .cross_validation() or .map_accuracy()) is also provided to perform model and map validation.

Contributing

The forestatrisk Python package is Open Source and released under the GNU GPL version 3 license. Anybody who is interested can contribute to the package development following our Community guidelines. Every contributor must agree to follow the project's Code of conduct.

forestatrisk's Issues

Histogram of probability values

Would be faster to use GDAL functions, eg. gdalinfo -hist

Remove calls to os.system()

Try to remove call to os.system(). See for example os.makedirs and calls to gdal. Try using GDAL/OGR Python API (http://gdal.org/python/).

CalledProcessError after running country_forest_run

Hi @ghislainv could you please give me some help here? I'm totally stuck with this error. I've tried several alternatives and haven't been able to resolve it. Any tips will help me a lot. Thanks.

far.data.country_forest_run(
iso3, proj="EPSG:4326",
output_dir="data_raw",
keep_dir=True,
fcc_source="gfc", perc=90,
gdrive_remote_rclone="gdrive_zwarg",
gdrive_folder="GEE-forestatrisk-notebooks")

---

CalledProcessError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_15332/1526280857.py in
6 fcc_source="gfc", perc=90,
7 gdrive_remote_rclone="gdrive_zwarg",
----> 8 gdrive_folder="GEE-forestatrisk-notebooks")

~\miniconda3\envs\conda-far\lib\site-packages\forestatrisk\build_data\data.py in country_forest_run(iso3, proj, output_dir, keep_dir, fcc_source, perc, gdrive_remote_rclone, gdrive_folder)
170 # Check data availability
171 data_availability = ee_gfc.check(gdrive_remote_rclone,
--> 172 gdrive_folder, iso3)
173 # If not available, run GEE
174 if data_availability is False:

~\miniconda3\envs\conda-far\lib\site-packages\forestatrisk\build_data\ee_gfc.py in check(gdrive_remote_rclone, gdrive_folder, iso3)
125 cmd = ["rclone", "lsf", "--include", pattern, remote_path]
126 cmd = " ".join(cmd)
--> 127 out = subprocess.check_output(cmd, shell=True).decode("utf-8")
128 # Filename to find
129 fname = "forest_" + iso3

~\miniconda3\envs\conda-far\lib\subprocess.py in check_output(timeout, *popenargs, **kwargs)
409
410 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
--> 411 **kwargs).stdout
412
413

~\miniconda3\envs\conda-far\lib\subprocess.py in run(input, capture_output, timeout, check, *popenargs, **kwargs)
510 if check and retcode:
511 raise CalledProcessError(retcode, process.args,
--> 512 output=stdout, stderr=stderr)
513 return CompletedProcess(process.args, retcode, stdout, stderr)
514

CalledProcessError: Command 'rclone lsf --include 'forest_BRA*.tif' gdrive_zwarg:GEE-forestatrisk-notebooks' returned non-zero exit status 1.

[JOSS Review] Missing links in Community Guidelines

Hi @ghislainv,

I'm reviewing your software for publication in the JOSS. The journal requires instructions for collaborators and contributors detailed in a community guideline and a code of conduct.
The ´README.md´ file includes a section for these items but the links to these resources seem broken.

Cheers,

Juan

Function to predict deforestation without spatial random effects

We can either:

• build an hSDM_binomial function
• or use a linear model from module statsmodels

Problem with gadm3.6 access

gadm server for gadm3.6 seems to be down since few days (at least, I can't access it, from different locations and computers).
gadm4.1 seems to be ok though. Any reason not to update to this version ?

forestatrisk/forestatrisk/data/download/download_gadm.py

Line 38 in d1c813c

"shp/gadm36_" + iso3 + "_shp.zip"

Replace patsy with formulaic

patsy is a dependency of forestatrisk. patsy is not maintained anymore and has been replaced with formulaic.

patsy should then be replaced by formulaic as a dependency for the forestatrisk package.

More info for formulaic and the migration here:

• https://github.com/pydata/patsy
• https://github.com/matthewwardrop/formulaic
• https://matthewwardrop.github.io/formulaic/migration/

Forest cover input

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

Make docs

Make a web site with the doc for the package. See these tools:

• https://www.mkdocs.org/
• https://brendanhasz.github.io/2019/01/05/sphinx.html

Test failed on GitHub but not locally.

Tests failed on GitHub but not locally. Might be a random seed problem. To be checked...

Overview forest cover 2050

Due to the way pyramids are computed, ares with forest loss do not appear in red (cf. difference with probability of deforestation). It should be better to make new figures and png files correcting that.

Save parameter estimates

Create a function to save parameter estimates

Update `Articles` in documentation

Due to a change in API and introduction of geefcc dependency, articles need to be updated.

Do not recall Google Earth Engine if results already in Google Cloud Storage

Write a function to avoid recalling Google Earth Engine if results already in Google Cloud Storage bucket.

Update country list for Africa

Remove countries without humid tropical forests. E.g. Senegal, Gambia

far.predict_raster_binomial_iCAR

Dear @ghislainv

Could you please explain what is the forest_t3.tif that you have been used in the far.predict_raster_binomial_iCAR? I'm assuming that t1 in your example is forest in 2010, and t2 forest in 2020. So, what about t3?

Thank you,

All the best

Ricardo.

Add an option to build data for Brazilian states automatically

See this related issue #74.

Replace sklearn in dependencies

Replace witk scikit-learn or statsmodels to avoid having multiple dependencies and to avoid having to install these three packages for the tutorial.

Argument to keep/remove data_raw

In function dfp.data.country, add an argument to keep data_raw folder.

change gdal dependency to pygdal ?

forestatrisk/setup.py

Line 80 in 20d077f

install_requires=["earthengine-api", "gdal", "numpy", "matplotlib",

You define as a dependency gdal which is the worst possible neighbore in a python environment. What would you say if instead we use the saythised ersion of it: pygdal ?

I don't know what gdal is used for but if it's possible to get fully rid of it it's even better (it's embeded in the wheel of rasterio, rioxarray, geopandas)

Compute the forest data but HTTP error and doesn’t recognize BRA-AC--reverts to Caledonia extent--error: no driver to open shape file for Brazil!?

Rconfig points to all the right folders in terminal...

Compute gee forest data

far.data.country_forest_run(
iso3, proj="EPSG:4326",
output_dir="data_raw",
keep_dir=True,
fcc_source="jrc", perc=50,
gdrive_remote_rclone="gdrive_xxx",
gdrive_folder="GEE-forestatrisk-notebooks")

HTTPError Traceback (most recent call last)
/var/folders/hv/rrx0n8mj291dq3qrfcb9x9bw0000gn/T/ipykernel_68748/465982169.py in <cell line: 2>()
1 # Compute gee forest data
----> 2 far.data.country_forest_run(
3 iso3, proj="EPSG:4326",
4 output_dir="data_raw",
5 keep_dir=True,

~/opt/miniconda3/envs/FAR/lib/python3.10/site-packages/forestatrisk/build_data/data.py in country_forest_run(iso3, proj, output_dir, keep_dir, fcc_source, perc, gdrive_remote_rclone, gdrive_folder)
137 fname = output_dir + "/" + iso3 + "shp.zip"
138 url = "http://biogeo.ucdavis.edu/data/gadm3.6/shp/gadm36" + iso3 + "_shp.zip"
--> 139 urlretrieve(url, fname)
140
141 # Extract files from zip

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in urlretrieve(url, filename, reporthook, data)
239 url_type, path = _splittype(url)
240
--> 241 with contextlib.closing(urlopen(url, data)) as fp:
242 headers = fp.info()
243

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
214 else:
215 opener = _opener
--> 216 return opener.open(url, data, timeout)
217
218 def install_opener(opener):

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in open(self, fullurl, data, timeout)
523 for processor in self.process_response.get(protocol, []):
524 meth = getattr(processor, meth_name)
--> 525 response = meth(req, response)
526
527 return response

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in http_response(self, request, response)
632 # request was successfully received, understood, and accepted.
633 if not (200 <= code < 300):
--> 634 response = self.parent.error(
635 'http', request, response, code, msg, hdrs)
636

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in error(self, proto, *args)
555 http_err = 0
556 args = (dict, proto, meth_name) + args
--> 557 result = self._call_chain(*args)
558 if result:
559 return result

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
494 for handler in handlers:
495 func = getattr(handler, meth_name)
--> 496 result = func(*args)
497 if result is not None:
498 return result

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in http_error_302(self, req, fp, code, msg, headers)
747 fp.close()
748
--> 749 return self.parent.open(new, timeout=req.timeout)
750
751 http_error_301 = http_error_303 = http_error_307 = http_error_302

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in open(self, fullurl, data, timeout)
523 for processor in self.process_response.get(protocol, []):
524 meth = getattr(processor, meth_name)
--> 525 response = meth(req, response)
526
527 return response

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in http_response(self, request, response)
632 # request was successfully received, understood, and accepted.
633 if not (200 <= code < 300):
--> 634 response = self.parent.error(
635 'http', request, response, code, msg, hdrs)
636

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in error(self, proto, *args)
555 http_err = 0
556 args = (dict, proto, meth_name) + args
--> 557 result = self._call_chain(*args)
558 if result:
559 return result

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
494 for handler in handlers:
495 func = getattr(handler, meth_name)
--> 496 result = func(*args)
497 if result is not None:
498 return result

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in http_error_302(self, req, fp, code, msg, headers)
747 fp.close()
748
--> 749 return self.parent.open(new, timeout=req.timeout)
750
751 http_error_301 = http_error_303 = http_error_307 = http_error_302

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in open(self, fullurl, data, timeout)
523 for processor in self.process_response.get(protocol, []):
524 meth = getattr(processor, meth_name)
--> 525 response = meth(req, response)
526
527 return response

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in http_response(self, request, response)
632 # request was successfully received, understood, and accepted.
633 if not (200 <= code < 300):
--> 634 response = self.parent.error(
635 'http', request, response, code, msg, hdrs)
636

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in error(self, proto, *args)
561 if http_err:
562 args = (dict, 'default', 'http_error_default') + orig_args
--> 563 return self._call_chain(*args)
564
565 # XXX probably also want an abstract factory that knows when it makes

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
494 for handler in handlers:
495 func = getattr(handler, meth_name)
--> 496 result = func(*args)
497 if result is not None:
498 return result

~/opt/miniconda3/envs/FAR/lib/python3.10/urllib/request.py in http_error_default(self, req, fp, code, msg, hdrs)
641 class HTTPDefaultErrorHandler(BaseHandler):
642 def http_error_default(self, req, fp, code, msg, hdrs):
--> 643 raise HTTPError(req.full_url, code, msg, hdrs, fp)
644
645 class HTTPRedirectHandler(BaseHandler):

HTTPError: HTTP Error 404: Not Found

Annual forest loss

Create a function to compute the historical mean annual deforestation rate (in ha).

Installation testing

I did the installation via conda as described here and when testing the installation by running "forestatrisk" I get this import error:

no display found. Using non-interactive Agg backend
RuntimeError: module compiled against API version 0xe but this version of numpy is 0xd
Traceback (most recent call last):`
  File "c:\users\felip\miniconda3\envs\conda-far\lib\runpy.py", line 193, in _run_module_as_main
    "__main__", mod_spec)
  File "c:\users\felip\miniconda3\envs\conda-far\lib\runpy.py", line 85, in _run_code
    exec(code, run_globals)
  File "C:\Users\felip\miniconda3\envs\conda-far\Scripts\forestatrisk.exe\__main__.py", line 4, in <module>
  File "c:\users\felip\miniconda3\envs\conda-far\lib\site-packages\forestatrisk\__init__.py", line 29, in <module>
    from .model import cellneigh, cellneigh_ctry
  File "c:\users\felip\miniconda3\envs\conda-far\lib\site-packages\forestatrisk\model\__init__.py", line 13, in <module>
    from .model_binomial_iCAR import model_binomial_iCAR
  File "c:\users\felip\miniconda3\envs\conda-far\lib\site-packages\forestatrisk\model\model_binomial_iCAR.py", line 24, in <module>
    from .. import hbm
ImportError: numpy.core.multiarray failed to import

[JOSS Review] Automated test failed

Hi @ghislainv,

The automated tests (located in tests) failed. I'm appending the obtained output.

(conda-far) otus:forestatrisk juan$ py.test
========================================== test session starts ==========================================
platform darwin -- Python 3.7.7, pytest-6.2.2, py-1.10.0, pluggy-0.13.1
rootdir: /.../JOSS/forestatrisk
collected 17 items

test/test_get_started.py ......FF........                                                         [ 94%]
test/test_import.py .                                                                             [100%]

=============================================== FAILURES ================================================
_______________________________________ test_model_binomial_iCAR ________________________________________

gstart = {'adj': array([ 1, 11, 12,  0,  2, 11, 12, 13,  1,  3, 12, 13, 14,  2,  4, 13, 14,
       15,  3,  5, 14, 15, 16,  4, ...       5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 5, 5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 5,
       3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 3]), ...}

    def test_model_binomial_iCAR(gstart):

        p = np.array([0.34388896, 0.29002158, 0.51594223, 0.48436339,
                      0.60838453, 0.61257058, 0.55034979, 0.58819568,
                      0.51087469, 0.58819568, 0.64149789, 0.57400436,
                      0.59570952, 0.63212285, 0.566676, 0.62562204,
                      0.55379459, 0.15644965, 0.61284327, 0.36638686,
                      0.55439297, 0.57325744, 0.62562204, 0.17995823,
                      0.4930868, 0.54641479, 0.59782004, 0.48159526,
                      0.62882886, 0.59831051, 0.76245777, 0.74576097,
                      0.77356767, 0.73863295, 0.78188891, 0.75056545,
                      0.60775752, 0.64978574, 0.74654465, 0.77378323,
                      0.53994416, 0.75852715, 0.77754366, 0.60053684,
                      0.71543739, 0.74565542, 0.7555028, 0.44598923,
                      0.76401273, 0.75953027, 0.49027142, 0.69610182,
                      0.75679461, 0.78543649, 0.76863321, 0.6209473,
                      0.77653139, 0.76182804, 0.78169681, 0.58816002,
                      0.50453473, 0.77980428, 0.76084413, 0.73054832,
                      0.78289747, 0.71858934, 0.78362842, 0.74702923,
                      0.67357571, 0.78940242, 0.75358937, 0.66791346,
                      0.75602843, 0.42494845, 0.77653139, 0.60509306,
                      0.60846943, 0.76187008, 0.73278992, 0.72792572,
                      0.47661681, 0.59456417, 0.71894598, 0.6731302,
                      0.74964489, 0.77247818, 0.78289747, 0.74200682,
                      0.78940242, 0.78508877, 0.73153419, 0.65636031,
                      0.78607775, 0.59738545, 0.72596162, 0.78216462,
                      0.75078253, 0.77527468, 0.69907386, 0.71991522])

>       assert np.allclose(gstart["pred_icar"][0:100], p)
E       assert False
E        +  where False = <function allclose at 0x10435e8c0>(array([0.38174266, 0.34110367, 0.49796017, 0.4767064 , 0.56710096,\n       0.57024475, 0.5418506 , 0.55263309, 0.512898...77, 0.67407704, 0.77655815, 0.57789   , 0.76288142,\n       0.77055253, 0.73327539, 0.76372531, 0.67347268, 0.68633676]), array([0.34388896, 0.29002158, 0.51594223, 0.48436339, 0.60838453,\n       0.61257058, 0.55034979, 0.58819568, 0.510874...19, 0.65636031, 0.78607775, 0.59738545, 0.72596162,\n       0.78216462, 0.75078253, 0.77527468, 0.69907386, 0.71991522]))
E        +    where <function allclose at 0x10435e8c0> = np.allclose

/../forestatrisk/test/test_get_started.py:132: AssertionError
_______________________________________________ test_rho ________________________________________________

gstart = {'adj': array([ 1, 11, 12,  0,  2, 11, 12, 13,  1,  3, 12, 13, 14,  2,  4, 13, 14,
       15,  3,  5, 14, 15, 16,  4, ...       5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 5, 5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 5,
       3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 3]), ...}

    def test_rho(gstart):

        r = np.array([-3.72569484e-02, -1.16871478e-01, -1.82400711e-01,
                      2.13446770e-01, -6.44591325e-01, -9.89850864e-02,
                      1.10439030e-01, -2.31551563e-02, -3.30273946e-01,
                      -2.66995061e-01, -3.84426210e-01, 5.73572517e-02,
                      -5.73353804e-02, -3.12497338e-01, -8.37127591e-01,
                      7.62072575e-02, 3.86361945e-01, 1.26487021e-02,
                      -8.22069815e-02, -3.60656850e-01, -5.46586761e-01,
                      -4.17346094e-01, 1.05212875e+00, -4.32508096e-02,
                      -4.49589533e-01, -6.89872259e-01, -4.91230799e-01,
                      -3.84040358e-01, 5.67299746e-01, -2.10071117e-01,
                      -1.07456253e+00, -6.69339978e-01, -6.21974970e-01,
                      2.15020267e+00, -7.16437085e-02, -4.46424607e-01,
                      -2.17259138e-01, -3.30043032e-01, -2.59613996e-01,
                      2.68845283e-01, -3.78046974e-01, -5.18108829e-01,
                      -6.18235133e-01, -7.59652734e-01, 1.51771355e+00,
                      1.75357016e+00, -8.01814048e-02, 1.99270623e-01,
                      -1.75157345e-01, -6.10561635e-02, -1.26099802e-01,
                      -1.77864133e-01, -3.03381214e-01, -5.29892286e-01,
                      -5.47125418e-01, 1.30320979e+00, 2.37670385e+00,
                      4.97829325e-01, 8.88668246e-01, 3.92682659e-01,
                      -6.56913949e-03, -2.95774565e-01, -5.15489012e-01,
                      -6.01407176e-01, -5.67695385e-01, -6.48479745e-01,
                      1.47482553e+00, 1.45848019e+00, 4.05321503e-01,
                      1.06327906e+00, 4.37780456e-01, -1.12202021e-01,
                      -7.22139489e-01, -7.33312519e-01, -6.68442058e-01,
                      -7.76218335e-01, -8.02763852e-01, 1.41620727e+00,
                      1.56564133e+00, 1.24252305e+00, 9.07095194e-01,
                      4.38959947e-01, -2.95546782e-01, -4.92024764e-01,
                      -9.62965263e-01, -8.93107795e-01, -9.80673724e-01,
                      -9.94878624e-01, 1.41460696e+00, 1.38942057e+00,
                      1.97092977e+00, 1.06797639e+00, 4.36803818e-01,
                      2.15296806e-03, -6.14110567e-01, -7.76157636e-01,
                      -9.47693103e-01, -1.05424592e+00, -1.12226096e+00])

>       assert np.allclose(gstart["rho"], r)
E       assert False
E        +  where False = <function allclose at 0x10435e8c0>(array([ 0.0624157 , -0.15277958, -0.38287764,  0.23009508, -0.850125  ,\n       -0.24314522, -0.19327062, -0.58031306, ...36883 ,  1.16946349,  0.45945857, -0.15399801, -0.61728853,\n       -0.6482013 , -0.57797415, -0.38994641, -0.25915156]), array([-3.72569484e-02, -1.16871478e-01, -1.82400711e-01,  2.13446770e-01,\n       -6.44591325e-01, -9.89850864e-02,  1...3818e-01,  2.15296806e-03, -6.14110567e-01, -7.76157636e-01,\n       -9.47693103e-01, -1.05424592e+00, -1.12226096e+00]))
E        +    where <function allclose at 0x10435e8c0> = np.allclose

/../forestatrisk/test/test_get_started.py:171: AssertionError
=========================================== warnings summary ============================================
../../../miniconda3/envs/conda-far/lib/python3.7/site-packages/osgeo/__init__.py:8
  /../miniconda3/envs/conda-far/lib/python3.7/site-packages/osgeo/__init__.py:8: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses
    import imp

../../../miniconda3/envs/conda-far/lib/python3.7/site-packages/patsy/constraint.py:13
  /../miniconda3/envs/conda-far/lib/python3.7/site-packages/patsy/constraint.py:13: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3,and in 3.9 it will stop working
    from collections import Mapping

-- Docs: https://docs.pytest.org/en/stable/warnings.html
======================================== short test summary info ========================================
FAILED test/test_get_started.py::test_model_binomial_iCAR - assert False
FAILED test/test_get_started.py::test_rho - assert False
=============================== 2 failed, 15 passed, 2 warnings in 51.58s ===============================
:100.0%, mean accept. rates= beta:0.350, rho:0.443

Is it my side?
I don't have much idea of what could it be. However, the getting started tutorial worked 100% and the installation process showed no problems.

Cheers,

Juan

Confusion regarding Cell column in the sample.txt file and the data in forest_t3.tif

I was trying to follow the steps in https://ecology.ghislainv.fr/forestatrisk/notebooks/get_started.html#1.-Data. I am confused about the Cell column in the sample.txt file generated.

Also, what is the forest_t3.tif file used while prediction? Is it the list of pixels for which predictions needs to be done?

Remove `numpy.distutils` dependency

forestatrisk currently uses numpy.distutils which is now deprecated

Follow the migration advice to get rid of it, and make forestatrisk work with numpy>=1.23 and python>=3.12.

Organize package in submodules

With several init.py for plot, data, models, sample, predict and validation.

ReadAsArray and NoData

In function emissions.py, there is a problem when setting NoData to zero. The NoData value in the raster has not the exact same value as the NoData value in the array obtained with ReadAsArray. Maybe a type problem.

Sample number of observations proportionnally to forest area

Write a function to sample the number of points as a function of the forest area with a minimum of 20000 points.

Remove non-significant effect

Write a function for model selection and for removing variable with statistically insignificant effect.

Forest pixel with NA

When a forest pixel has NA as value for environmental variables, the spatial probability of deforestation cannot be computed for that forest pixel. It might happen for example for forest pixels close to the sea (mangroves) where the SRTM data is not provided. Write a function to empty holes with neighboring probability values.

Eg. Coast for Equatorial Guinea. Aspect and slope are not computed on the 90m along the coast. Also some forest pixels have NA values for aspect when slope=0.

Parallel computing for each country

Write a function for parallel computing. One CPU for each country.
Help can be found here: https://wiki.python.org/moin/ParallelProcessing

Carbon emissions

Create a function to compute carbon emissions associated to future deforestation.

drop gdal

As you can imagine I would like to set this as a module in SEPAL (and outside as well) but to do that I would need to drop the gdal requirements of the lib.

I can start the modifications in a PR but i would need some support to make sure I don't break anything.

Would you have some time for small support ?

Forest Cover Change Dataset

I have been trying to replicate this methodology to my study area, and I have had a few challenges, and a help from you will be great.

My challenge lies around the Land cover change datasets. From creating sample points in line 4 of your code right up till interpolate spatial random effects in line 9, the land cover change datasets used is "fcc23.tif". But in line 10 of your code in compute prediction the input_forest_raster makes use of "forest_t3.tif" which is different from the previous dataset you have been using, and it got me all confused.

So I wish to ask if you can please provide me with the methodology you used in creating the forest cover change dataset. Because what I did was created a supervised classification land use land cover map of my area for the year 2000 and 2021 then performed post-classification analysis to get the land cover change dataset for my area. But from the datasets you used, there are different types of forest cover change datasets.

Replace gsutil with google cloud storage Python API in ee_hansen.download()

Tor remove call to os.system(gsutil...).

Reduce computation type for gdalwrap -cutline

It might be more efficient to rasterize and then use gdal_calc to accelerate computation.

csize impact in iCAR model

Dear @ghislainv

How are you?

I want to bring here a question about the csize parameter in the functions far.sample and far.cellneigh, and their relation with the performance of the iCAR model.

To be more specific, when I run the iCAR model by setting the csize equal to 10 in the functons cited above, wich is the same value that you have been work with in the Get Started example, the model crashes taking a long time to respond. So, after several minutes I decided to interrompe the model and set the csize equal to 0.0001, and it solved the problem. But I'm still confuse what should be a acceptable csize according to the size of my area of interest (in this case ~ 1,5 millions hectares), and what is the impact of the csize in the model outputs.

Thank you very much @ghislainv

Using csize equal to 10 I got:

Compute number of 10 x 10 km spatial cells
... 1 cells (1 x 1)
Identify adjacent cells and compute number of neighbors

Changing the csize from 10 to 0.0001:

Compute number of 0.0001 x 0.0001 km spatial cells
... 144 cells (12 x 12)
Identify adjacent cells and compute number of neighbors