<reponame>Alpine-Land-Information-Services/PyVF<gh_stars>0
import numpy as np
import pandas as pd
import pyvf

import pkg_resources
species_references=pd.read_csv(pkg_resources.resource_stream(__name__,'data/Species_reference.csv'))
species_references_vf=pd.read_csv(pkg_resources.resource_stream(__name__,'data/Species_reference_vf.csv'))
deciduous_sp_dict = {'BO': 'BOS', 'BM': 'BMS', 'LO': 'LOS', 'QA': 'QAS', 'RA': 'RAS'}

def __get_models_dict():
    import io
    f=pkg_resources.resource_stream(__name__,'data/Tree_models.csv')
    """Feed in models csv with species code, dbh breakpoint, and model name. Return dictionary of species(dbh breakpoint,model name)"""
    models_dict = dict()


    record = f.readline()  # skip header
    while True:
        record = f.readline().decode('UTF-8').strip()
        if len(record) < 1:
            break

        record = str(record).split(',')
        species = record[0]
        breakpoint = record[1]

        try:
            models_dict[species].append((float(record[1]), record[2]))
        except:
            models_dict[species] = list()
            models_dict[species].append((float(record[1]), record[2]))

    return models_dict

models_dict=__get_models_dict()


def convert_sp_deciduous(vf_species):
    vf_species = vf_species.upper()

    if vf_species in deciduous_sp_dict.keys():
        return deciduous_sp_dict[vf_species]
    else:
        return vf_species


def get_model(vf_species, dbh):
    vf_species = convert_sp_deciduous(vf_species)

    for record in models_dict[vf_species]:
        if dbh < record[0]:
            return record[1]+'_A'

    return record[1]+'_A'




def CreateSqGrid(qty_per_ac, variance=0, bbox=[29.2,29.2,158.8,158.8]):
    """Return numpy array of x,y coordinates organized as a square grid with random variance

    Parameters
    qty_per_ac (flt): Point density as number of points per acre
    variance (flt): Set to 0 for a perfect grid, set >1 for very random
    bounds (list): [xmin, ymin, xmax, ymax]

    """
    x_meters = bbox[2] - bbox[0]
    y_meters = bbox[3] - bbox[1]

    # Figure out grid spacing
    acres = ((x_meters * y_meters) * 3.2808399 ** 2) / 43560
    qty = qty_per_ac * acres
    area_per_object = (acres / qty) * 43560 / 3.2808399 ** 2  ## output square meters
    spacing = area_per_object ** .5  ## output meters

    # Figure out number of grid points per axis
    points_x = int(x_meters / spacing)
    points_y = int(y_meters / spacing)

    # Define storage variables
    points = []

    xvals = np.linspace(start=bbox[0], stop=bbox[2], num=(round(x_meters / spacing) + 1))
    yvals = np.linspace(start=bbox[1], stop=bbox[3], num=(round(y_meters / spacing) + 1))
    points = np.array(np.meshgrid(xvals, yvals)).T.reshape(-1, 2)


    # Create an adjustments array with the same shape as the coordinates array. This array represents a random adjustment in meters based on variance as a percent of spacing.

    if variance > .01:
        rand_adj = variance * spacing
        rands = np.random.randint(rand_adj * -1000, rand_adj * 1000, (len(points), 2))
        rands = rands / 1000

        # Apply adjustments
        points = points + rands

    return (points)


def CreateHexGrid(qty_per_ac, variance=0, bbox=[29.2, 29.2, 158.8, 158.8]):
    """Return numpy array of x,y coordinates organized as a hexagonal grid with random variance

    Parameters
    qty_per_ac (flt): Point density as number of points per acre
    variance (flt): Set to 0 for a perfect grid, set >1 for very random
    bounds (list): [xmin, ymin, xmax, ymax]


    """
    x_meters = bbox[2] - bbox[0]
    y_meters = bbox[3] - bbox[1]

    # Figure out hex grid spacing
    acres = ((x_meters * y_meters) * 3.2808399 ** 2) / 43560
    qty = qty_per_ac * acres
    area_per_object = (acres / qty) * 43560 / 3.2808399 ** 2  # Output square meters

    side_radius = (
                              area_per_object / 2.598076) ** .5  # Constant is equivalent to 1.5*(3**.5). Radius, side length, and "size" are all the same thing.
    apothem = side_radius * 0.8660254037844386  # Constant is equivalent to 3**.5/2

    x_spacing = 2 * apothem
    y_spacing = 1.5 * side_radius

    # Figure out number of grid points per axis
    points_x = round(x_meters / x_spacing)
    points_y = round(y_meters / y_spacing)

    yvals = np.linspace(bbox[1], bbox[3], points_y)
    xvals = np.linspace(bbox[0], bbox[2], points_x)

    points = np.array(np.meshgrid(xvals, yvals)).T.reshape(-1, 2)

    filt = np.isin(points[:, 1], yvals[::2])
    points[filt, 0] = points[filt, 0] + apothem

    # Create an adjustments array with the same shape as the coordinates array. This array represents a random adjustment in meters based on variance as a percent of spacing.

    if variance >= .01:
        rand_adj = variance * y_spacing
        rands = np.random.randint(rand_adj * -1000, rand_adj * 1000, (len(points), 2))
        rands = rands / 1000

        # Apply adjustments
        points = points + rands

    return (points)


