__title__ = "Code for resampling Geoscape building and surface cover data"
__version__ = "2022-28-01"
__author__ = "<NAME>"
__email__ = "<EMAIL>"

'''
Associated with the manuscript: A transformation in city-descriptive input data for urban climate models
Manuscript authors:  <NAME>, <NAME>,  <NAME>,  <NAME>, <NAME>

Copyright 2022 <NAME>

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

Instructions: 
1. Get sample data from Geoscape: https://geoscape.com.au/get-sample/
2. Install necessary python dependancies
3. Update project and data paths as required
4. Run script
'''

import os
import numpy as np
import pandas as pd
import xarray as xr
import rioxarray as rxr
from rioxarray import merge

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

import shapely
import geopandas as gpd
import cartopy.geodesic as cgeo

################################################################################
# user input

domain = 'sample'
output_v = 'v0.97'

projpath = f'.'                   # update path to working folder
datapath = f'{projpath}/data'     # update path to geoscape data
outpath = f'{projpath}/outputs'   # netcdf/tif outputs
plotpath = f'{projpath}/figures'  # figures output

res_list = [30,100,300]         # list of resolutions to process in m

split_by_grid = False   # if true, splits buildings on grid edges (long processing)
plot_comp = True        # if true, plots comparison between different resolutions (for sample)

missing_float = -999.

################################################################################

def main():

    for res in res_list:
        print(f'processing {res} m')

        main_resample_cover(res)
        main_calculate_morphology(res)

    if plot_comp:
        main_plot_sample()

    return

def main_resample_cover(res):
    '''main function for resampling geoscape surface cover data to a lower resolution and recategorise'''

    print(f'opening {domain} geotif')

    # read tif
    tif_list = []
    for tif_fpath in tif_fpaths:
        tif_list.append(rxr.open_rasterio(tif_fpath))
    print('merging tifs')
    rorig =  merge.merge_arrays(tif_list)

    # get band 1
    orig_ds = rorig.drop(['band']).squeeze()
    # remove fill value
    orig_ds = orig_ds.where(orig_ds!=orig_ds._FillValue)

    ################################################################################

    print('creating downscaled geoscape dataset')
    geoscape = create_dataset_from_geoscape(orig_ds,res)

    ################################################################################

    print('saving geoscape dataset')
    fname = f'{outpath}/{domain.capitalize()}_geoscape_surface_cover_{res}mx{res}m_{output_v}' 

    encoding = {var:{'zlib':True} for var in geoscape.data_vars}
    geoscape.to_netcdf(f'{fname}.nc',format='NETCDF4', encoding=encoding)

    geoscape = xr.open_dataset(f'{fname}.nc')
    # geoscape.rio.write_crs("epsg:4326", inplace=True)
    # geoscape.rio.to_raster(f'{fname}.tiff', compress='LZW')

    ################################################################################

    derived = calc_derived_fractions(geoscape,res)
    derived = set_attributes(derived,res)

    print('saving derived surface cover')

    fname = f'{outpath}/{domain.capitalize()}_derived_surface_cover_{res}mx{res}m_{output_v}'

    compression = {'zlib':True}
    encoding = {var:compression for var in derived.data_vars}

    derived.to_netcdf(f'{fname}.nc',format='NETCDF4', encoding=encoding)
    # derived.rio.write_crs("epsg:4326", inplace=True)
    # derived.rio.to_raster(f'{fname}.tiff', compress='LZW')

    ################################################################################

    assert_tests(geoscape,derived)

    return

def calc_fraction(key,ds,res,long_name=None):
    '''function to take categorical band and calculate its fraction at lower resolution'''

    print(f'calculating {key} fraction at {res} m')

    # get resolution factors
    orig_res = ds.rio.resolution()[0]
    sfactor = int(res/orig_res)

    da = ds.where(ds == key)

    # get all key values, convert to 1
    tmp = da.copy()
    tmp.values = np.where(da.values==key,1,np.nan)
    # calculate sums over coarsened grid
    tmp_sum = tmp.coarsen(x=sfactor,y=sfactor,boundary='pad').sum()

    # get all non-nan values, convert to 1
    tmp = da.copy()
    tmp.values = np.where(np.isnan(ds).values,ds.values,1)
    # calculate sums over coarsened grid
    tmp_all = tmp.coarsen(x=sfactor,y=sfactor,boundary='pad').sum()

    result = tmp_sum/tmp_all
    result = result.rename(str(key))

    result.attrs['long_name'] = long_name
    result.encoding['_FillValue'] =  missing_float

    return result

def create_dataset_from_geoscape(orig_ds,res):

    '''create netcdf attributes from geoscape surface cover data

    # Table 2: SURFACE_COVER_2M_TYPE_AUT Codes from 'Surface Cover 1.6 Product Description.pdf'

    2 Bare Earth                Includes sand dunes, desert, rock outcrops, bare soil other than bare agricultural land, and sparsely vegetated areas of grass and shrub. Non-vegetated strip mines and quarries except where covered by development or water.
    3 Road and Path             Roads and parking lots covered in a man-made material excluding hard packed dirt trails.
    4 Grass                     Grass and herbaceous areas. The category may include herbaceous wetlands if images are collected during dry season or periods of drought.
    5 Trees                     All trees including deciduous and evergreen woody vegetation.
    6 Unspecified Vegetation    Any other vegetative material not included within the Grass or Tree class. This may include, but is not limited to, shrub, sc rub, agriculture, and aquatic plants.
    7 Built-up Areas            Any areas of man-made environments and infrastructure excluding road and paths and buildings.
    8 Water                     Depending on the resolution quality of the imagery used, natural water will include streams, canals, ponds, lakes, reservoirs, estuaries and bays.
    9 Buildings                 Where the majority of a pixel intersects a Building, vector building polygon representation. 
    10 Cloud                    The area covered with cloud on Date of collection.
    11 Shadow                   The area covered with shadow on Date/time of collection.
    12 Swimming Pool            An area identified as a swimming pool.
    '''
    
    print('creating dataset')

    ds = xr.Dataset({

        'bare_earth': calc_fraction(2,orig_ds,res=res,
            long_name= 'Bare earth: includes bare soil, sand, desert, rock outcrops'),

        'road_path': calc_fraction(3,orig_ds,res=res,
            long_name = 'Road and path: includes roads, parking lots and hard packed dirt trails'),

        'grass': calc_fraction(4,orig_ds,res=res,
            long_name = 'Grass and herbaceous areas'),

        'trees': calc_fraction(5,orig_ds,res=res,
            long_name = 'All trees including deciduous and evergreen'),

        'other_veg': calc_fraction(6,orig_ds,res=res,
            long_name = 'Other vegetation including shrub, scrub, agriculture and aquatic plants'),

        'built_area': calc_fraction(7,orig_ds,res=res,
            long_name = 'Human-made environments and infrastructure excluding road, paths and buildings'),

        'water': calc_fraction(8,orig_ds,res=res,
            long_name = 'Open water including streams, canals, lakes, reservoirs and bays'),

        'buildings': calc_fraction(9,orig_ds,res=res,
            long_name = 'Buildings'),

        'cloud': calc_fraction(10,orig_ds,res=res,
            long_name = 'Area covered with cloud at time of collection'),

        'shadow': calc_fraction(11,orig_ds,res=res,
            long_name = 'Area covered with shadow at time of collection'),

        'swimming_pool': calc_fraction(12,orig_ds,res=res,
            long_name = 'Swimming pool'),
        })

    ds.attrs['title'] = f'{domain.capitalize()} surface cover at {res}m resolution derived from 2m resolution GeoScape product'
    ds.attrs['version'] = output_v
    ds.attrs['authors'] = '<NAME> <<EMAIL>>'
    ds.attrs['institution'] = 'ARC Centre of Excellence for Climate Extremes, UNSW Sydney, Australia'
    ds.attrs['source'] = 'Geoscape Surface Cover v1.6. PSMA Australia, 2020. https://geoscape.com.au/'
    ds.attrs['comment'] = 'The contents of this file are restricted, not to be used or distributed without permission from the authors.'

    geoscape = ds.rio.reproject(dst_crs=4326)
    geoscape = geoscape.rename({'y':'latitude','x':'longitude'}).drop('spatial_ref')

    return geoscape

def calc_derived_fractions(ds,res):

    '''create derived classes and attributes (normalising for cloud and shadow fractions)'''

    derived = xr.Dataset()

    derived['building_fraction'] = (ds['buildings'])/(1. - ds['cloud'] - ds['shadow'])
    derived['tree_fraction'] = (ds['trees'])/(1. - ds['cloud'] - ds['shadow'])
    derived['lowveg_fraction'] = (ds['grass'] + ds['other_veg'])/(1. - ds['cloud'] - ds['shadow'])
    derived['bareearth_fraction'] = (ds['bare_earth'])/(1. - ds['cloud'] - ds['shadow'])
    derived['roadpath_fraction'] = (ds['road_path'] + ds['built_area'])/(1. - ds['cloud'] - ds['shadow'])
    derived['water_fraction'] = (ds['water'] + ds['swimming_pool'])/(1. - ds['cloud'] - ds['shadow'])
    derived['total_built'] = (ds['road_path'] + ds['built_area'] + ds['buildings'])/(1. - ds['cloud'] - ds['shadow'])
    derived['total_pervious'] = (ds['grass'] + ds['trees'] + ds['other_veg'] + ds['bare_earth'] + ds['water'] + ds['swimming_pool'])/(1. - ds['cloud'] - ds['shadow'])

    return derived

def assert_tests(geoscape,derived):

    '''test fractions sum to 1'''

    # test geoscape classess sum to 1
    total = np.full_like(geoscape.buildings.values,0)
    for key in geoscape.keys():
        total = total + geoscape[key]
    np.testing.assert_allclose(total.to_series().dropna().values, 1, 1E-10)

    # test derived classes sum to 1
    total = np.full_like(derived.building_fraction.values,0)
    for key in ['building_fraction','tree_fraction','lowveg_fraction','bareearth_fraction','roadpath_fraction','water_fraction']:
        total = total + derived[key]
    np.testing.assert_allclose(total.to_series().dropna().values, 1, 1E-10)

    # test total classes sum to 1
    total = np.full_like(derived.building_fraction.values,0)
    for key in ['total_built','total_pervious']:
        total = total + derived[key]
    np.testing.assert_allclose(total.to_series().dropna().values, 1, 1E-10)

    # test total classes sum to 1
    total = np.full_like(derived.building_fraction.values,0)
    for key in ['roadpath_fraction','total_pervious','building_fraction']:
        total = total + derived[key]
    np.testing.assert_allclose(total.to_series().dropna().values, 1, 1E-10)

    return

################################################################################

