Diagramación de una Chacra para Forestería Análoga

© Equipo I+D SomosAZUCAR - Bajo licencia AGPLv3

Diagramador de Planos para Forestería Análoga
Copyright (C) 2015 Sebastian Silva

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Affero General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Affero General Public License for more details.

You should have received a copy of the GNU Affero General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

In [1]:
# Bibliotecas utilizadas para confeccionar el mapa
%matplotlib inline
import matplotlib.pyplot as plt
from descartes import PolygonPatch
import matplotlib.cm as cmx
import matplotlib.colors as colors
from shapely import geometry
from shapely import ops
import numpy as np
import pandas as pd
import random
import mpld3
from tqdm import tqdm

1.- DATOS INICIALES

© Arbio - licencia por definir

La base de datos consiste en una lista de especies, cantidades, diámetro, y características.


In [2]:
especies = pd.read_csv("db_final.csv") ## Es una base
especies = especies[(especies.Usos=="Medicinal")]
especies = especies.sort(['Altura.tot'], ascending=False)
#especies = especies.tail(10)

In [3]:
consolidado = pd.DataFrame({'especies':especies.groupby('Estrato').size(), 
                            'individuos':especies.groupby('Estrato')['ind'].sum()})
pd.concat([consolidado,pd.DataFrame(consolidado.sum(axis=0),columns=['Total']).T])


Out[3]:
especies individuos
Estrato
1 4 612
2 6 132
3 8 56
4 7 28
5 2 8
Total 27 836

2.- LINDEROS DE LA CHACRA

Los linderos están definidos como un polígono en metros, desde una de las esquinas.


In [4]:
chacra = geometry.Polygon([ [0,0],
                            [0, 93],
                            [55, 93],
                            [55, 86],
                            [35, 65],
                            [37.5, 47.5],
                            [55, 44],
                            [65, 52.5],
                            [65, 0] ])
chacra


Out[4]:

3.- Definir ubicaciones


In [5]:
def punto_aleatorio(poligono):
    (minx, miny, maxx, maxy) = poligono.bounds
    while True:
        p = geometry.Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
        if poligono.contains(p):
            return(p)

In [6]:
def get_colormap(N):
    '''Returns a function that maps each index in 0, 1, ... N-1 to a distinct 
    RGB color.'''
    color_norm  = colors.Normalize(vmin=0, vmax=N-1)
    scalar_map = cmx.ScalarMappable(norm=color_norm, cmap='hsv') 
    def map_index_to_rgb_color(index):
        return scalar_map.to_rgba(index)
    return map_index_to_rgb_color
# http://stackoverflow.com/questions/14720331/how-to-generate-random-colors-in-matplotlib

In [7]:
def regla_1(poblacion, planta):
    # Regla 1: La distancia al lindero debe ser mayor que 3/4 de la altura
    if planta.pos.distance(poblacion.poligono.exterior) < planta['Altura.tot'] * 0.7:        raise ValueError('regla1')
                
def regla_2(poblacion, planta):
    for anterior in poblacion.individuos.itertuples():
        if planta.pos.distance(anterior[2]) < planta.Distancia / 2:
            raise ValueError('regla2')
            
def regla_3(poblacion, planta):
    mi_estrato = int(especies.Estrato[planta.id])
    mas_cercanos = []
    for anterior in poblacion.individuos.itertuples():
        if not mas_cercanos:
            mas_cercanos = [anterior]
        for cercano in mas_cercanos:
            if planta.pos.distance(anterior[2]) < planta.pos.distance(cercano[2]):
                mas_cercanos.append(anterior)
                if len(mas_cercanos) == 2:
                    mas_cercanos.pop(0)
                
    for cercano in mas_cercanos:
        if abs(especies.Estrato[int(cercano[1])] - mi_estrato) > 1:
            raise ValueError('regla3')