def CreateClumpsGrid(qty_per_ac, variance=0, bbox=[29.2,29.2,158.8,158.8]):
    """Return numpy array of x,y coordinates organized as a square grid with random variance

    Arguments
    qty_per_ac: Point density as number of points per acre
    variance: Set to 0 for a perfect grid, set >1 for very random
    bounds: [xmin, ymin, xmax, ymax]


    """
    variance = variance + .3

    clump_tpa = .8 * qty_per_ac * (1 + variance) ** 1.3
    grid1_tpa = .8 * qty_per_ac / (1 + variance) ** .8
    clump_percent = .1
    clump_size = 5

    grid1 = CreateHexGrid(grid1_tpa, variance, bbox)

    np.random.shuffle(grid1)

    clump_grid = grid1

    grid_slice = int(clump_percent * grid1.shape[0])

    grid1 = grid1[0:grid_slice, :]

    for point in grid1:
        win = 0
        while win == 0 and clump_size < 100:
            try:

                clump_radius1 = clump_size * np.random.randint(6, 14) / 10
                clump_radius2 = clump_size * np.random.randint(6, 14) / 10

                clump_bbox = [point[0] - clump_radius1, point[1] - clump_radius2, point[0] + clump_radius1,
                              point[1] + clump_radius2]
                clump = CreateHexGrid(clump_tpa, variance, clump_bbox)
                win = 1

            except:
                clump_size = clump_size + 1

        clump = CreateHexGrid(clump_tpa, variance, clump_bbox)
        clump_grid = np.append(clump_grid, clump, 0)
        clump_grid = np.append(clump_grid, clump, 0)

    return (clump_grid)

def treelist_from_stand_table(path,skiprows=0,custom_species=True):
    import pandas as pd
    try:
        df=pd.read_csv(path,index_col=0,header=None,skiprows=skiprows)
    except:
        df=pd.DataFrame(path)
        df=pd.concat([pd.DataFrame([path.columns.values], columns=path.columns.values), path])

    df.drop(index=np.nan,inplace=True,errors='ignore')

    # Import and clean stand table. Remove totals column and totals row. Use DBH column as index. Fill NA with 0.
    df.fillna(0, inplace=True)

    for index in df.index:
        if 'total' in str(index).lower():
            df.drop(index=index, inplace=True, errors='ignore')

    for index in df.index:
        for column in df.columns:
            if 'total' in str(df[column][index]).lower():
                df.drop(columns=column, inplace=True, errors='ignore')


    trees=np.array([[],[],[]])

    if custom_species == True:
        i_vf_species=2
        index = 2

    else:
        i_vf_species=0
        index = 0

    while index+1 < df.index.shape[0]:
        index=index+1
        for column in range(0,df.columns.shape[0]):
            vf_species=df.iloc[i_vf_species,column]
            dbh=df.index[index]
            tpa=df.iloc[index,column]
            tree=np.array([[tpa],[vf_species],[dbh]])
            trees=np.append(trees,tree,1)

    trees = trees.transpose()

    from pyvf.treelist import Treelist
    out=Treelist(trees[:,0],trees[:,1],trees[:,2])
    return out



def convert_species(sp_id, to_convention,from_convention):
    """
    Arguments
    sp_id: a single id code to convert
    to_convention: ID convention to convert to. Valid options: 'FIA', 'PLANTS', or 'VF'
    """

    if from_convention=='VF':

        codes=species_references_vf[species_references_vf['VF_SPCD']==sp_id.upper()]

        if to_convention == 'FIA':
            code=str(codes['SPCD'].iloc[0])
        elif to_convention=='PLANTS':
            code=str(codes['SPECIES_SYMBOL'].iloc[0])

    elif from_convention=='PLANTS':
        codes=species_references[species_references['SPECIES_SYMBOL']==sp_id.upper()]

        if to_convention == 'FIA':
            code=str(codes['SPCD'].iloc[0])
        elif to_convention=='VF':
            code=str(codes['VF_SPCD'].iloc[0])

    elif from_convention=='FIA':
        codes = species_references[species_references['SPCD'] == sp_id]

        if to_convention == 'PLANTS':
            code=str(codes['SPECIES_SYMBOL'].iloc[0])
        elif to_convention=='VF':
            code=str(codes['VF_SPCD'].iloc[0])
    return code

def treelist_from_lat_lon(lat,lon,path=None):
    """"Return a treelist based on the USDA TreeMap dataset. Optionally, write to a csv file.
    Arguments:
        lat: latitude
        lon: longitude
        path: path to save csv file to disk
    """

    import requests
    import io

    lat=round(lat,2)
    lon = round(lon, 2)

    url = f'https://treemap-host.uw.r.appspot.com/LAT={lat},LON={lon}'

    response = requests.get(url)

    if path != None:
        with open(path, 'wb') as f:
            f.write(response.content)

            return


    else:
        tlist = pyvf.treelist_from_table(io.BytesIO(response.content))
        return tlist


