Deforestation tiles

Global Forest Change 2000-2018 v1.6 data

  • Download data.

  • First find the cumulative loss over all years.

  • Combine mosaic tiles.

  • Regrid.

  • Find the individual loss for each year.

import numpy as np
import gdal
import xarray as xr
from datetime import datetime
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from mpl_toolkits.basemap import Basemap
from matplotlib import colors
path = '/nfs/a68/earlacoa/deforestation/hansen_gfs_v1.6_lossyear/'
def regrid_mosaic_defor(res, year):
    """
    Description:
        The Hansen dataset provides the year of forest loss per grid cell between 2000 and 2018.
        For each year (e.g. 2005), convert the forest loss to binary by converting all forest loss up to that year
        (e.g. 2000-2005) to a 1 and all forest loss in future years to a 0 (e.g. 2006-2018).
        This is the cumulative forest loss.
        Then to obtain individual forest loss for each year after 2000, subtract the forest loss from the
        previous year (e.g. individual 2006 = cumulative 2006 - cumulative 2005).
    Args:
        res (float): desired resolution in degrees.
        year (int): year.
    Returns:
        hdat (2d array, float): regridded deforestation data.
        lat (1d array, float): latitude.
        lon (1d array, float): longitude.
    """
    # setup domain
    lon_extent_min = -180
    lon_extent_max = 180
    lat_extent_min = -60
    lat_extent_max = 80
    scale = int(res / abs(0.00025)) # fixed
    ydim = (int((lat_extent_max - lat_extent_min) / res))
    xdim1 = (int(10 / res)) # fixed
    xdim2 = (int((lon_extent_max - lon_extent_min) / res))
    vdat = np.empty((ydim, xdim1))
    hdat = np.empty((ydim, xdim2))
    fname = ''
    j = 0
    
    for nx in np.arange(-lon_extent_min, -lon_extent_max, -10):
        i = 0
        
        for ny in np.arange(lat_extent_max, lat_extent_min, -10):
            # define the filename based on lat and lon as uses north, south, east, and west in filenames
            if (nx > 0 or nx == 180) and ny >= 0:
                xstr = str(nx).zfill(3)
                ystr = str(ny).zfill(2)
                fname = 'Hansen_GFC-2018-v1.6_lossyear_' + ystr + 'N_' + xstr + 'W.tif'
            elif (nx > 0 or nx == 180) and ny < 0:
                xstr = str(nx).zfill(3)
                ystr = str(abs(ny)).zfill(2)
                fname = 'Hansen_GFC-2018-v1.6_lossyear_' + ystr + 'S_' + xstr + 'W.tif'
            elif nx <= 0 and ny >= 0:
                xstr = str(abs(nx)).zfill(3)
                ystr = str(ny).zfill(2)
                fname = 'Hansen_GFC-2018-v1.6_lossyear_' + ystr + 'N_' + xstr + 'E.tif' 
            elif nx <= 0 and ny < 0:
                ystr = str(abs(ny)).zfill(2)
                xstr = str(abs(nx)).zfill(3)
                fname = 'Hansen_GFC-2018-v1.6_lossyear_' + ystr + 'S_' + xstr + 'E.tif'
            print(fname)
            
            # read data
            ds = gdal.Open(path + fname)
            band = ds.GetRasterBand(1)
            loss_data = band.ReadAsArray()
            
            # create new lat and lon
            gt  = ds.GetGeoTransform()
            lon = np.linspace(gt[0], gt[0] + gt[1] * loss_data.shape[1], loss_data.shape[1])
            lat = np.linspace(gt[3], gt[3] + gt[5] * loss_data.shape[0], loss_data.shape[0])
            xx, yy = lon[::scale], lat[::scale] 
            
            # transform to binary for the year of interest
            loss_data[loss_data > (year - 2000)] = 0
            loss_data[(loss_data >= 1) & (loss_data <= (year - 2000))] = 1
            
            # equate coarse grid to high-res grid
            data_regrid = np.zeros((xdim1, xdim1))
            
            for ix in range(len(xx)):
                temp_lon_ix = np.where((lon == xx[ix]))[0]
                
                for iy in range(len(yy)):                   
                    temp_lat_iy = np.where((lat == yy[iy]))[0]
                    
                    # trim the data
                    data_trim = (loss_data[temp_lat_iy[0]:temp_lat_iy[0] + scale, :][:, temp_lon_ix[0]:temp_lon_ix[0] + scale])
                    
                    # find average across this trimmed data
                    mean_loss = np.nanmean(data_trim)
                    
                    # set this mean to the corresponding value on the coarser grid
                    data_regrid[iy, ix] = mean_loss
                    
            # add each tiles' regridded data to a new global data
            vdat[i * xdim1:(i + 1) * xdim1, :] = data_regrid
            i += 1
            
        hdat[:, j * xdim1:(j + 1) * xdim1] = vdat
        j += 1
        
    # create final global lat and lon
    lat = np.arange(lat_extent_max, lat_extent_min, -res)
    lon = np.arange(lon_extent_min, lon_extent_max, res)
    
    return hdat, lat, lon
