近來來研究了一下農業耕作問題，尋找求解農機作業面積的方法。多年后編程，一半抄襲，一半原創，可供參考

#coding:GB2312 from sklearn.datasets.samples_generator import make_blobs from sklearn.cluster import DBSCAN from sklearn.cluster import KMeans from sklearn import metrics from sklearn.preprocessing import StandardScaler from sklearn.metrics import silhouette_score, silhouette_samples import numpy as np import matplotlib.pyplot as plt from itertools import cycle ##python自帶的迭代器模塊 import numpy as np import pandas as pd import math##本程序解決經緯度聚集點的區域劃分問題，用來解決農業實際耕作土地的區域劃分，主要應用： ##1、通過聚類分析算法評價耕作區域數量，并根據數量劃分耕作區域 ##2、通過凸包算法計算耕作區域的面積和周長 class area_algo:def __init__(self,queue=None):self.queue =queuedef LatLon2XY(self, latitude, longitude):a = 6378137.0# b = 6356752.3142# c = 6399593.6258# alpha = 1 / 298.257223563e2 = 0.0066943799013# epep = 0.00673949674227#將經緯度轉換為弧度latitude2Rad = (math.pi / 180.0) * latitudebeltNo = int((longitude + 1.5) / 3.0) #計算3度帶投影度帶號L = beltNo * 3 #計算中央經線l0 = longitude - L #經差tsin = math.sin(latitude2Rad)tcos = math.cos(latitude2Rad)t = math.tan(latitude2Rad)m = (math.pi / 180.0) * l0 * tcoset2 = e2 * pow(tcos, 2)et3 = e2 * pow(tsin, 2)X = 111132.9558 * latitude - 16038.6496 * math.sin(2 * latitude2Rad) + 16.8607 * math.sin(4 * latitude2Rad) - 0.0220 * math.sin(6 * latitude2Rad)N = a / math.sqrt(1 - et3)x = X + N * t * (0.5 * pow(m, 2) + (5.0 - pow(t, 2) + 9.0 * et2 + 4 * pow(et2, 2)) * pow(m, 4) / 24.0 + (61.0 - 58.0 * pow(t, 2) + pow(t, 4)) * pow(m, 6) / 720.0)y = 500000 + N * (m + (1.0 - pow(t, 2) + et2) * pow(m, 3) / 6.0 + (5.0 - 18.0 * pow(t, 2) + pow(t, 4) + 14.0 * et2 - 58.0 * et2 * pow(t, 2)) * pow(m, 5) / 120.0)return x, y#獲取基準點的下標,基準點是p[k]def get_leftbottompoint(self,p):k = 0for i in range(1, len(p)):if p[i][1] < p[k][1] or (p[i][1] == p[k][1] and p[i][0] < p[k][0]):k = ireturn k#叉乘計算方法def multiply(self, p1, p2, p0):return (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p2[0] - p0[0]) * (p1[1] - p0[1])#獲取極角，通過求反正切得出，考慮pi/2的情況def get_arc(self, p1, p0):# 兼容sort_points_tan的考慮if (p1[0] - p0[0]) == 0:if ((p1[1] - p0[1])) == 0:return -1;else:return math.pi / 2tan = float((p1[1] - p0[1])) / float((p1[0] - p0[0]))arc = math.atan(tan)if arc >= 0:return arcelse:return math.pi + arc# 對極角進行排序,排序結果list不包含基準點def sort_points_tan(self, p, pk):p2 = []for i in range(0, len(p)):p2.append({"index": i, "arc": self.get_arc(p[i], pk)})#print('排序前:',p2)p2.sort(key=lambda k: (k.get('arc')))#print('排序后:',p2)p_out = []for i in range(0, len(p2)):p_out.append(p[p2[i]["index"]])return p_out#基于逆時針遍歷多邊形,https://www.jb51.net/article/184544.htmdef convex_hull(self, p):#print(p)#p=list(set(p))#print('全部點:',p)k = self.get_leftbottompoint(p)pk = p[k]p.remove(p[k])#print('排序前去除基準點的所有點:',p,'基準點:',pk)p_sort = self.sort_points_tan(p, pk) #按與基準點連線和x軸正向的夾角排序后的點坐標#print('其余點與基準點夾角排序:',p_sort)p_result = [pk,p_sort[0]]top = 2for i in range(1, len(p_sort)):######################################叉乘為正,向前遞歸刪點;叉乘為負,序列追加新點.#最初提供的算法是下面的，計算的面積更精確，但視覺上的輪廓線不清晰，周長也不符合實際，修正后面積有增長，周長較精確，沒有十全十美的，但可以選擇更合適的。#while(self.multiply(p_result[-2], p_sort[i],p_result[-1]) > 0): while(self.multiply(p_result[-2],p_result[-1], p_sort[i]) <= 0):p_result.pop()p_result.append(p_sort[i]) return p_result#測試def HeronGetAreaOfPolyGonbyVector(self, points):# 基于海倫公式計算多邊形面積area = 0if(len(points)<3):raise Exception("error")pb=((points[-1][0]+points[0][0])/2,(points[-1][1]+points[0][1])/2) #基準點選為第一個點和最后一個點連線邊上的中點for i in range(0,len(points)-1):p1 = points[i]p2 = points[i + 1]db1 = 1 #根據維度轉化成經緯度距離d12 = 1d2b = 1#print(db1,d12,d2b)hc = (db1+d12+d2b)/2 #db1是基準點和p1的距離，d12是p1和p2的距離，d2b是p2和基準點距離#print(hc, hc-db1, hc-d12, hc-d2b)triArea = math.sqrt(hc*(hc-db1)*(hc-d12)*(hc-d2b)) #print(triArea)area += triAreareturn areadef GetAreaOfPolyGonbyVector(self, points):# 基于向量叉乘計算多邊形面積area = 0if(len(points)<3):raise Exception("error")for i in range(0,len(points)-1):p1 = points[i]p2 = points[i + 1]triArea = (p1[0]*p2[1] - p2[0]*p1[1])/2#print(triArea)area += triAreafn=(points[-1][0]*points[0][1]-points[0][0]*points[-1][1])/2#print(fn)return abs(area+fn) def GetAreaOfCircle(self, points):# 基于向量叉乘計算多邊形面積circl = 0if(len(points)<3):raise Exception("error")for i in range(0,len(points)-1):p1 = points[i]p2 = points[i + 1]cirlen = math.sqrt((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2)circl += cirlencircl += math.sqrt((points[-1][0]-points[0][0])**2 + (points[-1][1]-points[0][1])**2)print(circl)return circl def area_show(self,X):##產生聚類數據的中心qcenters = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]#按輪廓系數最小值求取作業地塊，排除簇內誤差平方和超級最大的值 for j in range(2,10):#y_predict = KMeans(n_clusters=j, random_state=9).fit_predict(X)model_kmeans=KMeans(n_clusters=j,random_state=0) #建立模型對象model_kmeans.fit(X) #訓練聚類模型y_predict = model_kmeans.predict(X)silh = silhouette_score(X.reshape(-1,2), y_predict)n_s_bigger_than_zero = (silhouette_samples(X.reshape(-1,2), y_predict) > 0).sum()sse = model_kmeans.inertia_ # 簇內誤差平方和#randindex = metrics.adjusted_rand_score(X, y_predict)#amiindex = metrics.adjusted_mutual_info_score(X[:,1], y_predict)#homoindex = metrics.homogeneity_score(X[:,1], y_predict) calinskiindex = metrics.calinski_harabasz_score(X, y_predict) # 卡賓越大越好dbiindex = metrics.davies_bouldin_score(X, y_predict) # 大衛系數越小越好print('n_cluster %d, quolity %f, num %d,cabinindex %f,dbiindex %f,sse %f' % (j,silh,n_s_bigger_than_zero, calinskiindex, dbiindex, sse)) # randindex,homoindex,homoindex#記錄前三名，最終結果在前三名中產生if((silh > qcenters[2][1]) or ((silh == qcenters[2][1]) and (sse < qcenters[2][2]))):if(silh > qcenters[1][1]):qcenters[2] = qcenters[1]if(silh > qcenters[0][1]): qcenters[1] = qcenters[0]qcenters[0] = [j,silh,sse]else:qcenters[1] = [j,silh,sse]else:qcenters[2] = [j,silh,sse]print(qcenters)pred_cluster = qcenters[0][0]if(qcenters[1][1] != 0):if(abs(1-qcenters[0][1]/qcenters[1][1])<10*abs(1-qcenters[0][2]/qcenters[1][2])): # 輪廓系數最小時，sse不能太大，這里暫時取10，純屬個人拍腦袋pred_cluster = qcenters[1][0]y_pred = KMeans(n_clusters=pred_cluster, random_state=9).fit(X) # 對x預測lables = y_pred.labels_plt.figure(1)plt.clf()plt.title('k-means:k=%d' % pred_cluster)core_samples_mask = np.zeros_like(y_pred.labels_, dtype=bool)core_samples_mask[y_pred.labels_] = Trueunique_labels = set(lables)colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))area_total = 0for k, col in zip(unique_labels, colors):##-1表示噪聲點,這里的k表示黑色if k == -1:col = 'k'break##生成一個True、False數組，lables == k 的設置成Trueclass_member_mask = (lables == k)##兩個數組做&運算，找出即是核心點又等于分類k的值 markeredgecolor='k',xy = X[class_member_mask & core_samples_mask]plt.plot(xy[:, 0], xy[:, 1], 'o', c=col,markersize=10)'''1)~優先級最高，按位對core_samples_mask 求反，求出的是噪音點的位置2)& 于運算之后，求出雖然剛開始是噪音點的位置，但是重新歸類卻屬于k的點3)對核心分類之后進行的擴展'''xy = X[class_member_mask & ~core_samples_mask] xy_arr = []for i in range(0,len(xy)):xy_arr.append([xy[i,0],xy[i,1]]) result = self.convex_hull(xy_arr)print(min(xy[:, 0]),max(xy[:, 0]),min(xy[:, 1]),max(xy[:, 1]),len(result),self.GetAreaOfCircle(result),self.GetAreaOfPolyGonbyVector(result)/667)area_total = area_total + self.GetAreaOfPolyGonbyVector(result)/667# 修正地塊邊界輪廓daw_area =np.mat(result) plt.plot(daw_area[:,0], daw_area[:,1], 'bp-', c=col, linewidth = 1)plt.plot([daw_area[0,0],daw_area[-1,0]], [daw_area[0,1],daw_area[-1,1]], 'bp-', c=col, linewidth = 1)print('作業地塊 %d，作業區域大小 %f Mu' % (pred_cluster,area_total))plt.show()if __name__ == '__main__':newarea = area_algo()##產生的簇個數n_samples=8##生產數據:此實驗結果受cluster_std的影響，或者說受eps 和cluster_std差值影響data = pd.read_excel('data/五征全天計畝.xlsx') #柴油車計畝.xlsx 0,1;14,15 index_col = '時間',nrows=10000x_oldarr = data[(data['經度']>1) & (data['滾筒轉速']>0)].iloc[:,3].valuesy_oldarr = data[(data['經度']>1) & (data['滾筒轉速']>0)].iloc[:,2].valuesX = np.zeros((len(x_oldarr),2))X[0,0],X[0,1] = newarea.LatLon2XY(y_oldarr[0],x_oldarr[0])for i in range(0, len(x_oldarr)-1):X[i+1,0],X[i+1,1] = newarea.LatLon2XY(y_oldarr[i+1],x_oldarr[i+1])newarea.area_show(X)

?

總結

以上是生活随笔為你收集整理的农机计亩区域划分问题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。