def stand_table_from_generic(max_dbh, tpa, species_composition, distribution,species_max_dbh=None):
    """Create a stand table from scratch using a generic distribution
    Arguments
    max_dbh: diameter of largest tree in inches
    tpa: total trees per acre of resulting stand
    species composition: dictionary of species code and percentage. E.g. {"DF":.5,"WF":.3,"IC":.2}
    distribution: list of values representing the relative distribution of trees per acre across the dbh range. Values are unitless and relative to other values in the list. E.g. [100,50,25,12,6,3]
    species_max_dbh: list of values representing largest dbh possible for each species in dictionary

    """
    max_dbh = int(max_dbh + 1)
    tpa = float(tpa)
    distribution = np.array(distribution)

    if type(species_composition) is not dict:
        raise TypeError("Species composition must be dictionary")

    # Create starting arrays
    if species_max_dbh is None:
        species_max_dbh=[72]
    else:
        species_max_dbh=list(species_max_dbh)

    if len(species_max_dbh) > len(species_composition):
        raise Exception("Size of species_max_dbh must equal size of species dictionary keys")
    else:
        while len(species_max_dbh)<len(species_composition):
            species_max_dbh.append(72)

    species_max_dbh=np.array(species_max_dbh)
    species_max_dbh = species_max_dbh.clip(0, max_dbh)
    dbh_array = np.linspace(start=0.01, stop=1.01, num=distribution.shape[0])

    # Expand species array and fill
    species_array = np.zeros([len(species_composition.keys()), species_max_dbh.max()])

    j = 0

    for value in species_composition.values():
        for i in range(0, species_array.shape[1]):
            if i < species_max_dbh[j]:
                species_array[j, i] = value
        j = j + 1

    # Normalize species array such that all rows sum to 1.0

    for i in range(0, species_array.shape[1]):
        total = species_array[:, i].sum()
        multiplier = 1 / total
        for j in range(0, len(species_composition.keys())):
            species_array[j, i] = species_array[j, i] * multiplier

    # Scale DBH array using target max dbh. The dbh array could be any length. Need to interpolate values to match the length of the species array.

    dbh_array = dbh_array * max_dbh
    dbh_array = dbh_array.round(0)  # Round to the nearest int so that the values can be matched with the species array

    i = 0
    while i < dbh_array.shape[0] - 1:
        if dbh_array[i] == dbh_array[
            i + 1]:  # Must delete any duplicates that arose after rounding to prevent future issues
            dbh_array = np.delete(dbh_array, i)
            distribution = np.delete(distribution, i)
            i = 0
        i = i + 1

    # dbh_array[0]=1
    print(dbh_array)
    print(distribution)

    output_array = np.arange(0, species_array.shape[
        1])  # Create an ascending array with the same length as the species array
    output_array = np.stack(
        [output_array, np.zeros(species_array.shape[1])])  # Add a 'column' to contain interpolated distribution values

    input_index = -1
    dbh_range = 0

    output_array[1, 0] = distribution[0]

    for output_index in output_array[0]:
        output_index = int(output_index)

        if output_index in dbh_array:

            if input_index + 2 < dbh_array.shape[0]:
                input_index = input_index + 1

            dbh_0 = dbh_array[input_index]
            dbh_1 = dbh_array[input_index + 1]
            dbh_range = dbh_1 - dbh_0

            if dbh_range > 0:

                if output_index == 0:
                    step = (distribution[input_index + 1] - distribution[input_index]) / (dbh_range - 1)

                else:
                    step = (distribution[input_index + 1] - distribution[input_index]) / (dbh_range)
                    output_array[1, output_index] = output_array[1, output_index - 1] + step

            # print(dbh_0)
            # print(dbh_1)
            # print(dbh_range)
            # print(step)
            # print('')
            # print('')

        elif dbh_range > 0:
            output_array[1, output_index] = output_array[1, output_index - 1] + step

    # print(output_array)

    # Scale output to desired TPA
    multiplier = tpa / output_array[1].sum()
    output_array[1] = output_array[1] * multiplier

    # Apply tpa distribution to species distribution

    for j in range(0, species_array.shape[0]):
        species_array[j] = species_array[j] * output_array[1]

    # Format output as dataframe
    output = pd.DataFrame(species_array.transpose(), columns=species_composition.keys())
    # output['Total']=output.sum(axis=1)

    return output

def stem_map_from_table(table,on='headers'):
    """Create a stem map object from a csv or DataFrame containing matching headers or column order. Column names or
    column order must match x,y,vf_species,dbh,height,tree_id,structure,condition,sp_name,sp_id,cubic,scribner,biomass
    Arguments:
        table: path to csv file or DataFrame-like object
        on: associate columns using 'headers' or 'position'

    """
    import pandas as pd
    from pyvf.stem_map import Stem_Map

    try:
        table=pd.read_csv(table)
    except:
        pass

    for column in table.columns:
        table.rename(columns={column: column.lower()}, inplace=True)

    if on =='headers':

        try:
            stems = Stem_Map(table['u'], table['v'], table['vf_species'], table['dbh'])
        except:
            stems = Stem_Map(table['x'], table['y'], table['vf_species'], table['dbh'])

        try:
            stems.height=np.array(table['height'])
        except:
            pass
        try:
            stems.tree_id=np.array(table['tree_id'])
        except:
            pass
        try:
            stems.structure=np.array(table['structure'])
        except:
            pass
        try:
            stems.condition=np.array(table['condition'])
        except:
            pass
        try:
            stems.sp_name=np.array(table['sp_name'])
        except:
            pass
        try:
            stems.sp_id=np.array(table['sp_id'])
        except:
            pass
        try:
            stems.cubic=np.array(table['cubic'])
        except:
            pass
        try:
            stems.scribner=np.array(table['scribner'])
        except:
            pass
        try:
            stems.biomass=np.array(table['biomass'])
        except:
            pass

    else:

        stems=Stem_Map(table.iloc[0],table.iloc[1],table.iloc[2],table.iloc[3])

        try:
            stems.height=np.array(table.iloc[4])
        except:
            pass
        try:
            stems.tree_id=np.array(table.iloc[5])
        except:
            pass
        try:
            stems.structure=np.array(table.iloc[6])
        except:
            pass
        try:
            stems.condition=np.array(table.iloc[7])
        except:
            pass
        try:
            stems.sp_name=np.array(table.iloc[8])
        except:
            pass
        try:
            stems.sp_id=np.array(table.iloc[9])
        except:
            pass
        try:
            stems.cubic=np.array(table.iloc[10])
        except:
            pass
        try:
            stems.scribner=np.array(table.iloc[11])
        except:
            pass
        try:
            stems.biomass=np.array(table.iloc[12])
        except:
            pass
    return stems



