Estimação de Densidade Kernel de Distribuições de Espécies

Introdução

Este laboratório demonstra um exemplo de uma consulta baseada em vizinhos (em particular uma estimativa de densidade de kernel) em dados geoespaciais, utilizando uma Árvore de Bolas construída sobre a métrica de distância Haversine -- ou seja, distâncias entre pontos em latitude/longitude. O conjunto de dados é fornecido por Phillips et. al. (2006). O exemplo utiliza a biblioteca basemap para representar as linhas costeiras e as fronteiras nacionais da América do Sul.

Dicas da Máquina Virtual

Após o arranque da máquina virtual, clique no canto superior esquerdo para mudar para a aba Notebook para aceder ao Jupyter Notebook para a prática.

Por vezes, pode ser necessário esperar alguns segundos para o Jupyter Notebook terminar de carregar. A validação das operações não pode ser automatizada devido a limitações no Jupyter Notebook.

Se tiver problemas durante o aprendizado, não hesite em contactar o Labby. Forneça feedback após a sessão e resolveremos prontamente o problema para si.

Importação de Bibliotecas Necessárias

O primeiro passo é importar as bibliotecas necessárias para este laboratório. Neste laboratório, usaremos as bibliotecas numpy, matplotlib, fetch_species_distributions e KernelDensity.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_species_distributions
from sklearn.neighbors import KernelDensity

Carregar Dados

O próximo passo é carregar os dados fornecidos por Phillips et. al. (2006). O conjunto de dados contém duas espécies que usaremos para demonstrar a estimativa de densidade de kernel.

data = fetch_species_distributions()
species_names = ["Bradypus Variegatus", "Microryzomys Minutus"]

Preparar Dados

Agora, prepararemos os dados para a estimativa de densidade de kernel. Extrairemos as informações de latitude e longitude do conjunto de dados e as convertemos para radianos.

Xtrain = np.vstack([data["train"]["dd lat"], data["train"]["dd long"]]).T
ytrain = np.array(
    [d.decode("ascii").startswith("micro") for d in data["train"]["species"]],
    dtype="int",
)
Xtrain *= np.pi / 180.0  ## Converter lat/long para radianos

Construir Malhas

Agora, construiremos a malha do mapa a partir do objeto de lote. Usaremos a função construct_grids para isso.

def construct_grids(batch):
    """Construir a malha do mapa a partir do objeto de lote

    Parâmetros
    ----------
    batch : Objeto Batch
        O objeto retornado por :func:`fetch_species_distributions`

    Retorna
    -------
    (xgrid, ygrid) : arrays 1-D
        A malha correspondente aos valores em batch.coverages
    """
    ## Coordenadas x,y para células de canto
    xmin = batch.x_left_lower_corner + batch.grid_size
    xmax = xmin + (batch.Nx * batch.grid_size)
    ymin = batch.y_left_lower_corner + batch.grid_size
    ymax = ymin + (batch.Ny * batch.grid_size)

    ## Coordenadas x das células da malha
    xgrid = np.arange(xmin, xmax, batch.grid_size)
    ## Coordenadas y das células da malha
    ygrid = np.arange(ymin, ymax, batch.grid_size)

    return (xgrid, ygrid)

## Chamar a função e armazenar os resultados em xgrid e ygrid
xgrid, ygrid = construct_grids(data)

Preparar a Grade de Dados

Vamos configurar a grade de dados para o gráfico de contorno. Usaremos a função construct_grids para isso.

X, Y = np.meshgrid(xgrid[::5], ygrid[::5][::-1])
land_reference = data.coverages[6][::5, ::5]
land_mask = (land_reference > -9999).ravel()

xy = np.vstack([Y.ravel(), X.ravel()]).T
xy = xy[land_mask]
xy *= np.pi / 180.0

Plotar Mapa da América do Sul

Agora, plotaremos o mapa da América do Sul com as distribuições de cada espécie.

fig = plt.figure()
fig.subplots_adjust(left=0.05, right=0.95, wspace=0.05)

for i in range(2):
    plt.subplot(1, 2, i + 1)

    print(" - calculando KDE em coordenadas esféricas")
    kde = KernelDensity(
        bandwidth=0.04, metric="haversine", kernel="gaussian", algorithm="ball_tree"
    )
    kde.fit(Xtrain[ytrain == i])

    Z = np.full(land_mask.shape[0], -9999, dtype="int")
    Z[land_mask] = np.exp(kde.score_samples(xy))
    Z = Z.reshape(X.shape)

    levels = np.linspace(0, Z.max(), 25)
    plt.contourf(X, Y, Z, levels=levels, cmap=plt.cm.Reds)

    if basemap:
        print(" - plotar linhas costeiras usando basemap")
        m = Basemap(
            projection="cyl",
            llcrnrlat=Y.min(),
            urcrnrlat=Y.max(),
            llcrnrlon=X.min(),
            urcrnrlon=X.max(),
            resolution="c",
        )
        m.drawcoastlines()
        m.drawcountries()
    else:
        print(" - plotar linhas costeiras a partir da cobertura")
        plt.contour(
            X, Y, land_reference, levels=[-9998], colors="k", linestyles="solid"
        )
        plt.xticks([])
        plt.yticks([])

    plt.title(species_names[i])

plt.show()

Resumo

Neste laboratório, aprendemos como realizar a estimativa de densidade de kernel em dados geoespaciais. Usamos o conjunto de dados de Phillips et al. (2006) para demonstrar essa técnica. Também aprendemos a plotar o mapa da América do Sul com as distribuições de cada espécie utilizando a biblioteca basemap.