vulnerability_map.py

import os
import numpy as np
from osgeo import gdal
from PyQt5.QtCore import QObject, pyqtSignal
import shutil

# GDAL exceptions
gdal.UseExceptions()

class VulnerabilityMap(QObject):
    progress_updated = pyqtSignal(int)
    def __init__(self):
        super(VulnerabilityMap, self).__init__()
        self.data_folder = None
        self.initial_directory = None

    def set_working_directory(self, directory):
        '''
        Set up the working directory
        :param directory: your local directory with all dat files
        '''
        self.progress_updated.emit(0)
        self.data_folder = directory
        os.chdir(self.data_folder)

    def image_to_array(self,image):
        # Set up a GDAL dataset
        in_ds = gdal.Open(image)
        # Set up a GDAL band
        in_band = in_ds.GetRasterBand(1)
        # Create Numpy Array1
        arr = in_band.ReadAsArray()
        return arr

    def nrt_calculation(self, in_fn, deforestation_hrp, mask):
        '''
        NRT calculation
        :param in_fn: map of distance from the forest eddge in CAL
        :param deforestation_hrp:deforestation binary map in HRP
        :param mask: mask of the non-excluded jurisdiction (binary map)
        :return: NRT: Negligible Risk Threshold
        '''
        # Convert image to NumPy array
        self.progress_updated.emit(10)
        distance_arr_cal = self.image_to_array(in_fn)
        self.progress_updated.emit(30)
        deforestation_hrp_arr = self.image_to_array(deforestation_hrp)
        self.progress_updated.emit(40)
        mask_arr = self.image_to_array(mask)
        self.progress_updated.emit(50)

        # Mask the distance arr within deforstation pixel and study area
        distance_arr_masked=distance_arr_cal*mask_arr*deforestation_hrp_arr
        self.progress_updated.emit(60)

        ## Calculate the histogram
        # Flatten the distance_arr_masked and expect 0 for np.histogram function
        # The np.histogram is computed over the flattened array
        distance_arr_masked_1d = distance_arr_masked.flatten()
        self.progress_updated.emit(80)
        distance_arr_masked_1d = distance_arr_masked_1d[distance_arr_masked_1d != 0]

        ## Calculate the histogram
        # Set up bin width as spatial resolution
        in_ds = gdal.Open(in_fn)
        P = in_ds.GetGeoTransform()[1]
        bin_width = P
        # Calculate the histogram
        hist, bin_edges = np.histogram(distance_arr_masked_1d, bins=np.arange(distance_arr_masked_1d.min(),
                                                                              distance_arr_masked_1d.max() + bin_width,
                                                                               bin_width))
        # Calculate the cumulative proportion
        # Normalize the histogram to get probability
        hist_normalized = hist / np.sum(hist)

        # Compute cumulative distribution
        cumulative_prop = np.cumsum(hist_normalized)

        # # Find the index cumulative proportion >= 0.995
        index_995 = np.argmax(cumulative_prop >= 0.995)

        # Get the bin edges for the NRT bin
        nrt_bin_start = bin_edges[index_995]
        nrt_bin_end = bin_edges[index_995 + 1]
        self.progress_updated.emit(90)

        # Calculate the average of the NRT bin
        NRT = int((nrt_bin_start + nrt_bin_end) / 2)
        self.progress_updated.emit(100)
        return NRT

    def geometric_classification(self, in_fn, NRT, n_classes, mask):
        '''
        geometric classification
        :param in_fn: map of distance from the forest eddge
        :param NRT:Negligible Risk Threshold
        :param n_classes:number of classes
        :param mask:mask of the non-excluded jurisdiction (binary map)
        :return: mask_arr: result array with mask larger than NRT
        '''
        # Convert in_fn to NumPy array
        # Set up a GDAL dataset
        in_ds = gdal.Open(in_fn)
        # Set up a GDAL band
        in_band = in_ds.GetRasterBand(1)
        # Create Numpy Array
        arr = in_band.ReadAsArray()

        # The lower limit of the highest class = spatial resolution (the minimum distance possible without being in non-forest)
        LL = int(in_ds.GetGeoTransform()[1])

        self.progress_updated.emit(10)
        # The upper limit of the lowest class = the Negligible Risk Threshold
        UL = NRT = int(NRT)
        n_classes = int(n_classes)

        # Calculate common ratio(r)=(LLmax/LLmin)^1/n_classes
        r = np.power(LL / UL, 1/n_classes)

        # Create 2D class_array for the areas within the NRT
        class_array = np.array([[i, i + 1] for i in range(n_classes)])

        # Calculate UL and LL value for the areas within the NRT
        x= np.power(r, class_array)
        risk_class=np.multiply(UL,x)
        risk_class[n_classes-1][1] = LL
        risk_class[0][0] = NRT

        self.progress_updated.emit(20)

        # Multiple Distance to Non-Forest to
        mask_arr0 = self.image_to_array(mask)
        mask_arr=arr * mask_arr0
        # Create mask: areas beyond the NRT, assign class 1
        mask_arr[mask_arr >= NRT] = 1

        self.progress_updated.emit(30)
        for i in range(n_classes):
            lower = risk_class[i][0]
            upper = risk_class[i][1]
            mask_arr[(lower > mask_arr) & (mask_arr >= upper)] = i + 2
            if i % 5 == 0:
                progress = 40 + (10 * (i // 5))
                self.progress_updated.emit(progress)

        return mask_arr

    def geometric_classification_alternative(self, in_fn, n_classes, mask, fmask):
        '''
        geometric classification for alternative vulnerability map
        :param in_fn: Empirical vulnerability map [0.0,1.0] range
        :param n_classes:number of classes
        :param mask: mask of the non-excluded jurisdiction (binary map)
        :param fmask: mask of the forest areas (binary map)
        :return: mask_arr: result array with mask larger than NRT
        '''
        # Convert in_fn to NumPy array
        # Set up a GDAL dataset
        in_ds = gdal.Open(in_fn)
        # Set up a GDAL band
        in_band = in_ds.GetRasterBand(1)
        # Create Numpy Array
        arr = in_band.ReadAsArray()

        self.progress_updated.emit(10)

        max_value = in_band.GetMaximum()
        if max_value is None:
            # Calculate max_value for raster images without metadata
            max_value = np.max(arr)

        # Rescaled empirical vulnerability map to a [1.0–2.0] range
        arr_rescale = 1+arr*1/max_value

        # The lower limit of the highest class = 1
        LL = int(1)

        # The upper limit of the lowest class = 2
        UL = int(2)
        n_classes = int(n_classes)

        # Calculate common ratio(r)=(LLmax/LLmin)^1/n_classes
        r = np.power(UL / LL, 1 / n_classes)

        # Create 2D class_array
        class_array = np.array([[i, i + 1] for i in range(n_classes-1, -1, -1)])
        x = np.power(r, class_array)
        risk_class = LL+(UL-(np.multiply(LL, x)))

        # Explicitly set llmin
        risk_class[0][1] = LL

        self.progress_updated.emit(20)

        # Mask jurisdiction and forest area
        # Create jurisdiction mask array
        in_ds1 = gdal.Open(mask)
        in_band1 = in_ds1.GetRasterBand(1)
        mask_arr = in_band1.ReadAsArray()
        # Create forest area array
        in_ds2 = gdal.Open(fmask)
        in_band2 = in_ds2.GetRasterBand(1)
        fmask_arr = in_band2.ReadAsArray()
        # Array multiple mask and fmask
        mask_arr = arr_rescale*mask_arr*fmask_arr

        self.progress_updated.emit(30)

        for i in range(n_classes):
            upper = risk_class[i][0]
            lower = risk_class[i][1]
            mask_arr[(upper > mask_arr) & (mask_arr >= lower)] = i + 1
            if i % 5 == 0:
                progress = 40 + (10 * (i // 5))
                self.progress_updated.emit(progress)

        return mask_arr

    def array_to_image(self, in_fn, out_fn, data, data_type, nodata=None):
        '''
         Create image from array
        :param in_fn: datasource to copy projection and geotransform from
        :param out_fn: path to the file to create
        :param data: NumPy array containing data to write
        :param data_type: output data type
        :param nodata: optional NoData value
        :return:
        '''
        in_ds = gdal.Open(in_fn)
        output_format = out_fn.split('.')[-1].upper()
        if (output_format == 'TIF'):
            output_format = 'GTIFF'
        elif (output_format == 'RST'):
            output_format = 'rst'
        driver = gdal.GetDriverByName(output_format)
        out_ds = driver.Create(out_fn, in_ds.RasterXSize, in_ds.RasterYSize, 1, data_type, options=["BigTIFF=YES"])
        out_ds.SetProjection(in_ds.GetProjection().encode('utf-8', 'backslashreplace').decode('utf-8'))
        out_ds.SetGeoTransform(in_ds.GetGeoTransform())
        out_band = out_ds.GetRasterBand(1)
        if nodata is not None:
            out_band.SetNoDataValue(nodata)
        out_band.WriteArray(data)
        return

    def replace_ref_system(self, in_fn, out_fn):
        '''
         RST raster format: correct reference system name in rdc file
         :param in_fn: datasource to copy correct projection name
         :param out_fn: rst raster file
        '''
        if out_fn.split('.')[-1] == 'rst':
            if in_fn.split('.')[-1] == 'rst':
                read_file_name, _ = os.path.splitext(in_fn)
                write_file_name, _ = os.path.splitext(out_fn)
                temp_file_path = 'rdc_temp.rdc'

                with open(read_file_name + '.rdc', 'r') as read_file:
                    for line in read_file:
                        if line.startswith("ref. system :"):
                            correct_name = line
                            break

                if correct_name:
                    with open(write_file_name + '.rdc', 'r') as read_file, open(temp_file_path, 'w') as write_file:
                        for line in read_file:
                            if line.startswith("ref. system :"):
                                write_file.write(correct_name)
                            else:
                                write_file.write(line)

                    # Move the temp file to replace the original
                    shutil.move(temp_file_path, write_file_name + '.rdc')
                    self.progress_updated.emit(100)

            elif in_fn.split('.')[-1] == 'tif':
                # Read projection information from the .tif file using GDAL
                dataset = gdal.Open(in_fn)
                projection = dataset.GetProjection()
                dataset = None

                # Extract the reference system name from the wkt projection
                ref_system_name = projection.split('PROJCS["')[1].split('"')[0]

                write_file_name, _ = os.path.splitext(out_fn)
                temp_file_path = 'rdc_temp.rdc'

                with open(write_file_name + '.rdc', 'r') as read_file, open(temp_file_path, 'w') as write_file:
                    for line in read_file:
                        if line.startswith("ref. system :"):
                            write_file.write(f"ref. system : {ref_system_name}\n")
                        else:
                            write_file.write(line)

                shutil.move(temp_file_path, write_file_name + '.rdc')
                self.progress_updated.emit(100)
            else:
                self.progress_updated.emit(100)