def treelist_from_table(table,on='headers'):
    """Create a treelist object from a csv or DataFrame containing matching headers or column order. Column names or
    column order must match tpa,vf_species,dbh,height,tree_id,structure,condition,sp_name,sp_id,cubic,scribner,biomass
    Arguments:
        table: path to csv file or DataFrame-like object
        on: associate columns using 'headers' or 'position'

    """
    import pandas as pd
    from pyvf.treelist import Treelist

    try:
        table=pd.read_csv(table)
    except:
        pass

    for column in table.columns:
        table.rename(columns={column: column.lower()}, inplace=True)

    if on =='headers':


        tlist=Treelist(table['tpa'],table['vf_species'],table['dbh'])

        try:
            tlist.height=np.array(table['height'])
        except:
            pass
        try:
            tlist.tree_id=np.array(table['tree_id'])
        except:
            pass
        try:
            tlist.structure=np.array(table['structure'])
        except:
            pass
        try:
            tlist.condition=np.array(table['condition'])
        except:
            pass
        try:
            tlist.sp_name=np.array(table['sp_name'])
        except:
            pass
        try:
            tlist.sp_id=np.array(table['sp_id'])
        except:
            pass
        try:
            tlist.cubic=np.array(table['cubic'])
        except:
            pass
        try:
            tlist.scribner=np.array(table['scribner'])
        except:
            pass
        try:
            tlist.biomass=np.array(table['biomass'])
        except:
            pass

    else:

        tlist=Treelist(table.iloc[0],table.iloc[1],table.iloc[2])

        try:
            tlist.height=np.array(table.iloc[3])
        except:
            pass
        try:
            tlist.tree_id=np.array(table.iloc[4])
        except:
            pass
        try:
            tlist.structure=np.array(table.iloc[5])
        except:
            pass
        try:
            tlist.condition=np.array(table.iloc[6])
        except:
            pass
        try:
            tlist.sp_name=np.array(table.iloc[7])
        except:
            pass
        try:
            tlist.sp_id=np.array(table.iloc[8])
        except:
            pass
        try:
            tlist.cubic=np.array(table.iloc[9])
        except:
            pass
        try:
            tlist.scribner=np.array(table.iloc[10])
        except:
            pass
        try:
            tlist.biomass=np.array(table.iloc[11])
        except:
            pass
    return tlist

def stem_map_from_vl(path):
    from pyvf.vl_file import VL_File
    vl=VL_File(path)
    return vl.to_stem_map()

def convert_sp_deciduous(vf_species):

    vf_species=vf_species.upper()

    if vf_species in deciduous_sp_dict.keys():
        return deciduous_sp_dict[vf_species]
    else:
        return vf_species



import numpy as np
import pandas as pd
import pyvf

import pkg_resources
species_references=pd.read_csv(pkg_resources.resource_stream(__name__,'data/Species_reference.csv'))
species_references_vf=pd.read_csv(pkg_resources.resource_stream(__name__,'data/Species_reference_vf.csv'))
deciduous_sp_dict = {'BO': 'BOS', 'BM': 'BMS', 'LO': 'LOS', 'QA': 'QAS', 'RA': 'RAS'}

def __get_models_dict():
    import io
    f=pkg_resources.resource_stream(__name__,'data/Tree_models.csv')
    """Feed in models csv with species code, dbh breakpoint, and model name. Return dictionary of species(dbh breakpoint,model name)"""
    models_dict = dict()


    record = f.readline()  # skip header
    while True:
        record = f.readline().decode('UTF-8').strip()
        if len(record) < 1:
            break

        record = str(record).split(',')
        species = record[0]
        breakpoint = record[1]

        try:
            models_dict[species].append((float(record[1]), record[2]))
        except:
            models_dict[species] = list()
            models_dict[species].append((float(record[1]), record[2]))

    return models_dict

models_dict=__get_models_dict()


def convert_sp_deciduous(vf_species):
    vf_species = vf_species.upper()

    if vf_species in deciduous_sp_dict.keys():
        return deciduous_sp_dict[vf_species]
    else:
        return vf_species


def get_model(vf_species, dbh):
    vf_species = convert_sp_deciduous(vf_species)

    for record in models_dict[vf_species]:
        if dbh < record[0]:
            return record[1]+'_A'

    return record[1]+'_A'