In [8]:
class Poblacion:
    def __init__(self, poligono, especies):
        self.poligono = poligono
        self.especies = especies
        individuos_columns = ['id', 'pos', 'x', 'y', 'color', 'Diametro.dosel', 
                              'Nombre', 'Estrato', 'Diametro.punto']
        self.individuos = pd.DataFrame( columns=individuos_columns, )
        self.n_individuos = 0
        
        self.excluidos = pd.DataFrame( columns=individuos_columns, )
        self.inubicables = 0
        self.errores = []
        
        colores = get_colormap(len(especies))
        
        i=0
        for especie in tqdm(list(especies.itertuples())):
            i=i+1
            n_especie = especie[0]
            color_especie = colores(i)
            planta_tipo = especies.loc[n_especie]
            for n_individuo in range(especies.ind[n_especie]):
                self.n_individuos = self.n_individuos + 1
                planta = planta_tipo.copy()
                planta['id'] = n_especie
                planta = self.ubicar(planta)
                if planta['pos']:
                    planta['color'] = color_especie
                    planta['Diametro.dosel'] = float(planta['Diametro.dosel'])
                    planta['Diametro.punto'] = float(planta['Diametro.dosel']*10)
                    self.individuos.loc[self.n_individuos] = planta
                else:
                    self.excluidos.loc[self.inubicables] = planta
                self.individuos = self.individuos.dropna()
        
        #self.individuos = self.individuos.sort(['Estrato'], ascending=False)

        if self.inubicables:
            print ("No se pudieron ubicar " + str(self.inubicables) + " individuos.")
            #self.inubicables = pd.DataFrame(self.inubicables)
                        
    def ubicar(self, planta):
        intentos = 500
        while intentos:
            intentos = intentos - 1
            distancia_min = 1 # 1m
                
            pos = punto_aleatorio(self.poligono)       
            planta['pos'] = pos
            planta['x'], planta['y'] = pos.xy
            
            try:
                regla_1(self, planta)
                regla_2(self, planta)
                regla_3(self, planta)
            except ValueError as e:
                self.errores.append(str(e))
                next
            else:
                return planta
            
        self.inubicables = self.inubicables + 1
        planta['pos'] = None          
        return planta

In [9]:
# Demora bastante
pob = Poblacion(chacra, especies)


                                               
No se pudieron ubicar 5 individuos.


In [111]:
print ("Se ubicaron %s de %s individuos." % 
        ( len(pob.individuos), especies.ind.sum(axis=0)))


Se ubicaron 830 de 836 individuos.

In [112]:
errores = pd.DataFrame({'error':pob.errores})
pd.DataFrame({'instancias':errores.groupby('error').size()})


Out[112]:
instancias
error
regla1 5581
regla2 925
regla3 372

In [113]:
pd.DataFrame({'excluidos':pob.excluidos.groupby('Nombre').size()})


Out[113]:
excluidos
Nombre
Aguaje 2
Coco 1
Copaiba 1
Mango 1
Pan de arbol 1

4.- Graficar el mapa

Empezamos por los linderos de la chacra.


In [10]:
fig = plt.figure(figsize=(6, 6))
chacra_x, chacra_y = chacra.exterior.xy
ax = fig.add_subplot(111)
ax.set_ylim(-10,100)
ax.set_xlim(-20,90)
ax.plot( *chacra.exterior.xy )

circulos = []
colores = []
puntos_x = []
puntos_y = []
for planta in pob.individuos.iterrows():
    circulos.append(planta[1]['pos'].buffer(planta[1]['Diametro.dosel']))
    colores.append(planta[1]['color'])
    puntos_x.append(planta[1]['x'])
    puntos_y.append(planta[1]['y'])

i=0
for poligono in circulos:
    patch = PolygonPatch(poligono, fc=colores[i-1], alpha=0.2, zorder=1)
    ax.add_patch(patch)
    i = i+1
    
ax.grid(color='gray', alpha=0.9)



In [11]:
ax.scatter(list(pob.individuos.x), 
        list(pob.individuos.y), 
        color=list(pob.individuos.color), 
        s=20, 
        alpha=0.6, marker="+")

# Gráfico interactivo
mpld3.display(fig)


Out[11]:

In [147]:
fig = plt.figure(figsize=(36, 36))
chacra_x, chacra_y = chacra.exterior.xy
ax = fig.add_subplot(111)
ax.set_ylim(-10,100)
ax.set_xlim(-20,90)
ax.plot( *chacra.exterior.xy )

i=0
for poligono in circulos:
    patch = PolygonPatch(poligono, fc=colores[i-1], alpha=0.4, zorder=1)
    ax.add_patch(patch)
    i = i+1

ax.scatter(list(pob.individuos.x), 
        list(pob.individuos.y), 
        color=[0,0,0], 
        s=50, 
        alpha=1, marker="+")
ax.grid(color='gray', alpha=1)
    
fig.savefig('mapa.png')
#
# *-* The End *-*



In [ ]:


In [ ]: