物种分布模型

Posted 2023-04-21

参考技术A 3、物种分布模型

from time import time

import numpy as np

import matplotlib.pyplot as plt

from sklearn.utils import Bunch

from sklearn.datasets import fetch_species_distributions

from sklearn import svm, metrics

try:

    from mpl_toolkits.basemap import Basemap

    basemap = True

except ImportError:

    basemap = False

def construct_grids(batch):

    # x,y角单元格坐标

    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)

    # 网格单元格的x坐标

    xgrid = np.arange(xmin, xmax, batch.grid_size)

    # 网格单元格的y坐标

    ygrid = np.arange(ymin, ymax, batch.grid_size)

    return (xgrid, ygrid)

def create_species_bunch(species_name, train, test, coverages, xgrid, ygrid):

    bunch = Bunch(name=" ".join(species_name.split("_")[:2]))

    species_name = species_name.encode("ascii")

    points = dict(test=test, train=train)

    for label, pts in points.items():

        # 选择与所需物种相关的点

        pts = pts[pts["species"] == species_name]

        bunch["pts_%s" % label] = pts

        #确定每个培训和测试点的复盖值

        ix = np.searchsorted(xgrid, pts["dd long"])

        iy = np.searchsorted(ygrid, pts["dd lat"])

        bunch["cov_%s" % label] = coverages[:, -iy, ix].T

    return bunch