def CreateSqGrid(qty_per_ac, variance=0, bbox=[29.2,29.2,158.8,158.8]):
    """Return numpy array of x,y coordinates organized as a square grid with random variance

    Parameters
    qty_per_ac (flt): Point density as number of points per acre
    variance (flt): Set to 0 for a perfect grid, set >1 for very random
    bounds (list): [xmin, ymin, xmax, ymax]

    """
    x_meters = bbox[2] - bbox[0]
    y_meters = bbox[3] - bbox[1]

    # Figure out grid spacing
    acres = ((x_meters * y_meters) * 3.2808399 ** 2) / 43560
    qty = qty_per_ac * acres
    area_per_object = (acres / qty) * 43560 / 3.2808399 ** 2  ## output square meters
    spacing = area_per_object ** .5  ## output meters

    # Figure out number of grid points per axis
    points_x = int(x_meters / spacing)
    points_y = int(y_meters / spacing)

    # Define storage variables
    points = []

    xvals = np.linspace(start=bbox[0], stop=bbox[2], num=(round(x_meters / spacing) + 1))
    yvals = np.linspace(start=bbox[1], stop=bbox[3], num=(round(y_meters / spacing) + 1))
    points = np.array(np.meshgrid(xvals, yvals)).T.reshape(-1, 2)


    # Create an adjustments array with the same shape as the coordinates array. This array represents a random adjustment in meters based on variance as a percent of spacing.

    if variance > .01:
        rand_adj = variance * spacing
        rands = np.random.randint(rand_adj * -1000, rand_adj * 1000, (len(points), 2))
        rands = rands / 1000

        # Apply adjustments
        points = points + rands

    return (points)


def CreateHexGrid(qty_per_ac, variance=0, bbox=[29.2, 29.2, 158.8, 158.8]):
    """Return numpy array of x,y coordinates organized as a hexagonal grid with random variance

    Parameters
    qty_per_ac (flt): Point density as number of points per acre
    variance (flt): Set to 0 for a perfect grid, set >1 for very random
    bounds (list): [xmin, ymin, xmax, ymax]


    """
    x_meters = bbox[2] - bbox[0]
    y_meters = bbox[3] - bbox[1]

    # Figure out hex grid spacing
    acres = ((x_meters * y_meters) * 3.2808399 ** 2) / 43560
    qty = qty_per_ac * acres
    area_per_object = (acres / qty) * 43560 / 3.2808399 ** 2  # Output square meters

    side_radius = (
                              area_per_object / 2.598076) ** .5  # Constant is equivalent to 1.5*(3**.5). Radius, side length, and "size" are all the same thing.
    apothem = side_radius * 0.8660254037844386  # Constant is equivalent to 3**.5/2

    x_spacing = 2 * apothem
    y_spacing = 1.5 * side_radius

    # Figure out number of grid points per axis
    points_x = round(x_meters / x_spacing)
    points_y = round(y_meters / y_spacing)

    yvals = np.linspace(bbox[1], bbox[3], points_y)
    xvals = np.linspace(bbox[0], bbox[2], points_x)

    points = np.array(np.meshgrid(xvals, yvals)).T.reshape(-1, 2)

    filt = np.isin(points[:, 1], yvals[::2])
    points[filt, 0] = points[filt, 0] + apothem

    # Create an adjustments array with the same shape as the coordinates array. This array represents a random adjustment in meters based on variance as a percent of spacing.

    if variance >= .01:
        rand_adj = variance * y_spacing
        rands = np.random.randint(rand_adj * -1000, rand_adj * 1000, (len(points), 2))
        rands = rands / 1000

        # Apply adjustments
        points = points + rands

    return (points)


def CreateClumpsGrid(qty_per_ac, variance=0, bbox=[29.2,29.2,158.8,158.8]):
    """Return numpy array of x,y coordinates organized as a square grid with random variance

    Arguments
    qty_per_ac: Point density as number of points per acre
    variance: Set to 0 for a perfect grid, set >1 for very random
    bounds: [xmin, ymin, xmax, ymax]


    """
    variance = variance + .3

    clump_tpa = .8 * qty_per_ac * (1 + variance) ** 1.3
    grid1_tpa = .8 * qty_per_ac / (1 + variance) ** .8
    clump_percent = .1
    clump_size = 5

    grid1 = CreateHexGrid(grid1_tpa, variance, bbox)

    np.random.shuffle(grid1)

    clump_grid = grid1

    grid_slice = int(clump_percent * grid1.shape[0])

    grid1 = grid1[0:grid_slice, :]

    for point in grid1:
        win = 0
        while win == 0 and clump_size < 100:
            try:

                clump_radius1 = clump_size * np.random.randint(6, 14) / 10
                clump_radius2 = clump_size * np.random.randint(6, 14) / 10

                clump_bbox = [point[0] - clump_radius1, point[1] - clump_radius2, point[0] + clump_radius1,
                              point[1] + clump_radius2]
                clump = CreateHexGrid(clump_tpa, variance, clump_bbox)
                win = 1

            except:
                clump_size = clump_size + 1

        clump = CreateHexGrid(clump_tpa, variance, clump_bbox)
        clump_grid = np.append(clump_grid, clump, 0)
        clump_grid = np.append(clump_grid, clump, 0)

    return (clump_grid)

