ふと思い立って「Geo AI」と呼ばれるものを試してみることにした。

Geo AI/Geo-DSVMとは

簡単にいうと簡単にいろいろできる環境をもう用意してあるから、勝手にインスタンスを作ってためしてみてくれよということか。

準備

Geo AI データサイエンス VM の作成

AZUREにサインインし、以下から作成する
https://ms.portal.azure.com/#create/microsoft-ads.geodsvmwindows

東日本にはGPUはないらしいが、とりあえずなのでそのまま作成

まあまあ時間がかかる。

出来上がったらリモートログイン

ArcGIS Proを起動

ArcGIS Proのライセンスを持っているのであれば、それでOK。
ない場合もトライアル版でとりあえず試してみることはできるっぽいのでそれを作成。

Geo 人工知能データサイエンス仮想マシンの使用

チュートリアルの一発目を試してみる
https://docs.microsoft.com/ja-jp/azure/machine-learning/data-science-virtual-machine/use-geo-ai-dsvm
→Python を使用した地理空間分析

上から順に流してみようかと思ったら一発目でエラー

c:\Anaconda\lib\site-packages\sklearn\ensemble\weight_boosting.py:29: DeprecationWarning: numpy.core.umath_tests is an internal NumPy module and should not be imported. It will be removed in a future NumPy release.
  from numpy.core.umath_tests import inner1d

あまり気にしなくてもいいらしい？ので、突き進んでみたら、ちゃんとMAPも表示される

チュートリアルを全部流してみるとこんな結果に。

雰囲気的にはどっかの海岸の海草の存在予測の結果をもとに世界中の海岸に対して予測してみましたって感じなのかな。
とまあ、なんとなく転記しただけなので、まったく理解できていない。
なので、それぞれ理解してみる必要がある。

ただし、その前に インスタンスを停止しておくのを忘れないようにする

Code Reading

Importのあたりはどうってことないので読み飛ばす。と思ったのだが、どうも使っていないImportや同じものが入っていたりと絶妙に気持ち悪い。。

from sklearn.ensemble import RandomForestClassifier
import numpy as NUM
import arcpy as ARCPY
import arcpy.da as DA
import pandas as PD
import seaborn as SEA
import matplotlib.pyplot as PLOT
import pandas as PD
import arcgisscripting as ARC
import arcpy as ARCPY
import SSUtilities as UTILS
from arcgis.gis import *
import os as OS

いらなさそうな変数もある。このIFRAME用の関数もここでしか使っていないので意味があるのかないのかよくわからない。

# 何に使っているのかよくわからない変数
mygis = GIS()

# 指定されたURLのマップをマップビューアでIFRAMEとして出力する関数
from IPython.display import IFrame  
def show_web_map_by_url(URL):
    url = 'http://www.arcgis.com/home/webmap/viewer.html?url=' + URL
    return IFrame(url, width='100%', height=500)

# 解析対象となるマップ？
analysis_map_URL = 'https://services3.arcgis.com/oZfKvdlWHN1MwS48/ArcGIS/rest/services/MachineLearningSeagrass/FeatureServer/1&source=sd'

#Where the seagrasses are in florida
show_web_map_by_url(analysis_map_URL)

直にマップビューアで参照するとこんな感じ

上記はPublicなデータだが、自分の権限内で参照できるデータであれば、なんでもいいんだと思われる。
JoinCountをもっているあたり何かしらのセンサーとかのデータを空間結合して作ったデータなのかな？

解析につかうデータを取得

# 解析につかう属性をピックアップ
# salinity:塩分濃度
# temp:たぶん温度。水温かな
# phosphate:リン酸塩
# nitrate:硝酸塩
# silicate:ケイ酸塩
# dissO2:不明
# NameEMU:不明
# Present:海草の有無？

#Load the features from the shapefile stored in GIS Desktop
inputDir = 'C:/GISDemo/SeaGrass/SeaGrass.gdb/FloridaSeaGrass'
# データをここからとってきているけど、上記のURLの元データと同じってことなんだと思う。。

#Names of Prediction Variables
predictVars = ['salinity', 'temp', 'phosphate','nitrate', 
               'silicate', 'dissO2', 'NameEMU']
#Name of Classification Variable
classVar = ['Present']
#List of all Variables
allVars = predictVars + classVar

#GIS data to Pandas dataframe

# 緯度経度と上で指定した属性値を取得
#Get Data from ArcGIS Desktop
FlLayer = DA.FeatureClassToNumPyArray(inputDir, ["SHAPE@XY"] + allVars)

# 後で使うから測地系も取っておく
#Obtain Spatial Reference
spatRef = ARCPY.Describe(inputDir).spatialReference

# NumPyArrayになっている結果をDataFrameにする
#Define Main Dataframe
data = PD.DataFrame(FlLayer, columns = allVars)
#Display Portion of the Data Frame
data.head()

データのフィールド定義はこんな感じ

データの前処理

