【Python例項第9講】物種分佈模型
阿新 • • 發佈：2018-11-09
機器學習訓練營——機器學習愛好者的自由交流空間（qq 群號：696721295）
建立物種的地理分佈模型，是保護生物學的一個重要問題。在這個例子裡，我們將根據已有的歷史觀測和14個環境變數，建立兩個南美洲哺乳動物類的地理分佈模型。這兩個物種分別是：
- 褐喉樹懶(
Bradypus variegatus)
- 森林小稻鼠(
Microryzomys minutus)
資料集介紹
本例使用的物種分佈資料集由datasets庫函式fetch_species_distributions載入。它有兩個引數，其中，data_home指定資料集下載後的儲存資料夾，該引數的預設值為None
scikit_learn_data子目錄下。download_if_missing表示如果本地沒有可利用的資料，是否從原始網站下載資料。該引數是邏輯型，預設值為True,如果取False, 則在沒有找到資料時給出一個錯誤提示。
函式fetch_species_distributions返回一個Bunch資料物件，它有屬性：
- coverages: 陣列型，形狀[14, 1592, 1212]
它表示在地圖網格測量的14個特徵的值，其中的缺失值用-9999表示。
- train: 記錄陣列，形狀 (1624,)
它表示資料集的訓練點，每個點有三個域：
- train[‘species’]是物種名字
- train[‘dd long’]是經度
- train[‘dd lat’]是緯度
- test: 記錄陣列，形狀 (620,)
它表示資料的檢驗點，與訓練資料格式相同。
- Nx, Ny: 整型
它們分別表示格點的經度(x), 緯度(y)值。
- x_left_lower_corner, y_left_lower_corner: 浮點型
左下角的座標位置(x, y)
- grid_size: 浮點型
網格上點與點之間的間隔。
例項詳解
首先，載入必需的函式模組和庫。
# Authors: Peter Prettenhofer <[email protected]>
# Jake Vanderplas <[email protected]>
#
# License: BSD 3 clause
from __future__ import print_function
from time import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets.base import Bunch
from sklearn.datasets import fetch_species_distributions
from sklearn.datasets.species_distributions import construct_grids
from sklearn import svm, metrics
# if basemap is available, we'll use it.
# otherwise, we'll improvise later...
try:
from mpl_toolkits.basemap import Basemap
basemap = True
except ImportError:
basemap = False
print(__doc__)
函式 create_species_bunch()
函式create_species_bunch返回一個bunch物件，它描述一個特定物種的資訊。該函式包括一個物種名字的引數，這樣，使用test/train記錄陣列提取指定物種的資料。
def create_species_bunch(species_name, train, test, coverages, xgrid, ygrid):
"""Create a bunch with information about a particular organism
This will use the test/train record arrays to extract the
data specific to the given species name.
"""
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():
# choose points associated with the desired species
pts = pts[pts['species'] == species_name]
bunch['pts_%s' % label] = pts
# determine coverage values for each of the training & testing points
ix = np.searchsorted(xgrid, pts['dd long'])
iy = np.searchsorted(ygrid, pts['dd lat'])
bunch['cov_%s' % label] = coverages[:, -iy, ix].T
return bunch
函式 plot_species_distribution()
函式plot_species_distribution畫出兩個物種，即，Bradypus variegatus, Microryzomys minutus的地理分佈圖。
def plot_species_distribution(species=("bradypus_variegatus_0",
"microryzomys_minutus_0")):
"""
Plot the species distribution.
"""
if len(species) > 2:
print("Note: when more than two species are provided,"
" only the first two will be used")
t0 = time()
# Load the compressed data
data = fetch_species_distributions()
# Set up the data grid
xgrid, ygrid = construct_grids(data)
# The grid in x,y coordinates
X, Y = np.meshgrid(xgrid, ygrid[::-1])
# create a bunch for each species
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)
# background points (grid coordinates) for evaluation
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
# We'll make use of the fact that coverages[6] has measurements at all
# land points. This will help us decide between land and water.
land_reference = data.coverages[6]
# Fit, predict, and plot for each species.
for i, species in enumerate([BV_bunch, MM_bunch]):
print("_" * 80)
print("Modeling distribution of species '%s'" % species.name)
# Standardize features
mean = species.cov_train.mean(axis=0)
std = species.cov_train.std(axis=0)
train_cover_std = (species.cov_train - mean) / std
# Fit OneClassSVM
print(" - fit OneClassSVM ... ", end='')
clf = svm.OneClassSVM(nu=0.1, kernel="rbf", gamma=0.5)
clf.fit(train_cover_std)
print("done.")
# Plot map of South America
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")
# Predict species distribution using the training data
Z = np.ones((data.Ny, data.Nx), dtype=np.float64)
# We'll predict only for the land points.
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
# plot contours of the prediction
plt.contourf(X, Y, Z, levels=levels, cmap=plt.cm.Reds)
plt.colorbar(format='%.2f')
# scatter training/testing points
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')
# Compute AUC with regards to background points
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, 畫物種的地理分佈圖。
plot_species_distribution()
plt.show()
閱讀更多精彩內容，請關注微信公眾號：統計學習與大資料