def treelist_from_stand_table(path,skiprows=0,custom_species=True):
    import pandas as pd
    try:
        df=pd.read_csv(path,index_col=0,header=None,skiprows=skiprows)
    except:
        df=pd.DataFrame(path)
        df=pd.concat([pd.DataFrame([path.columns.values], columns=path.columns.values), path])

    df.drop(index=np.nan,inplace=True,errors='ignore')

    # Import and clean stand table. Remove totals column and totals row. Use DBH column as index. Fill NA with 0.
    df.fillna(0, inplace=True)

    for index in df.index:
        if 'total' in str(index).lower():
            df.drop(index=index, inplace=True, errors='ignore')

    for index in df.index:
        for column in df.columns:
            if 'total' in str(df[column][index]).lower():
                df.drop(columns=column, inplace=True, errors='ignore')


    trees=np.array([[],[],[]])

    if custom_species == True:
        i_vf_species=2
        index = 2

    else:
        i_vf_species=0
        index = 0

    while index+1 < df.index.shape[0]:
        index=index+1
        for column in range(0,df.columns.shape[0]):
            vf_species=df.iloc[i_vf_species,column]
            dbh=df.index[index]
            tpa=df.iloc[index,column]
            tree=np.array([[tpa],[vf_species],[dbh]])
            trees=np.append(trees,tree,1)

    trees = trees.transpose()

    from pyvf.treelist import Treelist
    out=Treelist(trees[:,0],trees[:,1],trees[:,2])
    return out



def convert_species(sp_id, to_convention,from_convention):
    """
    Arguments
    sp_id: a single id code to convert
    to_convention: ID convention to convert to. Valid options: 'FIA', 'PLANTS', or 'VF'
    """

    if from_convention=='VF':

        codes=species_references_vf[species_references_vf['VF_SPCD']==sp_id.upper()]

        if to_convention == 'FIA':
            code=str(codes['SPCD'].iloc[0])
        elif to_convention=='PLANTS':
            code=str(codes['SPECIES_SYMBOL'].iloc[0])

    elif from_convention=='PLANTS':
        codes=species_references[species_references['SPECIES_SYMBOL']==sp_id.upper()]

        if to_convention == 'FIA':
            code=str(codes['SPCD'].iloc[0])
        elif to_convention=='VF':
            code=str(codes['VF_SPCD'].iloc[0])

    elif from_convention=='FIA':
        codes = species_references[species_references['SPCD'] == sp_id]

        if to_convention == 'PLANTS':
            code=str(codes['SPECIES_SYMBOL'].iloc[0])
        elif to_convention=='VF':
            code=str(codes['VF_SPCD'].iloc[0])
    return code

def treelist_from_lat_lon(lat,lon,path=None):
    """"Return a treelist based on the USDA TreeMap dataset. Optionally, write to a csv file.
    Arguments:
        lat: latitude
        lon: longitude
        path: path to save csv file to disk
    """

    import requests
    import io

    lat=round(lat,2)
    lon = round(lon, 2)

    url = f'https://treemap-host.uw.r.appspot.com/LAT={lat},LON={lon}'

    response = requests.get(url)

    if path != None:
        with open(path, 'wb') as f:
            f.write(response.content)

            return


    else:
        tlist = pyvf.treelist_from_table(io.BytesIO(response.content))
        return tlist


def stand_table_from_generic(max_dbh, tpa, species_composition, distribution,species_max_dbh=None):
    """Create a stand table from scratch using a generic distribution
    Arguments
    max_dbh: diameter of largest tree in inches
    tpa: total trees per acre of resulting stand
    species composition: dictionary of species code and percentage. E.g. {"DF":.5,"WF":.3,"IC":.2}
    distribution: list of values representing the relative distribution of trees per acre across the dbh range. Values are unitless and relative to other values in the list. E.g. [100,50,25,12,6,3]
    species_max_dbh: list of values representing largest dbh possible for each species in dictionary

    """
    max_dbh = int(max_dbh + 1)
    tpa = float(tpa)
    distribution = np.array(distribution)

    if type(species_composition) is not dict:
        raise TypeError("Species composition must be dictionary")

    # Create starting arrays
    if species_max_dbh is None:
        species_max_dbh=[72]
    else:
        species_max_dbh=list(species_max_dbh)

    if len(species_max_dbh) > len(species_composition):
        raise Exception("Size of species_max_dbh must equal size of species dictionary keys")
    else:
        while len(species_max_dbh)<len(species_composition):
            species_max_dbh.append(72)

    species_max_dbh=np.array(species_max_dbh)
    species_max_dbh = species_max_dbh.clip(0, max_dbh)
    dbh_array = np.linspace(start=0.01, stop=1.01, num=distribution.shape[0])

    # Expand species array and fill
    species_array = np.zeros([len(species_composition.keys()), species_max_dbh.max()])

    j = 0

    for value in species_composition.values():
        for i in range(0, species_array.shape[1]):
            if i < species_max_dbh[j]:
                species_array[j, i] = value
        j = j + 1

    # Normalize species array such that all rows sum to 1.0

    for i in range(0, species_array.shape[1]):
        total = species_array[:, i].sum()
        multiplier = 1 / total
        for j in range(0, len(species_composition.keys())):
            species_array[j, i] = species_array[j, i] * multiplier

    # Scale DBH array using target max dbh. The dbh array could be any length. Need to interpolate values to match the length of the species array.

    dbh_array = dbh_array * max_dbh
    dbh_array = dbh_array.round(0)  # Round to the nearest int so that the values can be matched with the species array

    i = 0
    while i < dbh_array.shape[0] - 1:
        if dbh_array[i] == dbh_array[
            i + 1]:  # Must delete any duplicates that arose after rounding to prevent future issues
            dbh_array = np.delete(dbh_array, i)
            distribution = np.delete(distribution, i)
            i = 0
        i = i + 1

    # dbh_array[0]=1
    print(dbh_array)
    print(distribution)

    output_array = np.arange(0, species_array.shape[
        1])  # Create an ascending array with the same length as the species array
    output_array = np.stack(
        [output_array, np.zeros(species_array.shape[1])])  # Add a 'column' to contain interpolated distribution values

    input_index = -1
    dbh_range = 0

    output_array[1, 0] = distribution[0]

    for output_index in output_array[0]:
        output_index = int(output_index)

        if output_index in dbh_array:

            if input_index + 2 < dbh_array.shape[0]:
                input_index = input_index + 1

            dbh_0 = dbh_array[input_index]
            dbh_1 = dbh_array[input_index + 1]
            dbh_range = dbh_1 - dbh_0

            if dbh_range > 0:

                if output_index == 0:
                    step = (distribution[input_index + 1] - distribution[input_index]) / (dbh_range - 1)

                else:
                    step = (distribution[input_index + 1] - distribution[input_index]) / (dbh_range)
                    output_array[1, output_index] = output_array[1, output_index - 1] + step

            # print(dbh_0)
            # print(dbh_1)
            # print(dbh_range)
            # print(step)
            # print('')
            # print('')

        elif dbh_range > 0:
            output_array[1, output_index] = output_array[1, output_index - 1] + step

    # print(output_array)

    # Scale output to desired TPA
    multiplier = tpa / output_array[1].sum()
    output_array[1] = output_array[1] * multiplier

    # Apply tpa distribution to species distribution

    for j in range(0, species_array.shape[0]):
        species_array[j] = species_array[j] * output_array[1]

    # Format output as dataframe
    output = pd.DataFrame(species_array.transpose(), columns=species_composition.keys())
    # output['Total']=output.sum(axis=1)

    return output