# One-Hotベクトル化:ダミー変数を作成
#Process Categorical Data for Analysis flatten the categorical data 
#Create Numeric Fields for One-Hot Encoding of the Categorical Variable
# 対象はNameEMUのみ
catVars = PD.get_dummies(data[predictVars[-1]])

#Remove raw Categories from Dataset
data = data.drop(predictVars[-1], axis = 1)
#Add Processed Categories Back into the Data Frame
data = data.join(catVars)
#Abbreviate Long Categorical Variable Names
newNames = ['c1','c2','c3']
for ind, name in enumerate(newNames):
    data.rename(columns={data.columns[len(predictVars)+ind]:name}, inplace=True)
#Update Predict Variable Names
predictVarsNew = predictVars[:-1] + newNames
#Display Portion of the Data Frame
data.head()

##EVALUATE CORRELATION BETWEEN PREDICTORS- EXCLUDE CATEGORICAL
#Calculate Correlation Coefficient between Prediction Variables
corr = data.drop(data.columns[-3:], axis = 1).astype('float64').corr()

# ここまでで作ったデータをマトリクスとして可視化
# https://qiita.com/hik0107/items/3dc541158fceb3156ee0
#Plot Correlation Matrix Between Prediction Variables
ax = SEA.heatmap(corr, cmap=SEA.diverging_palette(220, 10, as_cmap=True),
            square=True, annot = True, linecolor = 'k', linewidths = 1)
# これは必要なのだろうか？
PLOT.show()

seabornで可視化しているのに、そのあとにPLOT.show()を実行している。
これは意味があるんだろうか？

<参考>
https://note.nkmk.me/python-pandas-get-dummies/
https://qiita.com/haru1977/items/1d67b664f8b5ec48d507
https://blog.shikoan.com/pandas-get-dummies/

学習

# 解析に実際に食わせるデータを作成
##PERFORM RANDOM FOREST CLASSIFICATION
#Fraction of data to be used in Training
fracNum = 0.3
#Seperate the Data into Training and Test Datasets : 適当に3割ぐらいサンプル抽出
train_c1 = data[data[classVar[0]] == 1].sample(frac = fracNum)
train_c0 = data[data[classVar[0]] == 0].sample(frac = fracNum)

#Create Training Dataset
train_set = PD.concat([train_c0, train_c1])

#Create Testing Dataset
test_set = data.drop(train_set.index)
#Encode Seagrass Presence as Classes
indicator, _ = PD.factorize(train_set[classVar[0]])

#Print Test and Train Data Set Sizes
print('Training Data Size = ' + str(train_set.size))
print('Test Data Size = ' + str(test_set.size))


##Create Random Forest Classification Object
# 学習：交差検証ではなくOOB誤り率を使ってる
# n_estimators:決定機の数：500
# oob_score:OOB誤り率を使うかどうか
rfco = RandomForestClassifier(n_estimators = 500, oob_score = True)
#Perform Classification Using Training Set
rfco.fit(train_set[predictVars[:-1]], indicator)
#Plot Variable Importance
#Dummy Variables for Y-Ticks

<参考>
http://hatunina.hatenablog.com/entry/2018/05/27/020720
https://note.nkmk.me/python-pandas-sample/
https://www.haya-programming.com/entry/2018/04/16/234028

推測（とそれの可視化）

# 推測させる
#Predict Seagrass Occurance for the Test Dataset
seagrassPred = rfco.predict(test_set[predictVars[:-1]])
#Calculate Prediction Accuracy
test_seagrass = test_set['Present'].values.flatten()
#Calculate Estimation Error
error = NUM.sum(test_seagrass - seagrassPred)/len(seagrassPred) * 100
# Accuracyを算出。実際に海草がいた数と予測した数の差。本当はF値まで出した方がいいような気もするが。

#Print Accuracy Metrics
print('Accuracy = ' + str(100 - NUM.abs(error)) + ' % ')
print('Locations with Seagrass = ' + str(len(NUM.where(test_seagrass==1)[0])) )
print('Predicted Locations with Seagrass = ' + str(len(NUM.where(seagrassPred==1)[0])))

#Delete the feature  as ARCPY.da.NumPyArrayToFeatureClass does not overwrite 
#ARCPY.DeleteFeatures_management("C:/GISDemo/SeaGrass/SeaGrass.gdb\Florida_Seagrass_Prediction")
# 前に作ったフィーチャクラスが存在する場合はそれを削除する？
# NumPyArrayToFeatureClassは上書きができない？

# 予測した結果の中から海草がありうると予測された地点だけを抽出し保存
##BRING OUTPUT BACK INTO ArcGIS Desktop
#Allow overwriting Feature Classes
ARCPY.env.overwriteOutput = True
#Get Indexes for the Test Dataset
outputDir = r'C:/GISDemo/SeaGrass/SeaGrass.gdb'
nameFC = 'Florida_Seagrass_Prediction_Py'
#Locations with Seagrass
grassExists = FlLayer[["SHAPE@XY"]][test_set.index[NUM.where(seagrassPred==1)]]
# Write Locations with Seagrass to Feature Class
ARCPY.da.NumPyArrayToFeatureClass(grassExists, OS.path.join(outputDir, nameFC), ['SHAPE@XY'], spatRef)

