From 902476bc8b14e598e3c637f95ce3619d945f1d54 Mon Sep 17 00:00:00 2001 From: Maria Voigt <maria.voigt@idiv.de> Date: Thu, 10 Aug 2017 11:27:35 +0200 Subject: [PATCH] renaming the preparation of resource grid and future grid --- .../grid_mapping_for_resource_use.py | 214 ++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 src/validation/grid_mapping_for_resource_use.py diff --git a/src/validation/grid_mapping_for_resource_use.py b/src/validation/grid_mapping_for_resource_use.py new file mode 100644 index 0000000..bbc505f --- /dev/null +++ b/src/validation/grid_mapping_for_resource_use.py @@ -0,0 +1,214 @@ +# -*- coding: utf-8 -*- +""" +Spyder Editor + +# script to prepare the grid with the resource use information +# in the output file we will have a +# a 1 at the first position if pixel has to be discarded +# a 2 at the first position if the pixel is considered +# a 1 at the second position for IDN +# a 2 at the second position for MYS +# a 0 - 5 at the 3rd position depending on the resource use category +# the grid_id from the 4rth to the 10th position (6 positions) +# +# now assigning number for each resource use +# 0 - absence +# 1 - plantation +# 2 - deforestation +# 3 - landcover change +# 4 - logging +# 5 - other +""" + + +import numpy as np +import MacPyver as mp +import os + + +""" +# absence +os.system("gdal_rasterize \ +-l absence_shape_rep \ +-burn 1 \ +-ot Float32 \ +-te -1751798.359 1267454.241 -629798.359 2560454.241 -ts 1122 1293 \ +/homes/mv39zilo/work/Borneo/data/response/cleaned_data/absence_shape/absence_shape_rep.shp \ +/homes/mv39zilo/work/Borneo/data/response/cleaned_data/absence_shape/absence_shape_rep.tif") + +os.system("gdal_rasterize \ +-l absence_shape_rep \ +-burn 1 \ +-ot Float32 \ + -ts 1122 1293 \ +/homes/mv39zilo/work/Borneo/data/response/cleaned_data/absence_shape/absence_shape_rep.shp \ +/homes/mv39zilo/work/Borneo/data/response/cleaned_data/absence_shape/absence_shape_rep.tif") + +os.system("gdalwarp \ +-t_srs '+proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' \ +-ot Float64 -r near \ +-te -1751798.359 1267454.241 -629798.359 2560454.241 \ +-ts 1122 1293 \ +-overwrite \ +/homes/mv39zilo/work/Borneo/data/response/cleaned_data/absence_shape/absence_shape_rep.tif \ +/homes/mv39zilo/work/Borneo/data/response/cleaned_data/absence_shape/absence_shape_repro_res.tif" +) + +os.system("gdalwarp \ +-t_srs '+proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' \ +-ot Float64 -r near \ +-te -1751798.359 1267454.241 -629798.359 2560454.241 \ +-ts 1122 1293 \ +-overwrite \ +/homes/mv39zilo/work/Borneo/data/future/grid.tif \ +/homes/mv39zilo/work/Borneo/data/future/grid_repro_res.tif" +) + + + +""" + +absence_path = "/homes/mv39zilo/work/Borneo/data/response/cleaned_data/absence_shape/absence_shape_repro_res.tif" +absence = mp.tiff.read_tif(absence_path, 1) + +# input grid with grid_ids as values + +grid = mp.tiff.read_tif("/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/future/grid_with_id_repro_res.tif", 1) +np.unique(grid) + + + + + + +gaveau_defor_path = "/homes/mv39zilo/work/Borneo/data/predictors/raw_data/Remaining_forest_in_2015_and_deforestation_1973-2015/REGIONBorneo_FCDefDeg_1973to2015_CIFOR_repro_res.tif" +gaveau_defor_raw = mp.tiff.read_tif(gaveau_defor_path, 1) +gaveau_deforested = np.where((gaveau_defor_raw > 6 ) & + (gaveau_defor_raw < 10), 1, 0) +gaveau_logged = np.where(gaveau_defor_raw == 2, 1, 0) + +#plantation +plantation_path = "/homes/mv39zilo/work/Borneo/data/predictors/raw_data/Land_cover_trajectory_before_oil-palm_and_pulpwood_establishment/raster_plantation.tif" +plantations = mp.tiff.read_tif(plantation_path, 1) + +#non habitat +crisp_2000_path = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use//crisp_2000_repro_res.tif" +crisp_2000 = mp.tiff.read_tif(crisp_2000_path, 1) +crisp_2010_path = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use/crisp_2010_repro_res.tif" +crisp_2010 = mp.tiff.read_tif(crisp_2010_path, 1) +crisp_2015_path = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use/crisp_2015_repro_res.tif" +crisp_2015 = mp.tiff.read_tif(crisp_2015_path, 1) + +cover_change = np.where(((crisp_2000 > 2 ) & (crisp_2000 < 7)) & ((crisp_2010 > 6 ) |(crisp_2010 == 1 ) | (crisp_2015 > 6) | (crisp_2015 == 1 ) ), 1, 0) + + + + +# pixel is out +grid = np.where(absence == 1, 1000000000 + grid, grid) +# pixel is in +grid = np.where(absence == 0, 2000000000 + grid, grid) + + + + # now assigning number for each resource use +# 0 - absence +# 1 - plantation +# 2 - deforestation +# 3 - landcover change +# 4 - logging +# 5 - other + + +# we can save some computing time by only considering +# the grid and not everything, because abundance will only be +# in rid +# problem with overlap, need to do like with the table, otherwise I add up + +grid = np.where((plantations == 1), 10000000 + grid , grid) +grid = np.where((gaveau_deforested == 1), 20000000 + grid , grid) +grid = np.where((cover_change == 1), 30000000 + grid , grid) +grid = np.where((gaveau_logged == 1), 40000000 + grid , grid) +grid = np.where((plantations == 0) & + (gaveau_deforested == 0) & + (cover_change == 0) & + (gaveau_logged == 0), 50000000 + grid , grid) + + + + +np.unique(grid) + +mp.tiff.write_tif(file_with_srid = grid_layer_path, + full_output_name = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use/resource_use_grid.tif", + data = grid, + dtype = 4) + +# use this in R script "prepare_boot_grid.R" + + +""" +# test +indir = "/homes/mv39zilo/work/Borneo/analysis/model_prep_and_running/results/resource_use" +abundance_layer_path = indir + "/abundance_pred_2015_no_absence.tif" +# read the data +abundance_data = mp.tiff.read_tif(abundance_layer_path, 1) + +test = np.where((grid == 0) & + (abundance_data > 0), abundance_data, 0) +np.sum(test) +mp.tiff.write_tif(file_with_srid = grid_layer_path, + full_output_name = "/homes/mv39zilo/work/Borneo/data/bootstrap/test.tif", + data = test, + dtype = 4) +""" +# rasterize borneo with info on country or province + +""" +os.system("gdal_rasterize \ +-a ID_1 \ +-l Borneo \ +-ot Float32 \ +-te -1751798.359 1267454.241 -629798.359 2560454.241 -ts 1122 1293 \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo.shp \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo_province.tif") + +os.system("gdal_rasterize \ +-a ID_1 \ +-l Borneo \ +-ot Float32 \ + -ts 1122 1293 \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo.shp \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo_province.tif") + +os.system("gdalwarp \ +-r near -t_srs '+proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' \ +-te -1751798.359 1267454.241 -629798.359 2560454.241 -ts 1122 1293 \ +-overwrite \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo_province.tif \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo__province_repro_res.tif") + +os.system("gdal_rasterize \ +-a ID_0 \ +-l Borneo \ +-ot Float32 \ +-te -1751798.359 1267454.241 -629798.359 2560454.241 -ts 1122 1293 \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo.shp \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo_country.tif") + +os.system("gdal_rasterize \ +-a ID_0 \ +-l Borneo \ +-ot Float32 \ + -ts 1122 1293 \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo.shp \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo_country.tif") + +os.system("gdalwarp \ +-r near -t_srs '+proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' \ +-te -1751798.359 1267454.241 -629798.359 2560454.241 -ts 1122 1293 \ +-overwrite \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo_country.tif \ +/homes/mv39zilo/work/Borneo/data/auxiliary_additional_data/Borneo_shape/cleaned_data/Borneo_country_repro_res.tif") + +""" -- GitLab