def stem_map_from_table(table,on='headers'):
    """Create a stem map object from a csv or DataFrame containing matching headers or column order. Column names or
    column order must match x,y,vf_species,dbh,height,tree_id,structure,condition,sp_name,sp_id,cubic,scribner,biomass
    Arguments:
        table: path to csv file or DataFrame-like object
        on: associate columns using 'headers' or 'position'

    """
    import pandas as pd
    from pyvf.stem_map import Stem_Map

    try:
        table=pd.read_csv(table)
    except:
        pass

    for column in table.columns:
        table.rename(columns={column: column.lower()}, inplace=True)

    if on =='headers':

        try:
            stems = Stem_Map(table['u'], table['v'], table['vf_species'], table['dbh'])
        except:
            stems = Stem_Map(table['x'], table['y'], table['vf_species'], table['dbh'])

        try:
            stems.height=np.array(table['height'])
        except:
            pass
        try:
            stems.tree_id=np.array(table['tree_id'])
        except:
            pass
        try:
            stems.structure=np.array(table['structure'])
        except:
            pass
        try:
            stems.condition=np.array(table['condition'])
        except:
            pass
        try:
            stems.sp_name=np.array(table['sp_name'])
        except:
            pass
        try:
            stems.sp_id=np.array(table['sp_id'])
        except:
            pass
        try:
            stems.cubic=np.array(table['cubic'])
        except:
            pass
        try:
            stems.scribner=np.array(table['scribner'])
        except:
            pass
        try:
            stems.biomass=np.array(table['biomass'])
        except:
            pass

    else:

        stems=Stem_Map(table.iloc[0],table.iloc[1],table.iloc[2],table.iloc[3])

        try:
            stems.height=np.array(table.iloc[4])
        except:
            pass
        try:
            stems.tree_id=np.array(table.iloc[5])
        except:
            pass
        try:
            stems.structure=np.array(table.iloc[6])
        except:
            pass
        try:
            stems.condition=np.array(table.iloc[7])
        except:
            pass
        try:
            stems.sp_name=np.array(table.iloc[8])
        except:
            pass
        try:
            stems.sp_id=np.array(table.iloc[9])
        except:
            pass
        try:
            stems.cubic=np.array(table.iloc[10])
        except:
            pass
        try:
            stems.scribner=np.array(table.iloc[11])
        except:
            pass
        try:
            stems.biomass=np.array(table.iloc[12])
        except:
            pass
    return stems