# 実際と予測結果の差の可視化
# 予測値はXYしかないので位置しかない。予測につかった属性値も入れておいた方が親切な気もするが。
predicted_fl_map_url ='https://www.arcgis.com/home/webmap/viewer.html?webmap=2e65c37b14764e7a875efb5b5bad4633'
IFrame(predicted_fl_map_url, width='100%', height=500)

こんな感じで可視化される

対象範囲を広げてみる

# フロリダのデータからフロリダでの予測結果はある程度合ってそうなので、
# USA全土のデータから世界中に対する予測に対象を広げてみる
#Names of Prediction Variables#Import USA Seagrass Data 
inputDir = 'C:/GISDemo/SeaGrass/SeaGrass.gdb/USASeaGrass_training'
predictVars = ['salinity', 'temp', 'phosphate','nitrate', 
               'silicate', 'dissO2']

#Name of Classification Variable
classVar = ['Present']

#List of all Variables
allVars = predictVars + classVar

#Create a Data Object
USAlayer = DA.FeatureClassToNumPyArray(inputDir, ["SHAPE@XY"] + allVars)

#Obtain Spatial Reference
spatRefGlobal = ARCPY.Describe(inputDir).spatialReference

#Define Main Dataframe
USA_Train = PD.DataFrame(USAlayer, columns = allVars)

#Display Portion of the Data Frame
USA_Train.head()

##Train Random Forest Using USA Data
#Encode Seagrass Presence as Classes
indicatorUSA, _ = PD.factorize(USA_Train['Present'])

#Create Random Forest Classification Object
rfco = RandomForestClassifier(n_estimators = 500)

#Perform Classification Using Training Set
rfco.fit(USA_Train[predictVars[:-1]], indicatorUSA)
predictVars[:-1]

全米データから学習。

世界に対して推測

## Import Global Data for Prediction ##
##globalUrl = r'https://services3.arcgis.com/oZfKvdlWHN1MwS48/arcgis/rest/services/MachineLearningSeagrass/FeatureServer/2?token=W4CSlYAQxBZKNONRrS_MVm0qNa7j21xD5BNMByQOuXTQ6xhXxFydToukAXYrndK0RRrMNjvbFHcbCOioIVNAgIe6f0iV7rclZ4kNGqoI3i1Ldjhv_UgwzeBGc7KJacuAshcwYXkb6PIgrGipbwD3YOV3AGhhRmh9e5RzMSWB7HJS9fSYzU2l8ZL9vkLmWccKvo_DGJaqm9GIkz9HoPM2hv2yGDdf9SEtcoQvBsNkCWk.'
globalUrl = 'C:/GISDemo/SeaGrass/SeaGrass.gdb/Global_Seagrass_Prediction'

#Create a SS Data Object
globalData = DA.FeatureClassToNumPyArray(globalUrl, ["SHAPE@XY"] + predictVars[:-1])

#Obtain Spatial Reference
spatRefGlobal = ARCPY.Describe(globalUrl).spatialReference

#Define Dataframe
globalTrain = PD.DataFrame(globalData, columns = predictVars[:-1])

#Predict Global Seagrass Occurance
seagrassPredGlobal = rfco.predict(globalTrain)

#Bring Output Back to ArcGIS
nameFC = 'GlobSeagrass'
outputDir = r'C:/GISDemo/SeaGrass/SeaGrass.gdb'
nameFC = 'Global_Seagrass_Prediction_Done'
grassExists = globalData[["SHAPE@XY"]][globalTrain.index[NUM.where(seagrassPredGlobal==1)]]

# Write Locations with Seagrass to Feature Class
ARCPY.da.NumPyArrayToFeatureClass(grassExists, OS.path.join(outputDir, nameFC), ['SHAPE@XY'], spatRefGlobal)

predicted_gl_map_url='https://www.arcgis.com/home/webmap/viewer.html?webmap=503efd8acb7b48e3b92b4f4ac4d6e82e'
IFrame(predicted_gl_map_url, width='100%', height=500)

こんな感じで可視化される

本当であれば、このように位置を予測されたデータをさらに空間解析を行うことでGISのデータとして活用するんだろうなとは思うが、このチュートリアルではここまで。

感想

・・・丹念に読み解いてみた結果、これは全然「Geo AI」じゃなくね？？というのが率直な感想。
属性値をランダムフォレストにかけて、その属性が持っている位置情報（XY）で可視化しただけじゃないのかな。
※どちらかというと、この属性情報を作る際に空間結合などのGIS的な処理が必要になるような気がするけど、「GeoAI」として期待しているのは「位置情報」そのものも学習の対象にするってことなきがするんだよな。。。

チュートリアルをもう少し進めるとそういうことでも出てくるんだろうか？