def main_calculate_morphology(res):
    
    '''main function for calculating grid-level morphology from geoscape building and tree data
    split_by_grid=False: building morphology information placed at centroid point of building for faster processing
    split_by_grid=True splits buildings at grid edges so morphology data sits in multiple cells'''

    ################################################################################
    # format individual building information

    print('reading building shapefile')
    raw = gpd.read_file(shp_fname)
    raw.columns= raw.columns.str.lower()

    if domain == 'melbourne':
        # drop buildings without height info
        raw = raw.rename(columns={'vic_heig_4':'bld_hgt'})
        raw = raw[raw.bld_hgt.notna()]
    else:
        print('calculating average of roof and eave height per building')
        raw['bld_hgt'] = (raw['roof_hgt']+raw['eave_hgt'])/2
    
    buildings = raw[['bld_pid','bld_hgt','geometry']]

    print('creating grid from template')
    fname = f'{outpath}/{domain.capitalize()}_geoscape_surface_cover_{res}mx{res}m_{output_v}.nc'
    template = xr.open_dataset(fname)
    cell = create_grid_from_template(template)

    # split buildings where intersecting with grid (warning, very long processing)
    if split_by_grid:
        spg='_split'
        # split buildings into cells by clipping
        print('splitting geometry by grid')
        total_list = []
        for shape in cell.geometry:
            bounds =  gpd.GeoDataFrame(crs='EPSG:4283', geometry=[shape])
            try:
                total_list.append(gpd.clip(buildings,bounds))
            except Exception:
                pass
        buildings = pd.concat(total_list)
    else:
        spg=''

    print('calculating building area and centroid')
    x,y = get_midpoint(buildings)
    new_proj = buildings['geometry'].set_crs(epsg=4283).to_crs(f'+proj=cea +lat_0={y} +lon_0={x} +units=m')
    buildings = buildings.assign(area=new_proj.area)

    print('calculating perimeter per building')
    buildings['perimeter'] = buildings.geometry.apply(get_perimeter)

    print('calculating wall area per building')
    buildings['wall_area'] = buildings['perimeter']*buildings['bld_hgt']

    print('calculating frontal area per building')
    x_frontal = buildings.apply(get_dist_x,axis=1)*buildings['bld_hgt']
    y_frontal = buildings.apply(get_dist_y,axis=1)*buildings['bld_hgt']
    buildings['avg_frontal'] = (x_frontal + y_frontal)/2

    print('saving geometry as geopackage')
    buildings.to_file(f'{projpath}/geopackages/buildings{spg}_{domain}.gpkg', driver='GPKG')

    ################################################################################
    # create gridded data from buildings information

    print('reading tree file')
    tree_orig = rxr.open_rasterio(tree_fname,masked=True).squeeze(['band'],drop=True)
    tree_orig = tree_orig.where(tree_orig>0)

    print('creating grid from template')
    fname = f'{outpath}/{domain.capitalize()}_derived_surface_cover_{res}mx{res}m_{output_v}.nc'
    template = xr.open_dataset(fname)
    cell = create_grid_from_template(template)

    ################################################################################

    print('getting weighted gridded values')
    df = calc_all_grids_from_template(buildings,cell,template,domain)

    print('creating netcdf')
    ds = create_ds(df)

    print('regridding tree height')
    ds['tree_height'] = regrid_raster_mean(tree_orig,template,res=res).round(4)
    ds['tree_height_std'] = regrid_raster_sd(tree_orig,template,res=res).round(4)

    print('include surface cover from geoscape data')
    ds['building_fraction']  = template['building_fraction']
    ds['tree_fraction']      = template['tree_fraction']
    ds['lowveg_fraction']    = template['lowveg_fraction']
    ds['water_fraction']     = template['water_fraction']
    ds['bareearth_fraction'] = template['bareearth_fraction']
    ds['roadpath_fraction']  = template['roadpath_fraction']
    ds['total_built']        = template['total_built']
    ds['total_pervious']     = template['total_pervious']

    print('calculating height to width ratio')
    # from masson 2020: h/w = lambda_w/(2*(1-lambda_p))
    ds['height_to_width'] = ds['wall_density']/(2*(1.-ds['building_fraction'])).round(4)

    # replace undefined with nan (where building fraction = 1)
    ds['height_to_width'].values = np.where(ds['height_to_width'].values==np.inf,np.nan,ds['height_to_width'].values)

    print('calculating skyview factor')
    # from masson 2020: [(h/w)^2 + 1]^(1/2) - h/w
    ds['skyview_factor'] = ((ds['height_to_width']**2 + 1.)**(1./2) - ds['height_to_width']).round(4)

    ds['displacement_mac'], ds['roughness_mac'] = calc_site_roughness(
        bldhgt_ave=ds['building_height'],
        sigma_h=ds['building_height_std'],
        lambda_p=ds['building_fraction'],
        lambda_pv=ds['tree_fraction'],
        lambda_f=ds['frontal_density'],
        bldhgt_max=ds['building_height_max'],
        mode=0) # macdonald et al 1998 staggered grid

    ds['displacement_kanda'], ds['roughness_kanda'] = calc_site_roughness(
        bldhgt_ave=ds['building_height'],
        sigma_h=ds['building_height_std'],
        lambda_p=ds['building_fraction'],
        lambda_pv=ds['tree_fraction'],
        lambda_f=ds['frontal_density'],
        bldhgt_max=ds['building_height_max'],
        mode=3) # Kanda et al 2013

    ds['displacement_mac'] = ds['displacement_mac'].round(4)
    ds['roughness_mac'] = ds['roughness_mac'].round(4)
    ds['displacement_mac'] = ds['displacement_mac'].round(4)
    ds['roughness_mac'] = ds['roughness_mac'].round(4)
    ds['displacement_kanda'] = ds['displacement_kanda'].round(4)
    ds['roughness_kanda'] = ds['roughness_kanda'].round(4)

    ################################################################################

    ds = set_attributes(ds,res)

    encoding = {var:{'zlib':True}for var in ds.data_vars}

    for key in ds.keys():
        ds[key].encoding['_FillValue'] = missing_float

    fname = f'{outpath}/{domain.capitalize()}_surface_{res}m_{output_v}'
    print(f'saving output to {fname}')
    ds.to_netcdf(f'{fname}.nc',encoding=encoding)
    
    try:
        ds.rio.write_crs("epsg:4326", inplace=True)
        ds.rio.to_raster(f'{fname}.tiff', compress='LZW')
    except Exception as e:
        print('raster tiff save error:')
        print(e)

    # get all data into dataframe in form for geopandas
    df_all = ds.to_dataframe().sort_values(['longitude','latitude'], 
        ascending=[True,False]).reset_index()

    gf = gpd.GeoDataFrame(df_all, 
            geometry=cell.geometry,crs='EPSG:4283')

    gf.to_file(f'{projpath}/geopackages/g{res}_{domain}.gpkg', driver='GPKG')

    return

def calc_all_grids_from_template(buildings,cell,template,domain='sample'):

    '''calculates grid-level morphology from building data using template grid '''

    print(f'calculating all grids using template')

    # set building morphology as point at centroid of individual buildings
    print('\nNOTE: Warning regarding centroid being likely incorrect can be ignored as the small scales being processed will have negligible error\n')
    points = gpd.GeoDataFrame(buildings, crs = f'EPSG:4283', geometry=gpd.points_from_xy(buildings.geometry.centroid.x, buildings.geometry.centroid.y))
    points['bld_hgt_wgt'] = points['bld_hgt'] * points['area']
    # spatial join building points within grid cells
    building_points = gpd.sjoin(points, cell, how='left', op='within')
    grp = building_points.groupby('index_right')

    gridded = pd.DataFrame()
    gridded['longitude'] = np.meshgrid(template.latitude.values,template.longitude.values)[1].flatten()
    gridded['latitude'] = np.meshgrid(template.latitude.values,template.longitude.values)[0].flatten()

    print('calculating building heights')
    gridded['building_height'] = grp['bld_hgt_wgt'].sum()/grp['area'].sum() # use weighted calculation so larger buildings in grid have greater effect
    gridded['building_height_max'] = grp['bld_hgt'].max()
    gridded['building_height_std'] = grp['bld_hgt'].std()

    gridded['wall_density'] = grp['wall_area'].sum()/cell['cell_area']
    gridded['frontal_density'] = grp['avg_frontal'].sum()/cell['cell_area']

    ## replace nan with zero
    # keys = ['wall_density','frontal_density']
    # gridded[keys] = gridded[keys].fillna(0)

    #################################################################
    # # use building points building_fraction for morphology calculations, RESULTS IN VERY HIGH H/W RATIO
    # gridded['building_fraction'] = grp['area'].sum()/cell['cell_area']
    # gridded['building_fraction'] = np.where(gridded['building_fraction'].values > 0.99,0.99,gridded['building_fraction'].values)

    return gridded

def create_grid_from_template(template,crs='EPSG:4283'):

    '''calculate shapely polygon for each grid cell using netcdf template '''

    cell_size_x = template.longitude.values[1] - template.longitude.values[0]
    cell_size_y = template.latitude.values[1] - template.latitude.values[0]

    # create the cells in a loop
    grid_cells = []
    for x0 in template.longitude.values+cell_size_x/2:
        for y0 in template.latitude.values-cell_size_y/2:
            # bounds
            x1 = x0-cell_size_x
            y1 = y0+cell_size_y
            grid_cells.append( shapely.geometry.box(x0, y0, x1, y1)  )
    cell = gpd.GeoDataFrame(grid_cells, columns=['geometry'],crs=crs)

    cell['cell_area'] = cell.apply(get_dist_x,axis=1)*cell.apply(get_dist_y,axis=1)
                                     
    return cell

def get_midpoint(gdf):

    # total area for the grid
    xmin, ymin, xmax, ymax= gdf.total_bounds
    xmid = (xmin+xmax)/2
    ymid = (ymin+ymax)/2

    return xmid,ymid

def get_perimeter(shp):

    # convert shapely polygon from 3D to 2D (ignore height)
    geom = shapely.wkb.loads(
        shapely.wkb.dumps(shp, output_dimension=2))

    # calculate perimeter along the geodesic Earth
    perim = cgeo.Geodesic().geometry_length(geom)

    return perim

def get_dist_x(row):
    '''building frontal wall area in x direction'''

    bounds = row['geometry'].bounds

    spoint = (bounds[0],bounds[1])
    epoint = (bounds[2],bounds[1])

    x = cgeo.Geodesic().inverse(spoint,epoint).base[0][0]

    return x

def get_dist_y(row):

    bounds = row['geometry'].bounds

    spoint = (bounds[0],bounds[1])
    epoint = (bounds[0],bounds[3])

    y = cgeo.Geodesic().inverse(spoint,epoint).base[0][0]

    return y

def regrid_raster_mean(rorig,ds,res,method='nearest'):

    orig_res = rorig.rio.resolution()[0]
    sfactor = int(res/orig_res)
    new = rorig.coarsen(x=sfactor,y=sfactor,boundary='pad').mean()
    new = new.rio.reproject(dst_crs='+proj=longlat',nodata=np.nan)
    # new = new.where(new>0)
    new = new.rename({'x':'longitude','y':'latitude'}).drop('spatial_ref')

    da = new.interp_like(ds,method=method)

    return da

def regrid_raster_sd(rorig,ds,res,method='nearest'):

    orig_res = rorig.rio.resolution()[0]
    sfactor = int(res/orig_res)
    new = rorig.coarsen(x=sfactor,y=sfactor,boundary='pad').std()
    new = new.rio.reproject(dst_crs='+proj=longlat',nodata=np.nan)
    # new = new.where(new>0)
    new = new.rename({'x':'longitude','y':'latitude'}).drop('spatial_ref')

    da = new.interp_like(ds,method=method)

    return da

def create_ds(df):

    '''convert pandas dataframe to xarray dataset'''

    ds = df.set_index(['latitude', 'longitude']).sort_index().to_xarray()

    # attributes to coordinates
    ds.latitude.attrs = {
            'axis':           'Y',
            'long_name':      'latitude',
            'standard_name':  'latitude',
            'units':          'degrees_north',
        }
    ds.longitude.attrs = {
            'axis':           'X',
            'long_name':      'Longitude',
            'standard_name':  'longitude',
            'units':          'degrees_east',
        }

    return ds