def treelist_from_table(table,on='headers'):
    """Create a treelist object from a csv or DataFrame containing matching headers or column order. Column names or
    column order must match tpa,vf_species,dbh,height,tree_id,structure,condition,sp_name,sp_id,cubic,scribner,biomass
    Arguments:
        table: path to csv file or DataFrame-like object
        on: associate columns using 'headers' or 'position'

    """
    import pandas as pd
    from pyvf.treelist import Treelist

    try:
        table=pd.read_csv(table)
    except:
        pass

    for column in table.columns:
        table.rename(columns={column: column.lower()}, inplace=True)

    if on =='headers':


        tlist=Treelist(table['tpa'],table['vf_species'],table['dbh'])

        try:
            tlist.height=np.array(table['height'])
        except:
            pass
        try:
            tlist.tree_id=np.array(table['tree_id'])
        except:
            pass
        try:
            tlist.structure=np.array(table['structure'])
        except:
            pass
        try:
            tlist.condition=np.array(table['condition'])
        except:
            pass
        try:
            tlist.sp_name=np.array(table['sp_name'])
        except:
            pass
        try:
            tlist.sp_id=np.array(table['sp_id'])
        except:
            pass
        try:
            tlist.cubic=np.array(table['cubic'])
        except:
            pass
        try:
            tlist.scribner=np.array(table['scribner'])
        except:
            pass
        try:
            tlist.biomass=np.array(table['biomass'])
        except:
            pass

    else:

        tlist=Treelist(table.iloc[0],table.iloc[1],table.iloc[2])

        try:
            tlist.height=np.array(table.iloc[3])
        except:
            pass
        try:
            tlist.tree_id=np.array(table.iloc[4])
        except:
            pass
        try:
            tlist.structure=np.array(table.iloc[5])
        except:
            pass
        try:
            tlist.condition=np.array(table.iloc[6])
        except:
            pass
        try:
            tlist.sp_name=np.array(table.iloc[7])
        except:
            pass
        try:
            tlist.sp_id=np.array(table.iloc[8])
        except:
            pass
        try:
            tlist.cubic=np.array(table.iloc[9])
        except:
            pass
        try:
            tlist.scribner=np.array(table.iloc[10])
        except:
            pass
        try:
            tlist.biomass=np.array(table.iloc[11])
        except:
            pass
    return tlist

def stem_map_from_vl(path):
    from pyvf.vl_file import VL_File
    vl=VL_File(path)
    return vl.to_stem_map()

def convert_sp_deciduous(vf_species):

    vf_species=vf_species.upper()

    if vf_species in deciduous_sp_dict.keys():
        return deciduous_sp_dict[vf_species]
    else:
        return vf_species
Feed in models csv with species code, dbh breakpoint, and model name. Return dictionary of species(dbh breakpoint,model name) # skip header Return numpy array of x,y coordinates organized as a square grid with random variance      Parameters     qty_per_ac (flt): Point density as number of points per acre     variance (flt): Set to 0 for a perfect grid, set >1 for very random     bounds (list): [xmin, ymin, xmax, ymax] # Figure out grid spacing ## output square meters ## output meters # Figure out number of grid points per axis # Define storage variables # Create an adjustments array with the same shape as the coordinates array. This array represents a random adjustment in meters based on variance as a percent of spacing. # Apply adjustments Return numpy array of x,y coordinates organized as a hexagonal grid with random variance      Parameters     qty_per_ac (flt): Point density as number of points per acre     variance (flt): Set to 0 for a perfect grid, set >1 for very random     bounds (list): [xmin, ymin, xmax, ymax] # Figure out hex grid spacing # Output square meters # Constant is equivalent to 1.5*(3**.5). Radius, side length, and "size" are all the same thing. # Constant is equivalent to 3**.5/2 # Figure out number of grid points per axis # Create an adjustments array with the same shape as the coordinates array. This array represents a random adjustment in meters based on variance as a percent of spacing. # Apply adjustments Return numpy array of x,y coordinates organized as a square grid with random variance      Arguments     qty_per_ac: Point density as number of points per acre     variance: Set to 0 for a perfect grid, set >1 for very random     bounds: [xmin, ymin, xmax, ymax] # Import and clean stand table. Remove totals column and totals row. Use DBH column as index. Fill NA with 0. Arguments     sp_id: a single id code to convert     to_convention: ID convention to convert to. Valid options: 'FIA', 'PLANTS', or 'VF' "Return a treelist based on the USDA TreeMap dataset. Optionally, write to a csv file.     Arguments:         lat: latitude         lon: longitude         path: path to save csv file to disk Create a stand table from scratch using a generic distribution     Arguments     max_dbh: diameter of largest tree in inches     tpa: total trees per acre of resulting stand     species composition: dictionary of species code and percentage. E.g. {"DF":.5,"WF":.3,"IC":.2}     distribution: list of values representing the relative distribution of trees per acre across the dbh range. Values are unitless and relative to other values in the list. E.g. [100,50,25,12,6,3]     species_max_dbh: list of values representing largest dbh possible for each species in dictionary # Create starting arrays # Expand species array and fill # Normalize species array such that all rows sum to 1.0 # Scale DBH array using target max dbh. The dbh array could be any length. Need to interpolate values to match the length of the species array. # Round to the nearest int so that the values can be matched with the species array # Must delete any duplicates that arose after rounding to prevent future issues # dbh_array[0]=1 # Create an ascending array with the same length as the species array # Add a 'column' to contain interpolated distribution values # print(dbh_0) # print(dbh_1) # print(dbh_range) # print(step) # print('') # print('') # print(output_array) # Scale output to desired TPA # Apply tpa distribution to species distribution # Format output as dataframe # output['Total']=output.sum(axis=1) Create a stem map object from a csv or DataFrame containing matching headers or column order. Column names or     column order must match x,y,vf_species,dbh,height,tree_id,structure,condition,sp_name,sp_id,cubic,scribner,biomass     Arguments:         table: path to csv file or DataFrame-like object         on: associate columns using 'headers' or 'position' Create a treelist object from a csv or DataFrame containing matching headers or column order. Column names or     column order must match tpa,vf_species,dbh,height,tree_id,structure,condition,sp_name,sp_id,cubic,scribner,biomass     Arguments:         table: path to csv file or DataFrame-like object         on: associate columns using 'headers' or 'position'