軌跡壓縮之Douglas-Peucker算法之Java實現
第一部分 問題描述
1.1 具體任務
本次作業任務是軌跡壓縮,給定一個GPS數據記錄文件,每條記錄包含經度和維度兩個坐標字段,所有記錄的經緯度坐標構成一條軌跡,要求采用合適的壓縮算法,使得壓縮后軌跡的距離誤差小于30m。
1.2 程序輸入
本程序輸入是一個GPS數據記錄文件。
1.3 數據輸出
輸出形式是文件,包括三部分,壓縮后點的ID序列及坐標、點的個數、平均距離誤差、壓縮率
第二部分 問題解答
根據問題描述,我們對問題進行求解,問題求解分為以下幾步:
2.1 數據預處理
本次程序輸入為GPS數據記錄文件,共有3150行記錄,每行記錄又分為若干個字段,根據題意,我們只需關注經度和緯度坐標字段即可,原始數據文件部分記錄如圖2.1所示:
圖2.1 原始數據文件部分記錄示意圖
如圖2.1所示,原始數據文件每條記錄中經緯度坐標字段數據的保存格式是典型的GPS坐標表達方式,即度分格式,形式為dddmm.mmmm,其中ddd表示度,mm.mmmm表示分,小數點前面表示分的整數部分,小數點后表示分的小數部分;本次數據預處理,為方便后面兩個坐標點之間距離的計算,我們需要將度分格式的經緯度坐標數據換算成度的形式,換算方法是ddd+mm.mmmm/60,此處我們保留小數點后6位數字,換算后的形式為ddd.xxxxxx。
我們以第一條記錄中經緯度坐標(11628.2491,3955.6535)為例,換算后的結果為(116.470818,39.927558),所有記錄中經緯度坐標都使用方法進行,并且可以為每一個轉換后的坐標點生成一個ID,進行唯一標識,壓縮后,我們只需輸出所有被保留點的ID即可。
2.2 Douglas-Peucker軌跡壓縮算法
軌跡壓縮算法分為兩大類,分別是無損壓縮和有損壓縮,無損壓縮算法主要包括哈夫曼編碼,有損壓縮算法又分為批處理方式和在線數據壓縮方式,其中批處理方式又包括DP(Douglas-Peucker)算法、TD-TR(Top-Down Time-Ratio)算法和Bellman算法,在線數據壓縮方式又包括滑動窗口、開放窗口、基于安全區域的方法等。
由于時間有限,本次軌跡壓縮,我們決定采用相對簡單的DP算法。
DP算法步驟如下:
(1)在軌跡曲線在曲線首尾兩點A,B之間連接一條直線AB,該直線為曲線的弦;
(2)遍歷曲線上其他所有點,求每個點到直線AB的距離,找到最大距離的點C,最大距離記為dmax;
(3)比較該距離dmax與預先定義的閾值Dmax大小,如果dmax<Dmax,則將該直線AB作為曲線段的近似,曲線段處理完畢;
(4)若dmax>=Dmax,則使C點將曲線AB分為AC和CB兩段,并分別對這兩段進行(1)~(3)步處理;
(5)當所有曲線都處理完畢時,依次連接各個分割點形成的折線,即為原始曲線的路徑。
2.3 點到直線的距離
DP算法中需要求點到直線的距離,該距離指的是垂直歐式距離,即直線AB外的點C到直線AB的距離d,此處A、B、C三點均為經緯度坐標;我們采用三角形面積相等法求距離d,具體方法是:A、B、C三點構成三角形,該三角形的面積有兩種求法,分別是普通方法(底x高/2)和海倫公式,海倫公式如下:
假設有一個三角形,邊長分別為a、b、c,三角形的面積S可由以下公式求得:,其中p為半周長:。
我們通過海倫公式求得三角形面積,然后就可以求得高的大小,此處高即為距離d。要想用海倫公式,必須求出A、B、C三點兩兩之間的距離,該距離公式是由老師給出的,直接調用距離函數即可。
注意:求出距離后,要加上絕對值,以防止距離為負數。
2.4 平均誤差求解
平均誤差指的是壓縮時忽略的那些點到對應線段的距離之和除以總點數得到的數值。
2.5 壓縮率求解
壓縮率的計算公式如下:
2.6 數據結果文件的生成
經過上面的處理和計算,我們將壓縮后剩余點的ID和點的個數、平均距離誤差、壓縮率等參數都寫入最終的結果文件中,問題解答完成。
第三部分 代碼實現
本程序采用Java語言編寫,開發環境為IntelliJ IDEA 14.0.2,代碼共分為兩個類,一個是ENPoint類,用于保存經緯度點信息,一個是TrajectoryCompressionMain類,用于編寫數據處理、DP算法、點到直線距離、求平均誤差等函數。
3.1 程序總流程
整個程序流程主要包括以下幾個步驟:
(1)定義相關ArrayList數組和File對象,其中ArrayList數組對象有三個,分別是原始經緯度坐標數組pGPSArryInit、過濾后的點坐標數組pGPSArrayFilter、過濾并排序后的點坐標數組pGPSArrayFilterSort;File文件對象共有五個,分別是原始數據文件對象fGPS、壓縮后的結果數據文件對象oGPS、保持轉換后的原始經緯度坐標點的數據文件fInitGPSPoint、仿真測試文件fTestInitPoint和fTestFilterPoint。
(2)獲取原始點坐標并將其寫入到文件中,主要包括讀文件和寫文件兩種操作;
(3)進行軌跡壓縮;
(4)對壓縮后的經緯度點坐標進行排序;
(5)生成仿真測試文件,并用R語言工具進行圖形繪制,得到最終的結果;
(6)求平均誤差和壓縮率,平均誤差通過函數求得,壓縮率直接計算獲得;
(7)將最終結果寫入結果文件中,包括過濾后的點的ID,點的個數、平均誤差和壓縮率;
3.2 具體實現代碼
(1)ENPoint.java
1 package cc.xidian.main; 2 3 import java.text.DecimalFormat; 4 5 /** 6 * Created by hadoop on 2015/12/20. 7 */ 8 public class ENPoint implements Comparable<ENPoint>{ 9 public int id;//點ID 10 public double pe;//經度 11 public double pn;//維度 12 13 public ENPoint(){}//空構造函數 14 public String toString(){ 15 //DecimalFormat df = new DecimalFormat("0.000000"); 16 return this.id+"#"+this.pn+","+this.pe; 17 } 18 public String getTestString(){ 19 DecimalFormat df = new DecimalFormat("0.000000"); 20 return df.format(this.pn)+","+df.format(this.pe); 21 } 22 public String getResultString(){ 23 DecimalFormat df = new DecimalFormat("0.000000"); 24 return this.id+"#"+df.format(this.pn)+","+df.format(this.pe); 25 } 26 @Override 27 public int compareTo(ENPoint other) { 28 if(this.id<other.id) return -1; 29 else if(this.id>other.id) return 1; 30 else 31 return 0; 32 } 33 }
(2)TrajectoryCompressionMain.java
1 package cc.xidian.main; 2 3 import java.io.*; 4 import java.text.DecimalFormat; 5 import java.util.*; 6 import java.util.List; 7 8 /** 9 * Created by hadoop on 2015/12/19. 10 */ 11 public class TrajectoryCompressionMain{ 12 public static void main(String[] args)throws Exception{ 13 //-----------------------1、相關ArrayList數組和File對象的聲明和定義-------------------------------------------------// 14 ArrayList<ENPoint> pGPSArrayInit = new ArrayList<ENPoint>();//原紀錄經緯度坐標數組 15 ArrayList<ENPoint> pGPSArrayFilter = new ArrayList<ENPoint>();//過濾后的經緯度坐標數組 16 ArrayList<ENPoint> pGPSArrayFilterSort = new ArrayList<ENPoint>();//過濾并排序后的經緯度坐標數組 17 File fGPS = new File("2007-10-14-GPS.log");//原始數據文件對象 18 File oGPS = new File("2015-12-25-GPS-Result.log");//過濾后的結果數據文件對象 19 //保持轉換成度后的原始經緯度數據文件,保持格式為“ID#經緯值,緯度值”,其中經度和維度單位為度,并保留小數點后6位數字 20 File fInitGPSPoint = new File("2007-10-14-GPS-ENPoint.log");//保持轉換后的原始經緯度坐標點的數據文件 21 File fTestInitPoint = new File("2007-10-14-GPS-InitTestPoint.log");//用于仿真的原始經緯度坐標點數據文件 22 File fTestFilterPoint = new File("2015-12-25-GPS-FilterTestPoint.log");//用于仿真的過濾后的經緯度坐標點數據文件 23 //-------------------------2、獲取原始點坐標并將其寫入到文件中-------------------------------------------------------// 24 pGPSArrayInit = getENPointFromFile(fGPS);//從原始數據文件中獲取轉換后的經緯度坐標點數據,存放到ArrayList數組中 25 writeInitPointToFile(fInitGPSPoint, pGPSArrayInit);//將轉換后的原始經緯度點數據寫入文件中 26 System.out.println(pGPSArrayInit.size());//輸出原始經緯度點坐標的個數 27 //-------------------------3、進行軌跡壓縮-----------------------------------------------------------------------// 28 double DMax = 30.0;//設定最大距離誤差閾值 29 pGPSArrayFilter.add(pGPSArrayInit.get(0));//獲取第一個原始經緯度點坐標并添加到過濾后的數組中 30 pGPSArrayFilter.add(pGPSArrayInit.get(pGPSArrayInit.size()-1));//獲取最后一個原始經緯度點坐標并添加到過濾后的數組中 31 ENPoint[] enpInit = new ENPoint[pGPSArrayInit.size()];//使用一個點數組接收所有的點坐標,用于后面的壓縮 32 Iterator<ENPoint> iInit = pGPSArrayInit.iterator(); 33 int jj=0; 34 while(iInit.hasNext()){ 35 enpInit[jj] = iInit.next(); 36 jj++; 37 }//將ArrayList中的點坐標拷貝到點數組中 38 int start = 0;//起始下標 39 int end = pGPSArrayInit.size()-1;//結束下標 40 TrajCompressC(enpInit,pGPSArrayFilter,start,end,DMax);//DP壓縮算法 41 System.out.println(pGPSArrayFilter.size());//輸出壓縮后的點數 42 //-------------------------4、對壓縮后的經緯度點坐標數據按照ID從小到大排序---------------------------------------------// 43 ENPoint[] enpFilter = new ENPoint[pGPSArrayFilter.size()];//使用一個點數組接收過濾后的點坐標,用于后面的排序 44 Iterator<ENPoint> iF = pGPSArrayFilter.iterator(); 45 int i = 0; 46 while(iF.hasNext()){ 47 enpFilter[i] = iF.next(); 48 i++; 49 }//將ArrayList中的點坐標拷貝到點數組中 50 Arrays.sort(enpFilter);//進行排序 51 for(int j=0;j<enpFilter.length;j++){ 52 pGPSArrayFilterSort.add(enpFilter[j]);//將排序后的點坐標寫到一個新的ArrayList數組中 53 } 54 //-------------------------5、生成仿真測試文件--------------------------------------------------------------------// 55 writeTestPointToFile(fTestInitPoint,pGPSArrayInit);//將原始經緯度數據點寫入仿真文件中,格式為“經度,維度” 56 writeTestPointToFile(fTestFilterPoint, pGPSArrayFilterSort);//將過濾后的經緯度數據點寫入仿真文件中,格式為“經度,維度” 57 //-------------------------6、求平均誤差-------------------------------------------------------------------------// 58 double mDError = getMeanDistError(pGPSArrayInit,pGPSArrayFilterSort);//求平均誤差 59 System.out.println(mDError); 60 //-------------------------7、求壓縮率--------------------------------------------------------------------------// 61 double cRate = (double)pGPSArrayFilter.size()/pGPSArrayInit.size()*100;//求壓縮率 62 System.out.println(cRate); 63 //-------------------------8、生成最終結果文件--------------------------------------------------------------------// 64 //將最終結果寫入結果文件中,包括過濾后的點的ID,點的個數、平均誤差和壓縮率 65 writeFilterPointToFile(oGPS,pGPSArrayFilterSort,mDError,cRate); 66 //------------------------------------------------------------------------------------------------------------// 67 } 68 69 /** 70 *函數功能:從源文件中讀出所以記錄中的經緯度坐標,并存入到ArrayList數組中,并將其返回 71 * @param fGPS:源數據文件 72 * @return pGPSArrayInit:返回保存所有點坐標的ArrayList數組 73 * @throws Exception 74 */ 75 public static ArrayList<ENPoint> getENPointFromFile(File fGPS)throws Exception{ 76 ArrayList<ENPoint> pGPSArray = new ArrayList<ENPoint>(); 77 if(fGPS.exists()&&fGPS.isFile()){ 78 InputStreamReader read = new InputStreamReader(new FileInputStream(fGPS)); 79 BufferedReader bReader = new BufferedReader(read); 80 String str; 81 String[] strGPS; 82 int i = 0; 83 while((str = bReader.readLine())!=null){ 84 strGPS = str.split(" "); 85 ENPoint p = new ENPoint(); 86 p.id = i; 87 i++; 88 p.pe = (dfTodu(strGPS[3])); 89 p.pn = (dfTodu(strGPS[5])); 90 pGPSArray.add(p); 91 } 92 bReader.close(); 93 } 94 return pGPSArray; 95 } 96 97 /** 98 * 函數功能:將過濾后的點的經緯度坐標、平均距離誤差、壓縮率寫到結果文件中 99 * @param outGPSFile:結果文件 100 * @param pGPSPointFilter:過濾后的點 101 * @param mDerror:平均距離誤差 102 * @param cRate:壓縮率 103 * @throws Exception 104 */ 105 public static void writeFilterPointToFile(File outGPSFile,ArrayList<ENPoint> pGPSPointFilter, 106 double mDerror,double cRate)throws Exception{ 107 Iterator<ENPoint> iFilter = pGPSPointFilter.iterator(); 108 RandomAccessFile rFilter = new RandomAccessFile(outGPSFile,"rw"); 109 while(iFilter.hasNext()){ 110 ENPoint p = iFilter.next(); 111 String sFilter = p.getResultString()+"\n"; 112 byte[] bFilter = sFilter.getBytes(); 113 rFilter.write(bFilter); 114 } 115 String strmc = "#"+Integer.toString(pGPSPointFilter.size())+","+ 116 Double.toString(mDerror)+","+Double.toString(cRate)+"%"+"#"+"\n"; 117 byte[] bmc = strmc.getBytes(); 118 rFilter.write(bmc); 119 120 rFilter.close(); 121 } 122 /** 123 * 函數功能:將轉換后的原始經緯度數據點存到文件中 124 * @param outGPSFile 125 * @param pGPSPointFilter 126 * @throws Exception 127 */ 128 public static void writeInitPointToFile(File outGPSFile,ArrayList<ENPoint> pGPSPointFilter)throws Exception{ 129 Iterator<ENPoint> iFilter = pGPSPointFilter.iterator(); 130 RandomAccessFile rFilter = new RandomAccessFile(outGPSFile,"rw"); 131 while(iFilter.hasNext()){ 132 ENPoint p = iFilter.next(); 133 String sFilter = p.toString()+"\n"; 134 byte[] bFilter = sFilter.getBytes(); 135 rFilter.write(bFilter); 136 } 137 rFilter.close(); 138 } 139 /** 140 * 函數功能:將數組中的經緯度點坐標數據寫入測試文件中,用于可視化測試 141 * @param outGPSFile:文件對象 142 * @param pGPSPointFilter:點數組 143 * @throws Exception 144 */ 145 public static void writeTestPointToFile(File outGPSFile,ArrayList<ENPoint> pGPSPointFilter)throws Exception{ 146 Iterator<ENPoint> iFilter = pGPSPointFilter.iterator(); 147 RandomAccessFile rFilter = new RandomAccessFile(outGPSFile,"rw"); 148 while(iFilter.hasNext()){ 149 ENPoint p = iFilter.next(); 150 String sFilter = p.getTestString()+"\n"; 151 byte[] bFilter = sFilter.getBytes(); 152 rFilter.write(bFilter); 153 } 154 rFilter.close(); 155 } 156 157 /** 158 * 函數功能:將原始經緯度坐標數據轉換成度 159 * @param str:原始經緯度坐標 160 * @return :返回對于的度數據 161 */ 162 public static double dfTodu(String str){ 163 int indexD = str.indexOf('.'); 164 String strM = str.substring(0,indexD-2); 165 String strN = str.substring(indexD-2); 166 double d = Double.parseDouble(strM)+Double.parseDouble(strN)/60; 167 return d; 168 } 169 /** 170 * 函數功能:保留一個double數的小數點后六位 171 * @param d:原始double數 172 * @return 返回轉換后的double數 173 */ 174 public static double getPointSix(double d){ 175 DecimalFormat df = new DecimalFormat("0.000000"); 176 return Double.parseDouble(df.format(d)); 177 } 178 /** 179 * 函數功能:使用三角形面積(使用海倫公式求得)相等方法計算點pX到點pA和pB所確定的直線的距離 180 * @param pA:起始點 181 * @param pB:結束點 182 * @param pX:第三個點 183 * @return distance:點pX到pA和pB所在直線的距離 184 */ 185 public static double distToSegment(ENPoint pA,ENPoint pB,ENPoint pX){ 186 double a = Math.abs(geoDist(pA, pB)); 187 double b = Math.abs(geoDist(pA, pX)); 188 double c = Math.abs(geoDist(pB, pX)); 189 double p = (a+b+c)/2.0; 190 double s = Math.sqrt(Math.abs(p*(p-a)*(p-b)*(p-c))); 191 double d = s*2.0/a; 192 return d; 193 } 194 195 /** 196 * 函數功能:用老師給的看不懂的方法求兩個經緯度點之間的距離 197 * @param pA:起始點 198 * @param pB:結束點 199 * @return distance:距離 200 */ 201 public static double geoDist(ENPoint pA,ENPoint pB) 202 { 203 double radLat1 = Rad(pA.pn); 204 double radLat2 = Rad(pB.pn); 205 double delta_lon = Rad(pB.pe - pA.pe); 206 double top_1 = Math.cos(radLat2) * Math.sin(delta_lon); 207 double top_2 = Math.cos(radLat1) * Math.sin(radLat2) - Math.sin(radLat1) * Math.cos(radLat2) * Math.cos(delta_lon); 208 double top = Math.sqrt(top_1 * top_1 + top_2 * top_2); 209 double bottom = Math.sin(radLat1) * Math.sin(radLat2) + Math.cos(radLat1) * Math.cos(radLat2) * Math.cos(delta_lon); 210 double delta_sigma = Math.atan2(top, bottom); 211 double distance = delta_sigma * 6378137.0; 212 return distance; 213 } 214 /** 215 * 函數功能:角度轉弧度 216 * @param d:角度 217 * @return 返回的是弧度 218 */ 219 public static double Rad(double d) 220 { 221 return d * Math.PI / 180.0; 222 } 223 224 /** 225 * 函數功能:根據最大距離限制,采用DP方法遞歸的對原始軌跡進行采樣,得到壓縮后的軌跡 226 * @param enpInit:原始經緯度坐標點數組 227 * @param enpArrayFilter:保持過濾后的點坐標數組 228 * @param start:起始下標 229 * @param end:終點下標 230 * @param DMax:預先指定好的最大距離誤差 231 */ 232 public static void TrajCompressC(ENPoint[] enpInit,ArrayList<ENPoint> enpArrayFilter, 233 int start,int end,double DMax){ 234 if(start < end){//遞歸進行的條件 235 double maxDist = 0;//最大距離 236 int cur_pt = 0;//當前下標 237 for(int i=start+1;i<end;i++){ 238 double curDist = distToSegment(enpInit[start],enpInit[end],enpInit[i]);//當前點到對應線段的距離 239 if(curDist > maxDist){ 240 maxDist = curDist; 241 cur_pt = i; 242 }//求出最大距離及最大距離對應點的下標 243 } 244 //若當前最大距離大于最大距離誤差 245 if(maxDist >= DMax){ 246 enpArrayFilter.add(enpInit[cur_pt]);//將當前點加入到過濾數組中 247 //將原來的線段以當前點為中心拆成兩段,分別進行遞歸處理 248 TrajCompressC(enpInit,enpArrayFilter,start,cur_pt,DMax); 249 TrajCompressC(enpInit,enpArrayFilter,cur_pt,end,DMax); 250 } 251 } 252 } 253 /** 254 * 函數功能:求平均距離誤差 255 * @param pGPSArrayInit:原始數據點坐標 256 * @param pGPSArrayFilterSort:過濾后的數據點坐標 257 * @return :返回平均距離 258 */ 259 public static double getMeanDistError( 260 ArrayList<ENPoint> pGPSArrayInit,ArrayList<ENPoint> pGPSArrayFilterSort){ 261 double sumDist = 0.0; 262 for(int i=1;i<pGPSArrayFilterSort.size();i++){ 263 int start = pGPSArrayFilterSort.get(i-1).id; 264 int end = pGPSArrayFilterSort.get(i).id; 265 for(int j=start+1;j<end;j++){ 266 sumDist += distToSegment( 267 pGPSArrayInit.get(start),pGPSArrayInit.get(end),pGPSArrayInit.get(j)); 268 } 269 } 270 double meanDist = sumDist/(pGPSArrayInit.size()); 271 return meanDist; 272 } 273 }
第四部分 程序結果
4.1 程序輸出結果
壓縮后的結果:
(1)總點數:140個點;(2)平均距離誤差:7.943786;(3)壓縮率:4.4444%
4.2 仿真結果
經過軌跡壓縮,我們將原始經緯度坐標點轉換為壓縮過濾后的經緯度坐標點,我們將這兩種點坐標數據分別寫入兩個文件中,然后根據這兩個文件使用R語言進行圖形繪制,分別畫出壓縮前和壓縮后的軌跡 ,進行對比,根據對比結果就可以看出我們軌跡壓縮算法是否有效,最終對比結果如圖4.1所示:
第五部分 總結
本程序編寫過程中,遇到了各種各樣的問題,也學到了很多編程經驗,下面就遇到的問題及解決方案做一個總結,最后對程序存在的不足提出改進建議。
5.1 遇到的問題及解決方案
問題1:經緯度坐標順序問題
解決:距離公式中的參數是緯度在前經度在后,需要調整下經緯度坐標點的順序。
問題2:距離不能為負值
解決:保證求出的距離不能為負值,加絕對值函數即可。
問題3:DP算法實現細節
解決:開始使用ArrayList數組解決下標問題,遞歸求解時出現巨大誤差,后來改用普通數組下標進行遞歸,結果好多了。
5.2 存在的不足與展望
(1)進行軌跡壓縮時,DP算法是最簡單的一種算法,并不是最優的,可選用一些效果好的算法再次進行軌跡壓縮;
(2)本次實驗數據的記錄為3150條,數據量不算大,如果有10億條數據,該怎么辦呢?我們可以從硬件、分布式、數據預處理、數據切分、性能好的數據倉庫等方面考慮。
注:原始數據文件見360云盤鏈接:https://yunpan.cn/cu8BiNsvH66SD 訪問密碼 7be0