tdicap/FitMLModel.py at master · berdaniera/tdicap · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
##################################3
from osgeo import gdal, gdalnumeric, ogr
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pa
import pickle
import re
import seaborn as sea
import simplejson
from sklearn.cross_validation import cross_val_score, KFold
from sklearn import neighbors
from sklearn import metrics
from sklearn.externals import joblib

# Bring in classified points
th = open("/home/aaron/Desktop/treehousedata.txt")
xo = simplejson.load(th)

# Bring in raster
fi = '/home/aaron/Desktop/2010448_1N1E08U4BAND/1n1e08u_4band.tif'
ras = gdal.Open(fi)
rasdat = ras.ReadAsArray()


# extract values
dat = pa.DataFrame(xo,columns=['type','x','y'])
dat.x = dat.x.astype(int)
dat.y = dat.y.astype(int)
sr = np.log(rasdat[3]*1./rasdat[0])
trees = sr[dat.x[dat.type=='T'],dat.y[dat.type=='T']]
houss = sr[dat.x[dat.type=='H'],dat.y[dat.type=='H']]

# plot histogram
bins = np.arange(-2,4,0.1)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(houss,bins=bins,alpha=0.6,label="Houses")
ax.hist(trees,bins=bins,alpha=0.6,label="Trees")
ax.set_xlim([-0.5,2.5])
ax.legend(loc='upper right',fontsize=20)
ax.set_xlabel('Spectral index (NIR / Red)',fontsize=20)
ax.tick_params(axis='both', which='major', labelsize=15)
#fig.show()
fig.savefig('/home/aaron/Desktop/outdata.pdf',format="pdf")

#sea.kdeplot(houss,shade=True,cut=0,bw=0.05)
#sea.kdeplot(trees,shade=True,cut=0,bw=0.05)
#sea.plt.show()

# KNN analysis
allval = sr[dat.x,dat.y]

nn = neighbors.KNeighborsClassifier(n_neighbors=10)
nn.fit(allval.reshape(-1,1),dat.type)
y_obs = dat.type
y_pred = nn.predict(allval.reshape(-1,1))

# save the model for loading later
joblib.dump(nn,'/home/aaron/Desktop/model/treemodel.pkl')
#load the model
# nn = joblib.load('/home/aaron/Desktop/model/treemodel.pkl')

scores = cross_val_score(nn,allval.reshape(-1,1),dat.type,cv=10,scoring='accuracy')
scores.mean()

metrics.confusion_matrix(y_obs,y_pred)

######################
# PREDICT WHOLE IMAGE
# Data for all points0
srr = sr.reshape(-1,1)
srr[np.isnan(srr)] = 0
srr[np.isinf(srr)] = 0
np.min(srr)
# k fold prediction

folds = np.array_split(srr,10)

classout = np.array([])
for f in range(len(folds)):
    classout = np.concatenate( (classout,nn.predict(folds[f])) )
    print (f+1.)/10

classout = classout.reshape(sr.shape)
classout.shape
#classout[:10,:10]

mask = np.zeros(classout.shape).astype(int)
mask[classout=='H'] = 255

# matrix to mask
realcol = np.array([rasdat[0,:,:],rasdat[1,:,:],rasdat[2,:,:]])
maskedcol = np.array([rasdat[0,:,:],rasdat[1,:,:],rasdat[2,:,:],mask])/255.

fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(realcol.T[1000:2500,:3000,:])
#fig.show()
fig.savefig('/home/aaron/Desktop/MapOrig.pdf',format="pdf")

fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(maskedcol.T[1000:2500,:3000,:])
#fig.show()
fig.savefig('/home/aaron/Desktop/MapMask.pdf',format="pdf")