def plot_species_distribution(

    species=("bradypus_variegatus_0", "microryzomys_minutus_0")

):

    #绘制物种分布图

    if len(species) > 2:

        print(

            "Note: when more than two species are provided,"

            " only the first two will be used"

        )

    t0 = time()

    #加载压缩数据

    data = fetch_species_distributions()

    # 设置数据网格

    xgrid, ygrid = construct_grids(data)

    # x，y坐标下的网格

    X, Y = np.meshgrid(xgrid, ygrid[::-1])

    # 创建一束物种

    BV_bunch = create_species_bunch(

        species[0], data.train, data.test, data.coverages, xgrid, ygrid

    )

    MM_bunch = create_species_bunch(

        species[1], data.train, data.test, data.coverages, xgrid, ygrid

    )

    #用于评估的背景点（网格坐标）

    np.random.seed(13)

    background_points = np.c_[

        np.random.randint(low=0, high=data.Ny, size=10000),

        np.random.randint(low=0, high=data.Nx, size=10000),

    ].T

    land_reference = data.coverages[6]

    # 对每个物种进行拟合、预测和绘图。

    for i, species in enumerate([BV_bunch, MM_bunch]):

        print("_" * 80)

        print("Modeling distribution of species '%s'" % species.name)

        # 标准化特征

        mean = species.cov_train.mean(axis=0)

        std = species.cov_train.std(axis=0)

        train_cover_std = (species.cov_train - mean) / std

        # 拟合 OneClassSVM

        print(" - fit OneClassSVM ... ", end="")

        clf = svm.OneClassSVM(nu=0.1, kernel="rbf", gamma=0.5)

        clf.fit(train_cover_std)

        print("done.")

        # 南美洲地形图

        plt.subplot(1, 2, i + 1)

        if basemap:

            print(" - plot coastlines using 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(" - plot coastlines from coverage")

            plt.contour(

                X, Y, land_reference, levels=[-9998], colors="k", linestyles="solid"

            )

            plt.xticks([])

            plt.yticks([])

        print(" - predict species distribution")

        # 利用训练数据预测物种分布

        Z = np.ones((data.Ny, data.Nx), dtype=np.float64)

        # 只预测陆地点

        idx = np.where(land_reference > -9999)

        coverages_land = data.coverages[:, idx[0], idx[1]].T

        pred = clf.decision_function((coverages_land - mean) / std)

        Z *= pred.min()

        Z[idx[0], idx[1]] = pred

        levels = np.linspace(Z.min(), Z.max(), 25)

        Z[land_reference == -9999] = -9999

        # 绘制预测的等高线

        plt.contourf(X, Y, Z, levels=levels, cmap=plt.cm.Reds)

        plt.colorbar(format="%.2f")

        # 分散（训练/测试）点

        plt.scatter(

            species.pts_train["dd long"],

            species.pts_train["dd lat"],

            s=2**2,

            c="black",

            marker="^",

            label="train",

        )

        plt.scatter(

            species.pts_test["dd long"],

            species.pts_test["dd lat"],

            s=2**2,

            c="black",

            marker="x",

            label="test",

        )

        plt.legend()

        plt.title(species.name)

        plt.axis("equal")

        # 计算关于背景点的AUC

        pred_background = Z[background_points[0], background_points[1]]

        pred_test = clf.decision_function((species.cov_test - mean) / std)

        scores = np.r_[pred_test, pred_background]

        y = np.r_[np.ones(pred_test.shape), np.zeros(pred_background.shape)]

        fpr, tpr, thresholds = metrics.roc_curve(y, scores)

        roc_auc = metrics.auc(fpr, tpr)

        plt.text(-35, -70, "AUC: %.3f" % roc_auc, ha="right")

        print("\n Area under the ROC curve : %f" % roc_auc)

    print("\ntime elapsed: %.2fs" % (time() - t0))

plot_species_distribution()

使用随机森林模型提高预测新数据的速度

【中文标题】使用随机森林模型提高预测新数据的速度

【英文标题】：Improving the speed of predicting new data using a Random Forest Model

【发布时间】：2014-07-28 06:49:41

【问题描述】：

我正在使用随机森林生成物种分布模型。这些模型试图根据各种环境属性预测一个物种发生的概率。对于大多数物种，我们的初始潜在预测变量集介于 10 到 25 个之间，每个预测变量都由一个包含 460,000,000 个像元的 GIS 栅格文件表示。由于训练数据的性质，我不会在这里详细介绍，我实际上是在基于数据子集构建多个 RF 模型（大约 10 到 100+），然后结合起来为每个物种创建我的整体模型.实际上构建模型花费的时间相对较少（通常几分钟或更短），但使用 predict 函数基于此模型生成预测概率的栅格层可能需要 20 多个小时。我怀疑这个漫长的过程大部分是由于读取/写入大型光栅文件，瓶颈可能是硬盘读取/写入速度。

为了提供更多细节...一旦我有了训练有素的模型，我将通过 raster 包创建代表预测层的层的栅格堆栈，然后使用 predict() 预测该堆栈光栅包中的函数。我有一个相当强大的桌面（Core i7，3.5GHz，w/32 GB RAM），输入和输出光栅文件在本地硬盘上，而不是通过网络移动。我看到 mbq 的回答 here 提供了关于使用 randomForest 加速模型生成的有用建议，并且正在寻找类似的建议来加速预测操作。我能想到许多可能有帮助的事情（例如，种植少量的树，使用其中一个库进行并行处理），我计划在时间允许的情况下测试这些，但它是如果问题主要是读写瓶颈，我不清楚这些是否会产生重大影响。我将不胜感激任何建议。

【问题讨论】：

我只是提供一个轶事，以这种方式在栅格上进行数学运算非常慢，所以这可能确实是您的问题。

怎么样：blog.mckuhn.de/2013/09/…

我确实看到了有关 parallelRandomForest 的帖子。它似乎只支持回归，而不是分类。此外，目前尚不清楚这是否真的会加快预测速度，因为重点似乎是在训练上。读/写速度可能仍然存在瓶颈。

【参考方案1】：

您可以查看mctune function here。这使用 e1071 包来找到最佳参数。但是，您可以对其进行调整以满足您的需求。

  source( './mctune.R')
  rf_ranges = list(ntree=c(seq(1,1000,100),seq(1000,8000,500)),
  mtry=seq(5,15,2))

  set.seed(10)
  tuned.rf = mctune(method = randomForest, train.x = formula1, 
  data = dataframe, tunecontrol = tune.control(sampling = "cross",cross = 5),   
  ranges=rf_ranges,mc.control=list(mc.cores=16, mc.preschedule=T),confusionmatrizes=T )
  save(tuned.rf, file = paste('./tuned_rf.RData',sep='') )

  tuned.rf$best.model
  plot(tuned.rf)

另一个选项可能是使用 doparallel 包中的 foreach (see here)。您可以将每个数据子集（用于新的 RF 模型）分配给每个内核：

RF_outputs = foreach(i=1:length(yourdatasubsets), .inorder=F, .package=c(randomForest)) %dopar% 
    set.seed(10)
    rf <- randomForest(formula, data=na.omit(yourdatasubsets[i]), ntree=2000, proximity=T)
    return(rf)
 

每个经过训练的 RF 模型都将作为 RF_outputs 列表的一部分返回给您。所以RF_outputs[[1]] 将是您的第一个训练有素的 RF 模型。