def plot(defor_xx, defor_yy, defor_array, year, values_max, label):
    """
    Description:
        Create a contoured plot of the global deforestation.
        Either cumulatively or individually per year.
    Args:
        defor_xx (2d array, float): longitude.
        defor_yy (2d array, float): latitude.
        defor_array (2d array, float): deforestation.
        year (int): year.
        values_max (float): maximum value for the colour bar.
        label (str): cumulative or individual.
    Returns:
        Plot displayed and saved to file.
    """
    fig = plt.figure(1, figsize=(14, 7))
    gs = gridspec.GridSpec(1, 1)
    ax = fig.add_subplot(gs[0])
    basemap = Basemap(
        llcrnrlon=-180, llcrnrlat=-90, urcrnrlon=180, urcrnrlat=90,
        projection='cyl', resolution='l'
    )
    basemap.drawcountries(linewidth=0.2)
    basemap.drawcoastlines(linewidth=0.2)
    basemap.fillcontinents(color='lightgrey', zorder=0)
    basemap.drawparallels(np.arange(-80., 81., 10.), labels=[1, 0, 0, 0], fontsize=14, linewidth=0)
    basemap.drawmeridians(np.arange(-180., 181., 30.), labels=[0, 0, 0, 1], fontsize=14, linewidth=0)
    cmap = 'viridis'
    norm = colors.Normalize(vmin=0, vmax=values_max)
    im = basemap.contourf(
        defor_xx, defor_yy, defor_array, 
        np.linspace(0, values_max, 11), cmap=cmap, norm=norm
    )
    sm = plt.cm.ScalarMappable(norm=norm, cmap=im.cmap)
    sm.set_array([])
    cb = basemap.colorbar(sm, "right", size="5%", pad='2%', norm=norm, cmap=cmap, ticks=im.levels)
    cb.set_label(
        'Hansen GFC 2018 v1.6, loss year, 0.25 degrees,\n fractional deforestation (0-1)', 
        size=14
    )
    cb.ax.tick_params(labelsize=14)
    plt.title(str(year), size=16)
    plt.tick_params(axis='both', which='major', labelsize=14)
    plt.tick_params(axis='both', which='minor', labelsize=14)
    plt.savefig(
        path + 'Hansen_GFC-2018-v1.6_lossyear-' + label + '_global_0.25deg_' + str(year) + '.png', 
        dpi=200, alpha=True, bbox_inches='tight'
    )
    plt.show()
years = np.linspace(2000, 2018, 19, dtype=int)
res = 0.25
defor_cumulative = {}
defor_individual = {}

for year in years:
    # run regridding and mosaicing function
    output, lat, lon = regrid_mosaic_defor(res, year)

    # create dataarray
    defor = xr.DataArray(
        data = output, 
        coords = [lat, lon],
        dims = ['lat', 'lon']
    )
    defor.name = 'deforestation'
    defor = defor.assign_coords({'time': datetime.strptime(str(year)[-2:], '%y')})
    defor = defor.expand_dims('time')
    defor.to_netcdf(path + 'Hansen_GFC-2018-v1.6_lossyear-cumulative_global_0.25deg_' + str(year) + '.nc')

    defor_array = defor.isel(time=0).values
    defor_xx, defor_yy = np.meshgrid(defor.lon.values, defor.lat.values)
    
    # plot cumulative
    plot(defor_xx, defor_yy, defor_array, year, 1, 'cumulative')
    
    # calculate individual yearly loss
    defor_cumulative.update({year: xr.open_dataset(
        path + 'Hansen_GFC-2018-v1.6_lossyear-cumulative_global_' + str(res) + 'deg_' + str(year) + '.nc'
    )['deforestation']})
    if year == 2000:
        pass
    else:
        defor_individual.update({
            year: defor_cumulative[year] - defor_cumulative[year - 1]
        })
    
    defor_individual[year].to_netcdf(
        path + 'Hansen_GFC-2018-v1.6_lossyear-individual_global_' + str(res) + 'deg_' + str(year) + '.nc'
    )
    
    # plot individual
    defor_array = defor_individual[year].isel(time=0).values
    plot(defor_xx, defor_yy, defor_array, year, 0.5, 'individual')