def set_attributes(ds,res):

    '''set dataset attributes '''

    ### surface cover ###
    if 'building_fraction' in ds.keys():
        ds['building_fraction'].attrs = {
                'long_name'   : 'building footprint fraction (lambda_p)',
                'units'       : '-',
                'description' : 'building footprint area as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'tree_fraction' in ds.keys():
        ds['tree_fraction'].attrs = {
                'long_name'   : 'tree canopy fraction (lambda_vt)',
                'units'       : '-',
                'description' : 'tree canopy plan area as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'lowveg_fraction' in ds.keys():
        ds['lowveg_fraction'].attrs = {
                'long_name'   : 'grass and low vegetation fraction (lambda_vl)',
                'units'       : '-',
                'description' : 'low vegetation (grass, shrubs, other vegetation) as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'water_fraction' in ds.keys():
        ds['water_fraction'].attrs = {
                'long_name'   : 'open water fraction (lambda_wa)',
                'units'       : '-',
                'description' : 'all open water (ocean, lakes, pools) as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'bareearth_fraction' in ds.keys():
        ds['bareearth_fraction'].attrs = {
                'long_name'   : 'bare earth, sand and rock fraction (lambda_be)',
                'units'       : '-',
                'description' : 'bare earth including construction sites, rock, sand and sparsely vegetated areas as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'roadpath_fraction' in ds.keys():
        ds['roadpath_fraction'].attrs = {
                'long_name'   : 'road, path and other hard surfaces on ground (lambda_if)',
                'units'       : '-',
                'description' : 'all hard surfaces on ground excluding buildings, defined as "impervious surface fraction" in Stewart and Oke, 2012, corrected for cloud and shadow fractions'
            }
    if 'total_built' in ds.keys():
        ds['total_built'].attrs = {
                'long_name'   : 'total built fraction (lambda_tb)',
                'units'       : '-',
                'description' : 'all impervious surfaces including buildings, roads, paths and other hard surfaces, corrected for cloud and shadow fractions'
            }
    if 'total_pervious' in ds.keys():
        ds['total_pervious'].attrs = {
                'long_name'   : 'total pervious fraction (lambda_tp)',
                'units'       : '-',
                'description' : 'all pervious surfaces including vegetation, water and bare earth, corrected for cloud and shadow fractions'
            }

    # morphology
    if 'frontal_density' in ds.keys():
        ds['frontal_density'].attrs = {
                'long_name'   : 'frontal area density (lambda_f)',
                'units'       : '-',
                'description' : 'sum of cardinally averaged building frontal area as fraction of grid area'
            }
    if 'wall_density' in ds.keys():
        ds['wall_density'].attrs = {
                'long_name'   : 'wall area density (lambda_w)',
                'units'       : '-',
                'description' : 'sum of building wall area as fraction of grid area'
            }

    if 'building_height' in ds.keys():
        ds['building_height'].attrs = {
                'long_name'   : 'building height mean (H_avg)',
                'units'       : 'm',
                'description' : 'mean building height in grid cell (avg. of geoscape roof and eave height)'
            }
    if 'building_height_max' in ds.keys():
        ds['building_height_max'].attrs = {
                'long_name'   : 'building height maximum (H_max)',
                'units'       : 'm',
                'description' : 'maximum building height in grid cell (avg. of geoscape roof and eave height)'
            }
    if 'building_height_std' in ds.keys():
        ds['building_height_std'].attrs = {
                'long_name'   : 'building height standard deviation (H_std)',
                'units'       : 'm',
                'description' : 'standard deviation of building height in grid cell (avg. of geoscape roof and eave height)'
            }

    if 'tree_height' in ds.keys():
        ds['tree_height'].attrs = {
                'long_name'   : 'tree height average (HT_ave)',
                'units'       : 'm',
                'description' : 'average vegetation canopy height in grid'
            }
    if 'tree_height_std' in ds.keys():
        ds['tree_height_std'].attrs = {
                'long_name'   : 'tree height standard deviation (HT_ave)',
                'units'       : 'm',
                'description' : 'standard deviation in vegetation canopy height in grid'
            }

    if 'roughness_mac' in ds.keys():
        ds['roughness_mac'].attrs = {
                'long_name'   : 'aerodynamic roughness length (Z_0,mac)',
                'units'       : 'm',
                'description' : 'roughness length for staggered arrays from Macdonald et al., 1998: https://doi.org/10.1016/S1352-2310(97)00403-2'
            }
    if 'displacement_mac' in ds.keys():
        ds['displacement_mac'].attrs = {
                'long_name'   : 'displacement height (Z_d,mac)',
                'units'       : 'm',
                'description' : 'zero-plane displacement height from Macdonald et al., 1998: https://doi.org/10.1016/S1352-2310(97)00403-2'
            }
    if 'roughness_kanda' in ds.keys():
        ds['roughness_kanda'].attrs = {
                'long_name'   : 'aerodynamic roughness length (Z_0,kan)',
                'units'       : 'm',
                'description' : 'roughness length for staggered arrays from Kanda et al., 2013: https://doi.org/10.1007/s10546-013-9818-x'
            }
    if 'displacement_kanda' in ds.keys():
        ds['displacement_kanda'].attrs = {
                'long_name'   : 'displacement height (Z_d,kan)',
                'units'       : 'm',
                'description' : 'zero-plane displacement height from Kanda et al., 2013: https://doi.org/10.1007/s10546-013-9818-x'
            }

    ### canyon attributes ###
    if 'height_to_width' in ds.keys():
        ds['height_to_width'].attrs = {
                'long_name'   : 'canyon height to width aspect ratio (h/w)',
                'units'       : '-',
                'description' : 'average aspect ratio assuming street canyon geometry using Eq 1 of Masson et al. 2020: https://doi.org/10.1016/j.uclim.2019.100536'
            }
    if 'skyview_factor' in ds.keys():
        ds['skyview_factor'].attrs = {
                'long_name'   : 'skyview factor (psi)',
                'units'       : '-',
                'description' : 'average skyview factor assuming street canyon geometry using Eq 2 of Masson et al. 2020: https://doi.org/10.1016/j.uclim.2019.100536'
            }

    # other
    if 'cell_area' in ds.keys():
        ds['cell_area'].attrs = {
                'long_name'   : 'cell area',
                'units'       : 'm2',
                'description' : f'area of grid cell, nominally at {res} m resolution'
            }
    if 'bld_hgt' in ds.keys():
        ds['bld_hgt'].attrs = {
                'long_name'   : 'building height average',
                'units'       : 'm',
                'description' : 'average for midpoint of roof and eave height in grid (no buiding area weighting)'
            }
    if 'bld_hgt_wgt' in ds.keys():
        ds['bld_hgt_wgt'].attrs = {
                'long_name'   : 'building height weighted average',
                'units'       : 'm',
                'description' : 'area weighted average for midpoint of roof and eave height in grid'
            }
    if 'roof_hgt' in ds.keys():
        ds['roof_hgt'].attrs = {
                'long_name'   : 'roof height average',
                'units'       : 'm',
                'description' : 'average roof height in grid (no buiding area weighting)'
            }
    if 'roof_hgt_wgt' in ds.keys():
        ds['roof_hgt_wgt'].attrs = {
                'long_name'   : 'roof height weighted average',
                'units'       : 'm',
                'description' : 'area weighted average roof height in grid'
            }
    if 'eave_hgt' in ds.keys():
        ds['eave_hgt'].attrs = {
                'long_name'   : 'eave height average',
                'units'       : 'm',
                'description' : 'average eave height in grid (no buiding area weighting)'
            }
    if 'eave_hgt_wgt' in ds.keys():
        ds['eave_hgt_wgt'].attrs = {
                'long_name'   : 'eave height weighted average',
                'units'       : 'm',
                'description' : 'area weighted average eave height in grid'
            }

    ds.attrs['title']       = f'Urban characteristics for {domain.capitalize()} at {res} m resolution derived from Geoscape datasets'
    ds.attrs['version']     = output_v
    ds.attrs['institution'] = 'ARC Centre of Excellence for Climate Extremes, UNSW Sydney, Australia'
    ds.attrs['source']      = 'Derived from Geoscape datasets, including Buildings v2.0, Trees v1.6 and Surface cover v1.6: PSMA Australia, 2020. https://geoscape.com.au/'
    ds.attrs['reference']  = 'A transformation in city-descriptive input data for urban climate models'
    ds.attrs['licence']  = 'The contents of this file are restricted, not to be used or distributed without permission from the authors.'

    return ds

def calc_site_roughness(bldhgt_ave,sigma_h,lambda_p,lambda_pv,lambda_f,bldhgt_max=None,mode=0):
    '''estimate urban site roughness length for momentum and zero-plane displacement based on various methods 
    as described in Grimmond and Oke (1999): Aerodynamic Properties of Urban Areas Derived from Analysis of Surface Form .

    modes: four different morphometric methods to calculate roughness and displacement:
        - 0: Macdonald et al. 1998: An improved method for the estimation of surface roughness of obstacle arrays
        - 1: Kent et al 2017: Aerodynamic roughness parameters in cities: Inclusion of vegetation
        - 2: Millward-Hopkins et al., (2011) per Kent et al. 2017 eq 11 and 12 # TYPO IN EQ 12 OF KENT PAPER
        - 3: Kanda et al., 2013

    Inputs
    ------
    bldhgt_ave [m] : building mean height
    sigma_h [m]    : building height standard deviation
    lambda_p [0-1] : building plan area fraction
    lambda_pv [0-1]: tree plan area fraction
    lambda_f [1]   : wall frontal area fraction
    bldhgt_max [m] : maximum building height (if none assume 1.25 times SD)

    Returns
    -------
    zd [m] zero-plane displacement
    z0 [m] rougness length for momentum

    Warning
    -------
    in mode=0 wall frontal area (assuming random canyon orientation) calculated from Porsen et al 2010 (DOI:10.1002/qj.668)
    in mode=1 (including vegetation), vegetation frontal area index assumed to be equal to tree plan fraction (from Table 1 in Kent et al 2017)
    in mode=3 (Kanda) maximum building height assumed to be 2 standard deviations
    '''

    vonkarm = 0.4
    dragbld = 1.2

    if mode==0:
        print('zd,z0 from Macdonald et al. 1998')

        alpha   = 4.43 # for staggered arrays
        beta    = 1.0  # for staggered arrays

        zd = (1. + alpha**(-lambda_p)*(lambda_p - 1.))*bldhgt_ave    # eq 23 from Mac1998
        z0 = ((1. - zd/bldhgt_ave) * np.exp( -1.*(0.5*beta*(dragbld/vonkarm**2)*(1.-zd/bldhgt_ave)*lambda_f)**(-0.5) ))*bldhgt_ave  # eq 26 from Mac1998

    if mode == 1: # with veg from Kent et al (2018) https://doi.org/10.1007/s11252-017-0710-1
        print('zd,z0 from Macdonald with vegetation from Kent et al. 2017 and 2018')
        # Aerodynamic roughness variation with vegetation: analysis in a suburban neighbourhood and a city park

        alpha   = 4.43 # for staggered arrays
        beta    = 1.0  # for staggered arrays

        lambda_fv = lambda_pv # estimated equality based on Table 1 of Kent et al 2017 for different sites
        P_3D = 0.4   # leaf-on porosity

        lambda_tot = (lambda_p + lambda_pv*(1. - P_3D)) # eq 4 from Kent et al (2018)
        Pv  = (-1.251*P_3D**2 + 0.489*P_3D + 0.803)/dragbld  # eq 7 from Kent et al (2018)
        weighted_frontal_area = lambda_f + lambda_fv*Pv

        zd = (1. + alpha**(-lambda_tot)*(lambda_tot-1.))*bldhgt_ave # eq 5
        z0 = (1. - zd/bldhgt_ave)*np.exp( -1.*(0.5*beta*(dragbld/vonkarm**2)*(1.-zd/bldhgt_ave)*weighted_frontal_area)**(-0.5) )*bldhgt_ave # eq 6

    if mode == 2: # Millward-Hopkins et al. (2011), per Kent et al. 2017 eq 11 and 12
        print('zd,z0 from Millward-Hopkins et al. 2011, per Kent et al. 2017 eq 11 and 12')

        MhoUzd_on_Hav_A = (19.2*lambda_p - 1. + np.exp(-19.2*lambda_p))/(19.2*lambda_p*(1.- np.exp(-19.2*lambda_p)))
        MhoUzd_on_Hav_B = (117*lambda_p + (187.2*lambda_p**3 - 6.1)*(1.- np.exp(-19.2*lambda_p)) )/( (1.+114*lambda_p + 187*lambda_p**3)*(1.-np.exp(-19.2*lambda_p)) )
        MhoUzd_on_Hav = np.where( lambda_p >= 0.19, MhoUzd_on_Hav_A,  MhoUzd_on_Hav_B )

        Mho_exp = np.exp( -1.*(0.5*dragbld*vonkarm**(-2)*lambda_f)**(-0.5) )
        MhoUz0_on_Hav = ( (1. - MhoUzd_on_Hav)*Mho_exp)

        zd = bldhgt_ave*(MhoUzd_on_Hav + ((0.2375 * np.log(lambda_p) + 1.1738)*(sigma_h/bldhgt_ave)) )

        # z0 = bldhgt_ave*(MhoUz0_on_Hav + (np.exp(0.8867*lambda_f) - 1.)*(sigma_h/bldhgt_ave)**(np.exp(2.3271*lambda_f))) # ERROR IN KENT ET AL PAPER
        z0 = bldhgt_ave*(MhoUz0_on_Hav + np.exp(0.8867*lambda_f - 1.)*(sigma_h/bldhgt_ave)**(np.exp(2.3271*lambda_f))) # correction based on UMEP

    if mode == 3:
        print('zd,z0 Kanda et al 2013')

        alpha   = 4.43 # for staggered arrays
        beta    = 1.0  # for staggered arrays

        if np.all(bldhgt_max)==None:
            bldhgt_max = 1.5*sigma_h + bldhgt_ave # assume bldhgt_max
            # bldhgt_max = 12.51*sigma_h**0.77 # eq 3 from Kanda et al 2013

        # first calculate macdonald et al 1998 (for later scaling)
        mac_zd = (1. + (lambda_p - 1.)*alpha**(-lambda_p))*bldhgt_ave
        mac_z0 = ((1. - mac_zd/bldhgt_ave) * np.exp( -1.*(0.5*beta*(dragbld/vonkarm**2)*(1.-mac_zd/bldhgt_ave)*lambda_f)**(-0.5) ))*bldhgt_ave

        # then Kanda et al 2013 scaling parameters
        a0,b0,c0 = 1.29, 0.36, -0.17
        a1,b1,c1 = 0.71, 20.21, -0.77

        X = (sigma_h + bldhgt_ave)/bldhgt_max  # eq. 10b
        Y = (lambda_p*sigma_h)/bldhgt_ave        # eq. 12b

        zd = bldhgt_max*(c0*X**2 + (a0*lambda_p**b0 - c0)*X) # eq 10a 
        z0 = mac_z0*(a1 + b1*Y**2 + c1*Y) # eq 12a

    zd = np.round(zd,4)
    zd = xr.where(zd<0,0,zd)
    z0 = np.round(z0,4)
    z0 = xr.where(z0<0,0,z0)

    return zd,z0

################################################################################

def main_plot_sample():

    buildings = gpd.read_file(f'{projpath}/geopackages/buildings_{domain}.gpkg')

    gf1,ds1 = load_files(res_list[0])
    gf2,ds2 = load_files(res_list[1])
    gf3,ds3 = load_files(res_list[2])

    key1='bld_hgt'
    keys=['building_height','building_height_max','building_height_std','tree_height',
            'wall_density','building_fraction','roadpath_fraction','total_pervious','total_built',
            'height_to_width','skyview_factor','roughness_mac','displacement_mac']

    for key2 in keys:
        print(f'plotting {key2}')
        plot_geoscape_and_grid(buildings,gf1,gf2,gf3,key1,key2)

    return

def load_files(res):

    ''' read building and derived data for plotting at particular resolution'''

    gf = gpd.read_file(f'{projpath}/geopackages/g{res}_{domain}.gpkg')
    ds = xr.open_dataset(f'{outpath}/{domain.capitalize()}_surface_{res}m_{output_v}.nc')

    return gf,ds

def plot_geoscape_and_grid(gdf,gf1,gf2,gf3,key1,key2):

    '''creat plot comparing original and derived (lower resolution) data'''
    
    print(f'plotting {key1} and {key2}')

    plt.close('all')
    fig, axes = plt.subplots(nrows=2,ncols=2,figsize=(9,6),sharey=True,sharex=True)

    ax = axes.flatten()

    plot_orig(gdf, ax[0],key1,cmap='plasma')
    plot_grid(gf1, ax[1],key2,res=f'{res_list[0]}m',cmap='plasma')
    plot_grid(gf2, ax[2],key2,res=f'{res_list[1]}m',cmap='plasma')
    plot_grid(gf3, ax[3],key2,res=f'{res_list[2]}m',cmap='plasma')

    fig.subplots_adjust(left=0.05 ,right=0.92,bottom=0.05,top=0.95,wspace=0.1,hspace=0.1)
    fig.savefig(f'{plotpath}/res_{key1}_{key2}_{output_v}.png', dpi=400, bbox_inches='tight')
    # plt.show()

    return

def plot_orig(gdf,ax,key,cmap='plasma'):
    ''' plot original geopandas with colorbar (no edgecolor)'''

    un,tt,vmax = get_labels(key)


    ax.set_title(f'Original data: {tt}', fontsize=11)
    div = make_axes_locatable(ax)
    cax = make_axes_locatable(ax).append_axes('right', size="3%", pad=-0.15)
    cax.text(1.10,1.06,un,fontsize=8,va='bottom',ha='left',transform=cax.transAxes)

    gdf.plot(column=key, ax=ax, cax=cax, vmax=vmax, legend=True, legend_kwds={'label':'', 'extend':'max'}, cmap=cmap)
    ax.tick_params(axis='both', labelsize=9)
    ax.locator_params(axis='x', nbins=5)
    ax.locator_params(axis='y', nbins=6)

    ax.text(0.01,0.99,f'min: {gdf[key].min():.2f}', fontsize=6, ha='left', va='top', transform=ax.transAxes)
    ax.text(0.01,0.95,f'max: {gdf[key].max():.2f}', fontsize=6, ha='left', va='top', transform=ax.transAxes)

    return ax

def plot_grid(gdf,ax,key,res,cmap='plasma'):
    ''' plot gridded geopandas with colorbar (with edgecolor)'''

    un,tt,vmax = get_labels(key)

    ax.set_title(f'Derived: {res} {tt}', fontsize=11)
    div = make_axes_locatable(ax)
    cax = make_axes_locatable(ax).append_axes('right', size="3%", pad=-0.15)
    cax.text(1.10,1.06,un,fontsize=8,va='bottom',ha='left',transform=cax.transAxes)

    gdf.plot(column=key, ax=ax, cax=cax, vmax=vmax, cmap=cmap, 
        legend=True, legend_kwds={'label': '', 'extend':'max'}, edgecolor='0.8',linewidth=0.3)
    ax.tick_params(axis='both', labelsize=9)
    ax.locator_params(axis='x', nbins=5)
    ax.locator_params(axis='y', nbins=6)

    ax.text(0.01,0.99,f'min: {gdf[key].min():.2f}', fontsize=6, ha='left', va='top', transform=ax.transAxes)
    ax.text(0.01,0.95,f'max: {gdf[key].max():.2f}', fontsize=6, ha='left', va='top', transform=ax.transAxes)

    return ax

def get_labels(key):

    '''get key labels for plotting '''

    units=''
    title = ''
    vmax=None

    if key=='area':
        title = 'building footprint area'
        units = '[$m^2$]'
    elif key=='geometry':
        title = 'building footprint'
        units = '[-]'
    elif key == 'bld_hgt':
        title = 'building height'
        units = '[m]'
        vmax = 100
    elif key == 'building_height':
        title = 'mean building height'
        units = '[m]'
        vmax = 30
    elif key == 'building_height_max':
        title = 'maximum building height'
        units = '[m]'
        vmax = 100
    elif key == 'building_height_std':
        title = 'standard deviation of building height'
        units = '[m]'
        vmax = 10
    elif key == 'tree_height':
        title = 'mean tree height'
        units = '[m]'
        vmax = 30
    elif key == 'wall_density':
        title = 'wall area density'
        units = '[-]'
        vmax = 2
    elif key == 'frontal_density':
        title = 'frontal area density'
        units = '[-]'
        vmax = 1
    elif key == 'building_fraction':
        title = 'building footprint fraction'
        units = '[-]'
        vmax = 1
    elif key == 'roadpath_fraction':
        title = 'roads, paths, hard surfaces fraction'
        units = '[-]'
        vmax = 1
    elif key == 'tree_fraction':
        title = 'tree fraction'
        units = '[-]'
        vmax = 1
    elif key == 'total_pervious':
        title = 'total pervious fraction'
        units = '[-]'
        vmax = 1
    elif key == 'total_vegetation':
        title = 'total vegetation fraction'
        units = '[-]'
        vmax = 1
    elif key == 'total_built':
        title = 'total built fraction'
        units = '[-]'
        vmax = 1
    elif key == 'skyview_factor':
        title = 'skyview factor'
        units = '[-]'
        vmax = 1
    elif key == 'height_to_width':
        title = 'canyon aspect ratio'
        units = '[-]'
        vmax = 2.5
    elif key == 'roughness_mac':
        title = 'roughness length (Mac 1998)'
        units = '[m]'
        vmax = 5
    elif key == 'displacement_mac':
        title = 'displacement height (Mac 1998)'
        units = '[m]'
        vmax = 10
    else:
        print('WARNING: key not recognised')

    return units,title,vmax

################################################################################

if __name__ == "__main__":

    ################################################################################
    # path setup

    for mkpath in [plotpath,outpath,f'{projpath}/geopackages']:
        if not os.path.exists(mkpath):
            os.makedirs(mkpath)
            print("created folder : ", mkpath)

    if domain == 'sample':
        tif_fpaths = [f'{datapath}/Geoscape_Sample_data_2021/06_Surface_Cover_Trees/NSW_SURFACECOVER_2M_Z56_16837.tif']
        tree_fname = f'{datapath}/Geoscape_Sample_data_2021/06_Surface_Cover_Trees/NSW_TREES_Z56_16837.tif'
        shp_fname =  f'{datapath}/Geoscape_Sample_data_2021/01_Buildings/SHP/buildings.shp'

    ################################################################################

    main()


__title__ = "Code for resampling Geoscape building and surface cover data"
__version__ = "2022-28-01"
__author__ = "<NAME>"
__email__ = "<EMAIL>"

'''
Associated with the manuscript: A transformation in city-descriptive input data for urban climate models
Manuscript authors:  <NAME>, <NAME>,  <NAME>,  <NAME>, <NAME>

Copyright 2022 <NAME>

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

Instructions: 
1. Get sample data from Geoscape: https://geoscape.com.au/get-sample/
2. Install necessary python dependancies
3. Update project and data paths as required
4. Run script
'''

import os
import numpy as np
import pandas as pd
import xarray as xr
import rioxarray as rxr
from rioxarray import merge

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

import shapely
import geopandas as gpd
import cartopy.geodesic as cgeo

################################################################################
# user input

domain = 'sample'
output_v = 'v0.97'

projpath = f'.'                   # update path to working folder
datapath = f'{projpath}/data'     # update path to geoscape data
outpath = f'{projpath}/outputs'   # netcdf/tif outputs
plotpath = f'{projpath}/figures'  # figures output

res_list = [30,100,300]         # list of resolutions to process in m

split_by_grid = False   # if true, splits buildings on grid edges (long processing)
plot_comp = True        # if true, plots comparison between different resolutions (for sample)

missing_float = -999.

################################################################################

def main():

    for res in res_list:
        print(f'processing {res} m')

        main_resample_cover(res)
        main_calculate_morphology(res)

    if plot_comp:
        main_plot_sample()

    return

def main_resample_cover(res):
    '''main function for resampling geoscape surface cover data to a lower resolution and recategorise'''

    print(f'opening {domain} geotif')

    # read tif
    tif_list = []
    for tif_fpath in tif_fpaths:
        tif_list.append(rxr.open_rasterio(tif_fpath))
    print('merging tifs')
    rorig =  merge.merge_arrays(tif_list)

    # get band 1
    orig_ds = rorig.drop(['band']).squeeze()
    # remove fill value
    orig_ds = orig_ds.where(orig_ds!=orig_ds._FillValue)

    ################################################################################

    print('creating downscaled geoscape dataset')
    geoscape = create_dataset_from_geoscape(orig_ds,res)

    ################################################################################

    print('saving geoscape dataset')
    fname = f'{outpath}/{domain.capitalize()}_geoscape_surface_cover_{res}mx{res}m_{output_v}' 

    encoding = {var:{'zlib':True} for var in geoscape.data_vars}
    geoscape.to_netcdf(f'{fname}.nc',format='NETCDF4', encoding=encoding)

    geoscape = xr.open_dataset(f'{fname}.nc')
    # geoscape.rio.write_crs("epsg:4326", inplace=True)
    # geoscape.rio.to_raster(f'{fname}.tiff', compress='LZW')

    ################################################################################

    derived = calc_derived_fractions(geoscape,res)
    derived = set_attributes(derived,res)

    print('saving derived surface cover')

    fname = f'{outpath}/{domain.capitalize()}_derived_surface_cover_{res}mx{res}m_{output_v}'

    compression = {'zlib':True}
    encoding = {var:compression for var in derived.data_vars}

    derived.to_netcdf(f'{fname}.nc',format='NETCDF4', encoding=encoding)
    # derived.rio.write_crs("epsg:4326", inplace=True)
    # derived.rio.to_raster(f'{fname}.tiff', compress='LZW')

    ################################################################################

    assert_tests(geoscape,derived)

    return

def calc_fraction(key,ds,res,long_name=None):
    '''function to take categorical band and calculate its fraction at lower resolution'''

    print(f'calculating {key} fraction at {res} m')

    # get resolution factors
    orig_res = ds.rio.resolution()[0]
    sfactor = int(res/orig_res)

    da = ds.where(ds == key)

    # get all key values, convert to 1
    tmp = da.copy()
    tmp.values = np.where(da.values==key,1,np.nan)
    # calculate sums over coarsened grid
    tmp_sum = tmp.coarsen(x=sfactor,y=sfactor,boundary='pad').sum()

    # get all non-nan values, convert to 1
    tmp = da.copy()
    tmp.values = np.where(np.isnan(ds).values,ds.values,1)
    # calculate sums over coarsened grid
    tmp_all = tmp.coarsen(x=sfactor,y=sfactor,boundary='pad').sum()

    result = tmp_sum/tmp_all
    result = result.rename(str(key))

    result.attrs['long_name'] = long_name
    result.encoding['_FillValue'] =  missing_float

    return result

def create_dataset_from_geoscape(orig_ds,res):

    '''create netcdf attributes from geoscape surface cover data

    # Table 2: SURFACE_COVER_2M_TYPE_AUT Codes from 'Surface Cover 1.6 Product Description.pdf'

    2 Bare Earth                Includes sand dunes, desert, rock outcrops, bare soil other than bare agricultural land, and sparsely vegetated areas of grass and shrub. Non-vegetated strip mines and quarries except where covered by development or water.
    3 Road and Path             Roads and parking lots covered in a man-made material excluding hard packed dirt trails.
    4 Grass                     Grass and herbaceous areas. The category may include herbaceous wetlands if images are collected during dry season or periods of drought.
    5 Trees                     All trees including deciduous and evergreen woody vegetation.
    6 Unspecified Vegetation    Any other vegetative material not included within the Grass or Tree class. This may include, but is not limited to, shrub, sc rub, agriculture, and aquatic plants.
    7 Built-up Areas            Any areas of man-made environments and infrastructure excluding road and paths and buildings.
    8 Water                     Depending on the resolution quality of the imagery used, natural water will include streams, canals, ponds, lakes, reservoirs, estuaries and bays.
    9 Buildings                 Where the majority of a pixel intersects a Building, vector building polygon representation. 
    10 Cloud                    The area covered with cloud on Date of collection.
    11 Shadow                   The area covered with shadow on Date/time of collection.
    12 Swimming Pool            An area identified as a swimming pool.
    '''
    
    print('creating dataset')

    ds = xr.Dataset({

        'bare_earth': calc_fraction(2,orig_ds,res=res,
            long_name= 'Bare earth: includes bare soil, sand, desert, rock outcrops'),

        'road_path': calc_fraction(3,orig_ds,res=res,
            long_name = 'Road and path: includes roads, parking lots and hard packed dirt trails'),

        'grass': calc_fraction(4,orig_ds,res=res,
            long_name = 'Grass and herbaceous areas'),

        'trees': calc_fraction(5,orig_ds,res=res,
            long_name = 'All trees including deciduous and evergreen'),

        'other_veg': calc_fraction(6,orig_ds,res=res,
            long_name = 'Other vegetation including shrub, scrub, agriculture and aquatic plants'),

        'built_area': calc_fraction(7,orig_ds,res=res,
            long_name = 'Human-made environments and infrastructure excluding road, paths and buildings'),

        'water': calc_fraction(8,orig_ds,res=res,
            long_name = 'Open water including streams, canals, lakes, reservoirs and bays'),

        'buildings': calc_fraction(9,orig_ds,res=res,
            long_name = 'Buildings'),

        'cloud': calc_fraction(10,orig_ds,res=res,
            long_name = 'Area covered with cloud at time of collection'),

        'shadow': calc_fraction(11,orig_ds,res=res,
            long_name = 'Area covered with shadow at time of collection'),

        'swimming_pool': calc_fraction(12,orig_ds,res=res,
            long_name = 'Swimming pool'),
        })

    ds.attrs['title'] = f'{domain.capitalize()} surface cover at {res}m resolution derived from 2m resolution GeoScape product'
    ds.attrs['version'] = output_v
    ds.attrs['authors'] = '<NAME> <<EMAIL>>'
    ds.attrs['institution'] = 'ARC Centre of Excellence for Climate Extremes, UNSW Sydney, Australia'
    ds.attrs['source'] = 'Geoscape Surface Cover v1.6. PSMA Australia, 2020. https://geoscape.com.au/'
    ds.attrs['comment'] = 'The contents of this file are restricted, not to be used or distributed without permission from the authors.'

    geoscape = ds.rio.reproject(dst_crs=4326)
    geoscape = geoscape.rename({'y':'latitude','x':'longitude'}).drop('spatial_ref')

    return geoscape

def calc_derived_fractions(ds,res):

    '''create derived classes and attributes (normalising for cloud and shadow fractions)'''

    derived = xr.Dataset()

    derived['building_fraction'] = (ds['buildings'])/(1. - ds['cloud'] - ds['shadow'])
    derived['tree_fraction'] = (ds['trees'])/(1. - ds['cloud'] - ds['shadow'])
    derived['lowveg_fraction'] = (ds['grass'] + ds['other_veg'])/(1. - ds['cloud'] - ds['shadow'])
    derived['bareearth_fraction'] = (ds['bare_earth'])/(1. - ds['cloud'] - ds['shadow'])
    derived['roadpath_fraction'] = (ds['road_path'] + ds['built_area'])/(1. - ds['cloud'] - ds['shadow'])
    derived['water_fraction'] = (ds['water'] + ds['swimming_pool'])/(1. - ds['cloud'] - ds['shadow'])
    derived['total_built'] = (ds['road_path'] + ds['built_area'] + ds['buildings'])/(1. - ds['cloud'] - ds['shadow'])
    derived['total_pervious'] = (ds['grass'] + ds['trees'] + ds['other_veg'] + ds['bare_earth'] + ds['water'] + ds['swimming_pool'])/(1. - ds['cloud'] - ds['shadow'])

    return derived

def assert_tests(geoscape,derived):

    '''test fractions sum to 1'''

    # test geoscape classess sum to 1
    total = np.full_like(geoscape.buildings.values,0)
    for key in geoscape.keys():
        total = total + geoscape[key]
    np.testing.assert_allclose(total.to_series().dropna().values, 1, 1E-10)

    # test derived classes sum to 1
    total = np.full_like(derived.building_fraction.values,0)
    for key in ['building_fraction','tree_fraction','lowveg_fraction','bareearth_fraction','roadpath_fraction','water_fraction']:
        total = total + derived[key]
    np.testing.assert_allclose(total.to_series().dropna().values, 1, 1E-10)

    # test total classes sum to 1
    total = np.full_like(derived.building_fraction.values,0)
    for key in ['total_built','total_pervious']:
        total = total + derived[key]
    np.testing.assert_allclose(total.to_series().dropna().values, 1, 1E-10)

    # test total classes sum to 1
    total = np.full_like(derived.building_fraction.values,0)
    for key in ['roadpath_fraction','total_pervious','building_fraction']:
        total = total + derived[key]
    np.testing.assert_allclose(total.to_series().dropna().values, 1, 1E-10)

    return

################################################################################

def main_calculate_morphology(res):
    
    '''main function for calculating grid-level morphology from geoscape building and tree data
    split_by_grid=False: building morphology information placed at centroid point of building for faster processing
    split_by_grid=True splits buildings at grid edges so morphology data sits in multiple cells'''

    ################################################################################
    # format individual building information

    print('reading building shapefile')
    raw = gpd.read_file(shp_fname)
    raw.columns= raw.columns.str.lower()

    if domain == 'melbourne':
        # drop buildings without height info
        raw = raw.rename(columns={'vic_heig_4':'bld_hgt'})
        raw = raw[raw.bld_hgt.notna()]
    else:
        print('calculating average of roof and eave height per building')
        raw['bld_hgt'] = (raw['roof_hgt']+raw['eave_hgt'])/2
    
    buildings = raw[['bld_pid','bld_hgt','geometry']]

    print('creating grid from template')
    fname = f'{outpath}/{domain.capitalize()}_geoscape_surface_cover_{res}mx{res}m_{output_v}.nc'
    template = xr.open_dataset(fname)
    cell = create_grid_from_template(template)

    # split buildings where intersecting with grid (warning, very long processing)
    if split_by_grid:
        spg='_split'
        # split buildings into cells by clipping
        print('splitting geometry by grid')
        total_list = []
        for shape in cell.geometry:
            bounds =  gpd.GeoDataFrame(crs='EPSG:4283', geometry=[shape])
            try:
                total_list.append(gpd.clip(buildings,bounds))
            except Exception:
                pass
        buildings = pd.concat(total_list)
    else:
        spg=''

    print('calculating building area and centroid')
    x,y = get_midpoint(buildings)
    new_proj = buildings['geometry'].set_crs(epsg=4283).to_crs(f'+proj=cea +lat_0={y} +lon_0={x} +units=m')
    buildings = buildings.assign(area=new_proj.area)

    print('calculating perimeter per building')
    buildings['perimeter'] = buildings.geometry.apply(get_perimeter)

    print('calculating wall area per building')
    buildings['wall_area'] = buildings['perimeter']*buildings['bld_hgt']

    print('calculating frontal area per building')
    x_frontal = buildings.apply(get_dist_x,axis=1)*buildings['bld_hgt']
    y_frontal = buildings.apply(get_dist_y,axis=1)*buildings['bld_hgt']
    buildings['avg_frontal'] = (x_frontal + y_frontal)/2

    print('saving geometry as geopackage')
    buildings.to_file(f'{projpath}/geopackages/buildings{spg}_{domain}.gpkg', driver='GPKG')

    ################################################################################
    # create gridded data from buildings information

    print('reading tree file')
    tree_orig = rxr.open_rasterio(tree_fname,masked=True).squeeze(['band'],drop=True)
    tree_orig = tree_orig.where(tree_orig>0)

    print('creating grid from template')
    fname = f'{outpath}/{domain.capitalize()}_derived_surface_cover_{res}mx{res}m_{output_v}.nc'
    template = xr.open_dataset(fname)
    cell = create_grid_from_template(template)

    ################################################################################

    print('getting weighted gridded values')
    df = calc_all_grids_from_template(buildings,cell,template,domain)

    print('creating netcdf')
    ds = create_ds(df)

    print('regridding tree height')
    ds['tree_height'] = regrid_raster_mean(tree_orig,template,res=res).round(4)
    ds['tree_height_std'] = regrid_raster_sd(tree_orig,template,res=res).round(4)

    print('include surface cover from geoscape data')
    ds['building_fraction']  = template['building_fraction']
    ds['tree_fraction']      = template['tree_fraction']
    ds['lowveg_fraction']    = template['lowveg_fraction']
    ds['water_fraction']     = template['water_fraction']
    ds['bareearth_fraction'] = template['bareearth_fraction']
    ds['roadpath_fraction']  = template['roadpath_fraction']
    ds['total_built']        = template['total_built']
    ds['total_pervious']     = template['total_pervious']

    print('calculating height to width ratio')
    # from masson 2020: h/w = lambda_w/(2*(1-lambda_p))
    ds['height_to_width'] = ds['wall_density']/(2*(1.-ds['building_fraction'])).round(4)

    # replace undefined with nan (where building fraction = 1)
    ds['height_to_width'].values = np.where(ds['height_to_width'].values==np.inf,np.nan,ds['height_to_width'].values)

    print('calculating skyview factor')
    # from masson 2020: [(h/w)^2 + 1]^(1/2) - h/w
    ds['skyview_factor'] = ((ds['height_to_width']**2 + 1.)**(1./2) - ds['height_to_width']).round(4)

    ds['displacement_mac'], ds['roughness_mac'] = calc_site_roughness(
        bldhgt_ave=ds['building_height'],
        sigma_h=ds['building_height_std'],
        lambda_p=ds['building_fraction'],
        lambda_pv=ds['tree_fraction'],
        lambda_f=ds['frontal_density'],
        bldhgt_max=ds['building_height_max'],
        mode=0) # macdonald et al 1998 staggered grid

    ds['displacement_kanda'], ds['roughness_kanda'] = calc_site_roughness(
        bldhgt_ave=ds['building_height'],
        sigma_h=ds['building_height_std'],
        lambda_p=ds['building_fraction'],
        lambda_pv=ds['tree_fraction'],
        lambda_f=ds['frontal_density'],
        bldhgt_max=ds['building_height_max'],
        mode=3) # Kanda et al 2013

    ds['displacement_mac'] = ds['displacement_mac'].round(4)
    ds['roughness_mac'] = ds['roughness_mac'].round(4)
    ds['displacement_mac'] = ds['displacement_mac'].round(4)
    ds['roughness_mac'] = ds['roughness_mac'].round(4)
    ds['displacement_kanda'] = ds['displacement_kanda'].round(4)
    ds['roughness_kanda'] = ds['roughness_kanda'].round(4)

    ################################################################################

    ds = set_attributes(ds,res)

    encoding = {var:{'zlib':True}for var in ds.data_vars}

    for key in ds.keys():
        ds[key].encoding['_FillValue'] = missing_float

    fname = f'{outpath}/{domain.capitalize()}_surface_{res}m_{output_v}'
    print(f'saving output to {fname}')
    ds.to_netcdf(f'{fname}.nc',encoding=encoding)
    
    try:
        ds.rio.write_crs("epsg:4326", inplace=True)
        ds.rio.to_raster(f'{fname}.tiff', compress='LZW')
    except Exception as e:
        print('raster tiff save error:')
        print(e)

    # get all data into dataframe in form for geopandas
    df_all = ds.to_dataframe().sort_values(['longitude','latitude'], 
        ascending=[True,False]).reset_index()

    gf = gpd.GeoDataFrame(df_all, 
            geometry=cell.geometry,crs='EPSG:4283')

    gf.to_file(f'{projpath}/geopackages/g{res}_{domain}.gpkg', driver='GPKG')

    return

def calc_all_grids_from_template(buildings,cell,template,domain='sample'):

    '''calculates grid-level morphology from building data using template grid '''

    print(f'calculating all grids using template')

    # set building morphology as point at centroid of individual buildings
    print('\nNOTE: Warning regarding centroid being likely incorrect can be ignored as the small scales being processed will have negligible error\n')
    points = gpd.GeoDataFrame(buildings, crs = f'EPSG:4283', geometry=gpd.points_from_xy(buildings.geometry.centroid.x, buildings.geometry.centroid.y))
    points['bld_hgt_wgt'] = points['bld_hgt'] * points['area']
    # spatial join building points within grid cells
    building_points = gpd.sjoin(points, cell, how='left', op='within')
    grp = building_points.groupby('index_right')

    gridded = pd.DataFrame()
    gridded['longitude'] = np.meshgrid(template.latitude.values,template.longitude.values)[1].flatten()
    gridded['latitude'] = np.meshgrid(template.latitude.values,template.longitude.values)[0].flatten()

    print('calculating building heights')
    gridded['building_height'] = grp['bld_hgt_wgt'].sum()/grp['area'].sum() # use weighted calculation so larger buildings in grid have greater effect
    gridded['building_height_max'] = grp['bld_hgt'].max()
    gridded['building_height_std'] = grp['bld_hgt'].std()

    gridded['wall_density'] = grp['wall_area'].sum()/cell['cell_area']
    gridded['frontal_density'] = grp['avg_frontal'].sum()/cell['cell_area']

    ## replace nan with zero
    # keys = ['wall_density','frontal_density']
    # gridded[keys] = gridded[keys].fillna(0)

    #################################################################
    # # use building points building_fraction for morphology calculations, RESULTS IN VERY HIGH H/W RATIO
    # gridded['building_fraction'] = grp['area'].sum()/cell['cell_area']
    # gridded['building_fraction'] = np.where(gridded['building_fraction'].values > 0.99,0.99,gridded['building_fraction'].values)

    return gridded

def create_grid_from_template(template,crs='EPSG:4283'):

    '''calculate shapely polygon for each grid cell using netcdf template '''

    cell_size_x = template.longitude.values[1] - template.longitude.values[0]
    cell_size_y = template.latitude.values[1] - template.latitude.values[0]

    # create the cells in a loop
    grid_cells = []
    for x0 in template.longitude.values+cell_size_x/2:
        for y0 in template.latitude.values-cell_size_y/2:
            # bounds
            x1 = x0-cell_size_x
            y1 = y0+cell_size_y
            grid_cells.append( shapely.geometry.box(x0, y0, x1, y1)  )
    cell = gpd.GeoDataFrame(grid_cells, columns=['geometry'],crs=crs)

    cell['cell_area'] = cell.apply(get_dist_x,axis=1)*cell.apply(get_dist_y,axis=1)
                                     
    return cell

def get_midpoint(gdf):

    # total area for the grid
    xmin, ymin, xmax, ymax= gdf.total_bounds
    xmid = (xmin+xmax)/2
    ymid = (ymin+ymax)/2

    return xmid,ymid

def get_perimeter(shp):

    # convert shapely polygon from 3D to 2D (ignore height)
    geom = shapely.wkb.loads(
        shapely.wkb.dumps(shp, output_dimension=2))

    # calculate perimeter along the geodesic Earth
    perim = cgeo.Geodesic().geometry_length(geom)

    return perim

def get_dist_x(row):
    '''building frontal wall area in x direction'''

    bounds = row['geometry'].bounds

    spoint = (bounds[0],bounds[1])
    epoint = (bounds[2],bounds[1])

    x = cgeo.Geodesic().inverse(spoint,epoint).base[0][0]

    return x

def get_dist_y(row):

    bounds = row['geometry'].bounds

    spoint = (bounds[0],bounds[1])
    epoint = (bounds[0],bounds[3])

    y = cgeo.Geodesic().inverse(spoint,epoint).base[0][0]

    return y

def regrid_raster_mean(rorig,ds,res,method='nearest'):

    orig_res = rorig.rio.resolution()[0]
    sfactor = int(res/orig_res)
    new = rorig.coarsen(x=sfactor,y=sfactor,boundary='pad').mean()
    new = new.rio.reproject(dst_crs='+proj=longlat',nodata=np.nan)
    # new = new.where(new>0)
    new = new.rename({'x':'longitude','y':'latitude'}).drop('spatial_ref')

    da = new.interp_like(ds,method=method)

    return da

def regrid_raster_sd(rorig,ds,res,method='nearest'):

    orig_res = rorig.rio.resolution()[0]
    sfactor = int(res/orig_res)
    new = rorig.coarsen(x=sfactor,y=sfactor,boundary='pad').std()
    new = new.rio.reproject(dst_crs='+proj=longlat',nodata=np.nan)
    # new = new.where(new>0)
    new = new.rename({'x':'longitude','y':'latitude'}).drop('spatial_ref')

    da = new.interp_like(ds,method=method)

    return da

def create_ds(df):

    '''convert pandas dataframe to xarray dataset'''

    ds = df.set_index(['latitude', 'longitude']).sort_index().to_xarray()

    # attributes to coordinates
    ds.latitude.attrs = {
            'axis':           'Y',
            'long_name':      'latitude',
            'standard_name':  'latitude',
            'units':          'degrees_north',
        }
    ds.longitude.attrs = {
            'axis':           'X',
            'long_name':      'Longitude',
            'standard_name':  'longitude',
            'units':          'degrees_east',
        }

    return ds

def set_attributes(ds,res):

    '''set dataset attributes '''

    ### surface cover ###
    if 'building_fraction' in ds.keys():
        ds['building_fraction'].attrs = {
                'long_name'   : 'building footprint fraction (lambda_p)',
                'units'       : '-',
                'description' : 'building footprint area as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'tree_fraction' in ds.keys():
        ds['tree_fraction'].attrs = {
                'long_name'   : 'tree canopy fraction (lambda_vt)',
                'units'       : '-',
                'description' : 'tree canopy plan area as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'lowveg_fraction' in ds.keys():
        ds['lowveg_fraction'].attrs = {
                'long_name'   : 'grass and low vegetation fraction (lambda_vl)',
                'units'       : '-',
                'description' : 'low vegetation (grass, shrubs, other vegetation) as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'water_fraction' in ds.keys():
        ds['water_fraction'].attrs = {
                'long_name'   : 'open water fraction (lambda_wa)',
                'units'       : '-',
                'description' : 'all open water (ocean, lakes, pools) as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'bareearth_fraction' in ds.keys():
        ds['bareearth_fraction'].attrs = {
                'long_name'   : 'bare earth, sand and rock fraction (lambda_be)',
                'units'       : '-',
                'description' : 'bare earth including construction sites, rock, sand and sparsely vegetated areas as fraction of grid area, corrected for cloud and shadow fractions'
            }
    if 'roadpath_fraction' in ds.keys():
        ds['roadpath_fraction'].attrs = {
                'long_name'   : 'road, path and other hard surfaces on ground (lambda_if)',
                'units'       : '-',
                'description' : 'all hard surfaces on ground excluding buildings, defined as "impervious surface fraction" in Stewart and Oke, 2012, corrected for cloud and shadow fractions'
            }
    if 'total_built' in ds.keys():
        ds['total_built'].attrs = {
                'long_name'   : 'total built fraction (lambda_tb)',
                'units'       : '-',
                'description' : 'all impervious surfaces including buildings, roads, paths and other hard surfaces, corrected for cloud and shadow fractions'
            }
    if 'total_pervious' in ds.keys():
        ds['total_pervious'].attrs = {
                'long_name'   : 'total pervious fraction (lambda_tp)',
                'units'       : '-',
                'description' : 'all pervious surfaces including vegetation, water and bare earth, corrected for cloud and shadow fractions'
            }

    # morphology
    if 'frontal_density' in ds.keys():
        ds['frontal_density'].attrs = {
                'long_name'   : 'frontal area density (lambda_f)',
                'units'       : '-',
                'description' : 'sum of cardinally averaged building frontal area as fraction of grid area'
            }
    if 'wall_density' in ds.keys():
        ds['wall_density'].attrs = {
                'long_name'   : 'wall area density (lambda_w)',
                'units'       : '-',
                'description' : 'sum of building wall area as fraction of grid area'
            }

    if 'building_height' in ds.keys():
        ds['building_height'].attrs = {
                'long_name'   : 'building height mean (H_avg)',
                'units'       : 'm',
                'description' : 'mean building height in grid cell (avg. of geoscape roof and eave height)'
            }
    if 'building_height_max' in ds.keys():
        ds['building_height_max'].attrs = {
                'long_name'   : 'building height maximum (H_max)',
                'units'       : 'm',
                'description' : 'maximum building height in grid cell (avg. of geoscape roof and eave height)'
            }
    if 'building_height_std' in ds.keys():
        ds['building_height_std'].attrs = {
                'long_name'   : 'building height standard deviation (H_std)',
                'units'       : 'm',
                'description' : 'standard deviation of building height in grid cell (avg. of geoscape roof and eave height)'
            }

    if 'tree_height' in ds.keys():
        ds['tree_height'].attrs = {
                'long_name'   : 'tree height average (HT_ave)',
                'units'       : 'm',
                'description' : 'average vegetation canopy height in grid'
            }
    if 'tree_height_std' in ds.keys():
        ds['tree_height_std'].attrs = {
                'long_name'   : 'tree height standard deviation (HT_ave)',
                'units'       : 'm',
                'description' : 'standard deviation in vegetation canopy height in grid'
            }

    if 'roughness_mac' in ds.keys():
        ds['roughness_mac'].attrs = {
                'long_name'   : 'aerodynamic roughness length (Z_0,mac)',
                'units'       : 'm',
                'description' : 'roughness length for staggered arrays from Macdonald et al., 1998: https://doi.org/10.1016/S1352-2310(97)00403-2'
            }
    if 'displacement_mac' in ds.keys():
        ds['displacement_mac'].attrs = {
                'long_name'   : 'displacement height (Z_d,mac)',
                'units'       : 'm',
                'description' : 'zero-plane displacement height from Macdonald et al., 1998: https://doi.org/10.1016/S1352-2310(97)00403-2'
            }
    if 'roughness_kanda' in ds.keys():
        ds['roughness_kanda'].attrs = {
                'long_name'   : 'aerodynamic roughness length (Z_0,kan)',
                'units'       : 'm',
                'description' : 'roughness length for staggered arrays from Kanda et al., 2013: https://doi.org/10.1007/s10546-013-9818-x'
            }
    if 'displacement_kanda' in ds.keys():
        ds['displacement_kanda'].attrs = {
                'long_name'   : 'displacement height (Z_d,kan)',
                'units'       : 'm',
                'description' : 'zero-plane displacement height from Kanda et al., 2013: https://doi.org/10.1007/s10546-013-9818-x'
            }

    ### canyon attributes ###
    if 'height_to_width' in ds.keys():
        ds['height_to_width'].attrs = {
                'long_name'   : 'canyon height to width aspect ratio (h/w)',
                'units'       : '-',
                'description' : 'average aspect ratio assuming street canyon geometry using Eq 1 of Masson et al. 2020: https://doi.org/10.1016/j.uclim.2019.100536'
            }
    if 'skyview_factor' in ds.keys():
        ds['skyview_factor'].attrs = {
                'long_name'   : 'skyview factor (psi)',
                'units'       : '-',
                'description' : 'average skyview factor assuming street canyon geometry using Eq 2 of Masson et al. 2020: https://doi.org/10.1016/j.uclim.2019.100536'
            }

    # other
    if 'cell_area' in ds.keys():
        ds['cell_area'].attrs = {
                'long_name'   : 'cell area',
                'units'       : 'm2',
                'description' : f'area of grid cell, nominally at {res} m resolution'
            }
    if 'bld_hgt' in ds.keys():
        ds['bld_hgt'].attrs = {
                'long_name'   : 'building height average',
                'units'       : 'm',
                'description' : 'average for midpoint of roof and eave height in grid (no buiding area weighting)'
            }
    if 'bld_hgt_wgt' in ds.keys():
        ds['bld_hgt_wgt'].attrs = {
                'long_name'   : 'building height weighted average',
                'units'       : 'm',
                'description' : 'area weighted average for midpoint of roof and eave height in grid'
            }
    if 'roof_hgt' in ds.keys():
        ds['roof_hgt'].attrs = {
                'long_name'   : 'roof height average',
                'units'       : 'm',
                'description' : 'average roof height in grid (no buiding area weighting)'
            }
    if 'roof_hgt_wgt' in ds.keys():
        ds['roof_hgt_wgt'].attrs = {
                'long_name'   : 'roof height weighted average',
                'units'       : 'm',
                'description' : 'area weighted average roof height in grid'
            }
    if 'eave_hgt' in ds.keys():
        ds['eave_hgt'].attrs = {
                'long_name'   : 'eave height average',
                'units'       : 'm',
                'description' : 'average eave height in grid (no buiding area weighting)'
            }
    if 'eave_hgt_wgt' in ds.keys():
        ds['eave_hgt_wgt'].attrs = {
                'long_name'   : 'eave height weighted average',
                'units'       : 'm',
                'description' : 'area weighted average eave height in grid'
            }

    ds.attrs['title']       = f'Urban characteristics for {domain.capitalize()} at {res} m resolution derived from Geoscape datasets'
    ds.attrs['version']     = output_v
    ds.attrs['institution'] = 'ARC Centre of Excellence for Climate Extremes, UNSW Sydney, Australia'
    ds.attrs['source']      = 'Derived from Geoscape datasets, including Buildings v2.0, Trees v1.6 and Surface cover v1.6: PSMA Australia, 2020. https://geoscape.com.au/'
    ds.attrs['reference']  = 'A transformation in city-descriptive input data for urban climate models'
    ds.attrs['licence']  = 'The contents of this file are restricted, not to be used or distributed without permission from the authors.'

    return ds

def calc_site_roughness(bldhgt_ave,sigma_h,lambda_p,lambda_pv,lambda_f,bldhgt_max=None,mode=0):
    '''estimate urban site roughness length for momentum and zero-plane displacement based on various methods 
    as described in Grimmond and Oke (1999): Aerodynamic Properties of Urban Areas Derived from Analysis of Surface Form .

    modes: four different morphometric methods to calculate roughness and displacement:
        - 0: Macdonald et al. 1998: An improved method for the estimation of surface roughness of obstacle arrays
        - 1: Kent et al 2017: Aerodynamic roughness parameters in cities: Inclusion of vegetation
        - 2: Millward-Hopkins et al., (2011) per Kent et al. 2017 eq 11 and 12 # TYPO IN EQ 12 OF KENT PAPER
        - 3: Kanda et al., 2013

    Inputs
    ------
    bldhgt_ave [m] : building mean height
    sigma_h [m]    : building height standard deviation
    lambda_p [0-1] : building plan area fraction
    lambda_pv [0-1]: tree plan area fraction
    lambda_f [1]   : wall frontal area fraction
    bldhgt_max [m] : maximum building height (if none assume 1.25 times SD)

    Returns
    -------
    zd [m] zero-plane displacement
    z0 [m] rougness length for momentum

    Warning
    -------
    in mode=0 wall frontal area (assuming random canyon orientation) calculated from Porsen et al 2010 (DOI:10.1002/qj.668)
    in mode=1 (including vegetation), vegetation frontal area index assumed to be equal to tree plan fraction (from Table 1 in Kent et al 2017)
    in mode=3 (Kanda) maximum building height assumed to be 2 standard deviations
    '''

    vonkarm = 0.4
    dragbld = 1.2

    if mode==0:
        print('zd,z0 from Macdonald et al. 1998')

        alpha   = 4.43 # for staggered arrays
        beta    = 1.0  # for staggered arrays

        zd = (1. + alpha**(-lambda_p)*(lambda_p - 1.))*bldhgt_ave    # eq 23 from Mac1998
        z0 = ((1. - zd/bldhgt_ave) * np.exp( -1.*(0.5*beta*(dragbld/vonkarm**2)*(1.-zd/bldhgt_ave)*lambda_f)**(-0.5) ))*bldhgt_ave  # eq 26 from Mac1998

    if mode == 1: # with veg from Kent et al (2018) https://doi.org/10.1007/s11252-017-0710-1
        print('zd,z0 from Macdonald with vegetation from Kent et al. 2017 and 2018')
        # Aerodynamic roughness variation with vegetation: analysis in a suburban neighbourhood and a city park

        alpha   = 4.43 # for staggered arrays
        beta    = 1.0  # for staggered arrays

        lambda_fv = lambda_pv # estimated equality based on Table 1 of Kent et al 2017 for different sites
        P_3D = 0.4   # leaf-on porosity

        lambda_tot = (lambda_p + lambda_pv*(1. - P_3D)) # eq 4 from Kent et al (2018)
        Pv  = (-1.251*P_3D**2 + 0.489*P_3D + 0.803)/dragbld  # eq 7 from Kent et al (2018)
        weighted_frontal_area = lambda_f + lambda_fv*Pv

        zd = (1. + alpha**(-lambda_tot)*(lambda_tot-1.))*bldhgt_ave # eq 5
        z0 = (1. - zd/bldhgt_ave)*np.exp( -1.*(0.5*beta*(dragbld/vonkarm**2)*(1.-zd/bldhgt_ave)*weighted_frontal_area)**(-0.5) )*bldhgt_ave # eq 6

    if mode == 2: # Millward-Hopkins et al. (2011), per Kent et al. 2017 eq 11 and 12
        print('zd,z0 from Millward-Hopkins et al. 2011, per Kent et al. 2017 eq 11 and 12')

        MhoUzd_on_Hav_A = (19.2*lambda_p - 1. + np.exp(-19.2*lambda_p))/(19.2*lambda_p*(1.- np.exp(-19.2*lambda_p)))
        MhoUzd_on_Hav_B = (117*lambda_p + (187.2*lambda_p**3 - 6.1)*(1.- np.exp(-19.2*lambda_p)) )/( (1.+114*lambda_p + 187*lambda_p**3)*(1.-np.exp(-19.2*lambda_p)) )
        MhoUzd_on_Hav = np.where( lambda_p >= 0.19, MhoUzd_on_Hav_A,  MhoUzd_on_Hav_B )

        Mho_exp = np.exp( -1.*(0.5*dragbld*vonkarm**(-2)*lambda_f)**(-0.5) )
        MhoUz0_on_Hav = ( (1. - MhoUzd_on_Hav)*Mho_exp)

        zd = bldhgt_ave*(MhoUzd_on_Hav + ((0.2375 * np.log(lambda_p) + 1.1738)*(sigma_h/bldhgt_ave)) )

        # z0 = bldhgt_ave*(MhoUz0_on_Hav + (np.exp(0.8867*lambda_f) - 1.)*(sigma_h/bldhgt_ave)**(np.exp(2.3271*lambda_f))) # ERROR IN KENT ET AL PAPER
        z0 = bldhgt_ave*(MhoUz0_on_Hav + np.exp(0.8867*lambda_f - 1.)*(sigma_h/bldhgt_ave)**(np.exp(2.3271*lambda_f))) # correction based on UMEP

    if mode == 3:
        print('zd,z0 Kanda et al 2013')

        alpha   = 4.43 # for staggered arrays
        beta    = 1.0  # for staggered arrays

        if np.all(bldhgt_max)==None:
            bldhgt_max = 1.5*sigma_h + bldhgt_ave # assume bldhgt_max
            # bldhgt_max = 12.51*sigma_h**0.77 # eq 3 from Kanda et al 2013

        # first calculate macdonald et al 1998 (for later scaling)
        mac_zd = (1. + (lambda_p - 1.)*alpha**(-lambda_p))*bldhgt_ave
        mac_z0 = ((1. - mac_zd/bldhgt_ave) * np.exp( -1.*(0.5*beta*(dragbld/vonkarm**2)*(1.-mac_zd/bldhgt_ave)*lambda_f)**(-0.5) ))*bldhgt_ave

        # then Kanda et al 2013 scaling parameters
        a0,b0,c0 = 1.29, 0.36, -0.17
        a1,b1,c1 = 0.71, 20.21, -0.77

        X = (sigma_h + bldhgt_ave)/bldhgt_max  # eq. 10b
        Y = (lambda_p*sigma_h)/bldhgt_ave        # eq. 12b

        zd = bldhgt_max*(c0*X**2 + (a0*lambda_p**b0 - c0)*X) # eq 10a 
        z0 = mac_z0*(a1 + b1*Y**2 + c1*Y) # eq 12a

    zd = np.round(zd,4)
    zd = xr.where(zd<0,0,zd)
    z0 = np.round(z0,4)
    z0 = xr.where(z0<0,0,z0)

    return zd,z0

################################################################################

def main_plot_sample():

    buildings = gpd.read_file(f'{projpath}/geopackages/buildings_{domain}.gpkg')

    gf1,ds1 = load_files(res_list[0])
    gf2,ds2 = load_files(res_list[1])
    gf3,ds3 = load_files(res_list[2])

    key1='bld_hgt'
    keys=['building_height','building_height_max','building_height_std','tree_height',
            'wall_density','building_fraction','roadpath_fraction','total_pervious','total_built',
            'height_to_width','skyview_factor','roughness_mac','displacement_mac']

    for key2 in keys:
        print(f'plotting {key2}')
        plot_geoscape_and_grid(buildings,gf1,gf2,gf3,key1,key2)

    return

def load_files(res):

    ''' read building and derived data for plotting at particular resolution'''

    gf = gpd.read_file(f'{projpath}/geopackages/g{res}_{domain}.gpkg')
    ds = xr.open_dataset(f'{outpath}/{domain.capitalize()}_surface_{res}m_{output_v}.nc')

    return gf,ds

def plot_geoscape_and_grid(gdf,gf1,gf2,gf3,key1,key2):

    '''creat plot comparing original and derived (lower resolution) data'''
    
    print(f'plotting {key1} and {key2}')

    plt.close('all')
    fig, axes = plt.subplots(nrows=2,ncols=2,figsize=(9,6),sharey=True,sharex=True)

    ax = axes.flatten()

    plot_orig(gdf, ax[0],key1,cmap='plasma')
    plot_grid(gf1, ax[1],key2,res=f'{res_list[0]}m',cmap='plasma')
    plot_grid(gf2, ax[2],key2,res=f'{res_list[1]}m',cmap='plasma')
    plot_grid(gf3, ax[3],key2,res=f'{res_list[2]}m',cmap='plasma')

    fig.subplots_adjust(left=0.05 ,right=0.92,bottom=0.05,top=0.95,wspace=0.1,hspace=0.1)
    fig.savefig(f'{plotpath}/res_{key1}_{key2}_{output_v}.png', dpi=400, bbox_inches='tight')
    # plt.show()

    return

def plot_orig(gdf,ax,key,cmap='plasma'):
    ''' plot original geopandas with colorbar (no edgecolor)'''

    un,tt,vmax = get_labels(key)


    ax.set_title(f'Original data: {tt}', fontsize=11)
    div = make_axes_locatable(ax)
    cax = make_axes_locatable(ax).append_axes('right', size="3%", pad=-0.15)
    cax.text(1.10,1.06,un,fontsize=8,va='bottom',ha='left',transform=cax.transAxes)

    gdf.plot(column=key, ax=ax, cax=cax, vmax=vmax, legend=True, legend_kwds={'label':'', 'extend':'max'}, cmap=cmap)
    ax.tick_params(axis='both', labelsize=9)
    ax.locator_params(axis='x', nbins=5)
    ax.locator_params(axis='y', nbins=6)

    ax.text(0.01,0.99,f'min: {gdf[key].min():.2f}', fontsize=6, ha='left', va='top', transform=ax.transAxes)
    ax.text(0.01,0.95,f'max: {gdf[key].max():.2f}', fontsize=6, ha='left', va='top', transform=ax.transAxes)

    return ax

def plot_grid(gdf,ax,key,res,cmap='plasma'):
    ''' plot gridded geopandas with colorbar (with edgecolor)'''

    un,tt,vmax = get_labels(key)

    ax.set_title(f'Derived: {res} {tt}', fontsize=11)
    div = make_axes_locatable(ax)
    cax = make_axes_locatable(ax).append_axes('right', size="3%", pad=-0.15)
    cax.text(1.10,1.06,un,fontsize=8,va='bottom',ha='left',transform=cax.transAxes)

    gdf.plot(column=key, ax=ax, cax=cax, vmax=vmax, cmap=cmap, 
        legend=True, legend_kwds={'label': '', 'extend':'max'}, edgecolor='0.8',linewidth=0.3)
    ax.tick_params(axis='both', labelsize=9)
    ax.locator_params(axis='x', nbins=5)
    ax.locator_params(axis='y', nbins=6)

    ax.text(0.01,0.99,f'min: {gdf[key].min():.2f}', fontsize=6, ha='left', va='top', transform=ax.transAxes)
    ax.text(0.01,0.95,f'max: {gdf[key].max():.2f}', fontsize=6, ha='left', va='top', transform=ax.transAxes)

    return ax

def get_labels(key):

    '''get key labels for plotting '''

    units=''
    title = ''
    vmax=None

    if key=='area':
        title = 'building footprint area'
        units = '[$m^2$]'
    elif key=='geometry':
        title = 'building footprint'
        units = '[-]'
    elif key == 'bld_hgt':
        title = 'building height'
        units = '[m]'
        vmax = 100
    elif key == 'building_height':
        title = 'mean building height'
        units = '[m]'
        vmax = 30
    elif key == 'building_height_max':
        title = 'maximum building height'
        units = '[m]'
        vmax = 100
    elif key == 'building_height_std':
        title = 'standard deviation of building height'
        units = '[m]'
        vmax = 10
    elif key == 'tree_height':
        title = 'mean tree height'
        units = '[m]'
        vmax = 30
    elif key == 'wall_density':
        title = 'wall area density'
        units = '[-]'
        vmax = 2
    elif key == 'frontal_density':
        title = 'frontal area density'
        units = '[-]'
        vmax = 1
    elif key == 'building_fraction':
        title = 'building footprint fraction'
        units = '[-]'
        vmax = 1
    elif key == 'roadpath_fraction':
        title = 'roads, paths, hard surfaces fraction'
        units = '[-]'
        vmax = 1
    elif key == 'tree_fraction':
        title = 'tree fraction'
        units = '[-]'
        vmax = 1
    elif key == 'total_pervious':
        title = 'total pervious fraction'
        units = '[-]'
        vmax = 1
    elif key == 'total_vegetation':
        title = 'total vegetation fraction'
        units = '[-]'
        vmax = 1
    elif key == 'total_built':
        title = 'total built fraction'
        units = '[-]'
        vmax = 1
    elif key == 'skyview_factor':
        title = 'skyview factor'
        units = '[-]'
        vmax = 1
    elif key == 'height_to_width':
        title = 'canyon aspect ratio'
        units = '[-]'
        vmax = 2.5
    elif key == 'roughness_mac':
        title = 'roughness length (Mac 1998)'
        units = '[m]'
        vmax = 5
    elif key == 'displacement_mac':
        title = 'displacement height (Mac 1998)'
        units = '[m]'
        vmax = 10
    else:
        print('WARNING: key not recognised')

    return units,title,vmax

################################################################################

if __name__ == "__main__":

    ################################################################################
    # path setup

    for mkpath in [plotpath,outpath,f'{projpath}/geopackages']:
        if not os.path.exists(mkpath):
            os.makedirs(mkpath)
            print("created folder : ", mkpath)

    if domain == 'sample':
        tif_fpaths = [f'{datapath}/Geoscape_Sample_data_2021/06_Surface_Cover_Trees/NSW_SURFACECOVER_2M_Z56_16837.tif']
        tree_fname = f'{datapath}/Geoscape_Sample_data_2021/06_Surface_Cover_Trees/NSW_TREES_Z56_16837.tif'
        shp_fname =  f'{datapath}/Geoscape_Sample_data_2021/01_Buildings/SHP/buildings.shp'

    ################################################################################

    main()


Associated with the manuscript: A transformation in city-descriptive input data for urban climate models Manuscript authors:  <NAME>, <NAME>,  <NAME>,  <NAME>, <NAME>  Copyright 2022 <NAME>  Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at      http://www.apache.org/licenses/LICENSE-2.0  Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.  Instructions:  1. Get sample data from Geoscape: https://geoscape.com.au/get-sample/ 2. Install necessary python dependancies 3. Update project and data paths as required 4. Run script ################################################################################ # user input # update path to working folder # update path to geoscape data # netcdf/tif outputs # figures output # list of resolutions to process in m # if true, splits buildings on grid edges (long processing) # if true, plots comparison between different resolutions (for sample) ################################################################################ main function for resampling geoscape surface cover data to a lower resolution and recategorise # read tif # get band 1 # remove fill value ################################################################################ ################################################################################ # geoscape.rio.write_crs("epsg:4326", inplace=True) # geoscape.rio.to_raster(f'{fname}.tiff', compress='LZW') ################################################################################ # derived.rio.write_crs("epsg:4326", inplace=True) # derived.rio.to_raster(f'{fname}.tiff', compress='LZW') ################################################################################ function to take categorical band and calculate its fraction at lower resolution # get resolution factors # get all key values, convert to 1 # calculate sums over coarsened grid # get all non-nan values, convert to 1 # calculate sums over coarsened grid create netcdf attributes from geoscape surface cover data      # Table 2: SURFACE_COVER_2M_TYPE_AUT Codes from 'Surface Cover 1.6 Product Description.pdf'      2 Bare Earth                Includes sand dunes, desert, rock outcrops, bare soil other than bare agricultural land, and sparsely vegetated areas of grass and shrub. Non-vegetated strip mines and quarries except where covered by development or water.     3 Road and Path             Roads and parking lots covered in a man-made material excluding hard packed dirt trails.     4 Grass                     Grass and herbaceous areas. The category may include herbaceous wetlands if images are collected during dry season or periods of drought.     5 Trees                     All trees including deciduous and evergreen woody vegetation.     6 Unspecified Vegetation    Any other vegetative material not included within the Grass or Tree class. This may include, but is not limited to, shrub, sc rub, agriculture, and aquatic plants.     7 Built-up Areas            Any areas of man-made environments and infrastructure excluding road and paths and buildings.     8 Water                     Depending on the resolution quality of the imagery used, natural water will include streams, canals, ponds, lakes, reservoirs, estuaries and bays.     9 Buildings                 Where the majority of a pixel intersects a Building, vector building polygon representation.      10 Cloud                    The area covered with cloud on Date of collection.     11 Shadow                   The area covered with shadow on Date/time of collection.     12 Swimming Pool            An area identified as a swimming pool. create derived classes and attributes (normalising for cloud and shadow fractions) test fractions sum to 1 # test geoscape classess sum to 1 # test derived classes sum to 1 # test total classes sum to 1 # test total classes sum to 1 ################################################################################ main function for calculating grid-level morphology from geoscape building and tree data     split_by_grid=False: building morphology information placed at centroid point of building for faster processing     split_by_grid=True splits buildings at grid edges so morphology data sits in multiple cells ################################################################################ # format individual building information # drop buildings without height info # split buildings where intersecting with grid (warning, very long processing) # split buildings into cells by clipping ################################################################################ # create gridded data from buildings information ################################################################################ # from masson 2020: h/w = lambda_w/(2*(1-lambda_p)) # replace undefined with nan (where building fraction = 1) # from masson 2020: [(h/w)^2 + 1]^(1/2) - h/w # macdonald et al 1998 staggered grid # Kanda et al 2013 ################################################################################ # get all data into dataframe in form for geopandas calculates grid-level morphology from building data using template grid # set building morphology as point at centroid of individual buildings # spatial join building points within grid cells # use weighted calculation so larger buildings in grid have greater effect ## replace nan with zero # keys = ['wall_density','frontal_density'] # gridded[keys] = gridded[keys].fillna(0) ################################################################# # # use building points building_fraction for morphology calculations, RESULTS IN VERY HIGH H/W RATIO # gridded['building_fraction'] = grp['area'].sum()/cell['cell_area'] # gridded['building_fraction'] = np.where(gridded['building_fraction'].values > 0.99,0.99,gridded['building_fraction'].values) calculate shapely polygon for each grid cell using netcdf template # create the cells in a loop # bounds # total area for the grid # convert shapely polygon from 3D to 2D (ignore height) # calculate perimeter along the geodesic Earth building frontal wall area in x direction # new = new.where(new>0) # new = new.where(new>0) convert pandas dataframe to xarray dataset # attributes to coordinates set dataset attributes ### surface cover ### # morphology ### canyon attributes ### # other estimate urban site roughness length for momentum and zero-plane displacement based on various methods      as described in Grimmond and Oke (1999): Aerodynamic Properties of Urban Areas Derived from Analysis of Surface Form .      modes: four different morphometric methods to calculate roughness and displacement:         - 0: Macdonald et al. 1998: An improved method for the estimation of surface roughness of obstacle arrays         - 1: Kent et al 2017: Aerodynamic roughness parameters in cities: Inclusion of vegetation         - 2: Millward-Hopkins et al., (2011) per Kent et al. 2017 eq 11 and 12 # TYPO IN EQ 12 OF KENT PAPER         - 3: Kanda et al., 2013      Inputs     ------     bldhgt_ave [m] : building mean height     sigma_h [m]    : building height standard deviation     lambda_p [0-1] : building plan area fraction     lambda_pv [0-1]: tree plan area fraction     lambda_f [1]   : wall frontal area fraction     bldhgt_max [m] : maximum building height (if none assume 1.25 times SD)      Returns     -------     zd [m] zero-plane displacement     z0 [m] rougness length for momentum      Warning     -------     in mode=0 wall frontal area (assuming random canyon orientation) calculated from Porsen et al 2010 (DOI:10.1002/qj.668)     in mode=1 (including vegetation), vegetation frontal area index assumed to be equal to tree plan fraction (from Table 1 in Kent et al 2017)     in mode=3 (Kanda) maximum building height assumed to be 2 standard deviations # for staggered arrays # for staggered arrays # eq 23 from Mac1998 # eq 26 from Mac1998 # with veg from Kent et al (2018) https://doi.org/10.1007/s11252-017-0710-1 # Aerodynamic roughness variation with vegetation: analysis in a suburban neighbourhood and a city park # for staggered arrays # for staggered arrays # estimated equality based on Table 1 of Kent et al 2017 for different sites # leaf-on porosity # eq 4 from Kent et al (2018) # eq 7 from Kent et al (2018) # eq 5 # eq 6 # Millward-Hopkins et al. (2011), per Kent et al. 2017 eq 11 and 12 # z0 = bldhgt_ave*(MhoUz0_on_Hav + (np.exp(0.8867*lambda_f) - 1.)*(sigma_h/bldhgt_ave)**(np.exp(2.3271*lambda_f))) # ERROR IN KENT ET AL PAPER # correction based on UMEP # for staggered arrays # for staggered arrays # assume bldhgt_max # bldhgt_max = 12.51*sigma_h**0.77 # eq 3 from Kanda et al 2013 # first calculate macdonald et al 1998 (for later scaling) # then Kanda et al 2013 scaling parameters # eq. 10b # eq. 12b # eq 10a # eq 12a ################################################################################ read building and derived data for plotting at particular resolution creat plot comparing original and derived (lower resolution) data # plt.show() plot original geopandas with colorbar (no edgecolor) plot gridded geopandas with colorbar (with edgecolor) get key labels for plotting ################################################################################ ################################################################################ # path setup ################################################################################