CN105355042A - Road network extraction method based on taxi GPS - Google Patents

Abstract

The invention discloses a road network extraction method based on a taxi GPS, and the method comprises seven steps: S1, extracting a taxi GPS database of one workday; S2, limiting and recognizing an city boundary; S3, carrying out data cleaning; S4, dividing an observation range; S5, judging an intersection; S6, carrying out road segment recognition based on the intersection; S7, outputting a road attribute list. The method is newer and more accurate in extracted road network information, is high in calculation speed, reduces the manpower cost, and reduces the errors.

Description

A kind of road network extraction method based on taxi GPS

Technical field

The present invention relates to road network extraction method, particularly a kind of road network extraction method based on taxi GPS.

Background technology

Urban Transportation Development in China is swift and violent, and road construction is maked rapid progress.Domestic still shortage, can be issued for the authoritative road network data of traffic analysis for traffic planning designer or traffic study person, therefore, the road network information of the previous exploration that domestic communication planning or the present situation road network studied often adopt planning department to provide, or planning team spends a large amount of manpower and materials to investigate to upgrade and rectification road network information on the spot.The road network information that the former obtains is seriously delayed; Latter also will spend a large amount of time, and the data deficiency of acquisition is ageing.

Summary of the invention

Goal of the invention: the object of this invention is to provide a kind of low cost, the high efficiency road network extraction method based on taxi GPS.

Road network extraction method based on taxi GPS of the present invention, comprises following step:

S1: extract a workaday GPS data from taxi storehouse: database comprises 8 row, i.e. car number, trigger event, operation state, gps time, GPS longitude, GPS latitude, GPS speed and GPS direction;

S2: limit and identify urban border: extract westernmost end, border, municipal administration region, most the east, northernmost with node coordinate southernmost, node coordinate comprises longitude and latitude; The longitude coordinate arranging westernmost end is the longitude minimum value LON of urban border _min, the longitude coordinate that arranges most the east is the longitude maximal value LON of urban border _max, the latitude maximal value LAT that latitude coordinate is northernmost urban border is set _max, the latitude coordinate arranged southernmost is the latitude minimum value LAT of urban border _min, the region that four boundary values enclose formation is observation scope; Reject the data line not meeting formula (1) in the GPS data from taxi storehouse obtained in step S1;

LON _min<GPS longitude <LON _max, and LAT _min<GPS latitude <LAT _max(1)

S3: data cleansing: by the GPS data from taxi storehouse after step S2 process successively according to gps time, the arrangement of car number Two Variables ascending order, make the data line of same vehicle in the new database formed connect together and in chronological sequence arrange; Delete non-rice habitats data line and malposition data line, obtain effective GPS data from taxi storehouse; Setting global iterative variable i teration is 1; Wherein, the data line that non-rice habitats data line is cancellation mark box-like (2) is deleted:

Rule of judgment 1: operation state=2 or 3, and GPS speed=0 (2)

Delete position abnormal data is capable is delete D _{x, x+1}meet (x+1)th row data time formula (3):

Rule of judgment 2:D _{x, x+1}>V _max(T _x+1-T _x) (3)

D _{x, x+1}shown in (4):

$D_{x,x+1}=2r\&CenterDot;\arcsin \left(\sqrt{{\sin }^2\left(\frac{{\mathrm{\&phi;}}_{x+1}-{\mathrm{\&phi;}}_x}{2}\right)+cos\left({\mathrm{\&phi;}}_x\right)cos\left({\mathrm{\&phi;}}_{x+1}\right){\sin }^2\left(\frac{{\mathrm{\&lambda;}}_{x+1}-{\mathrm{\&lambda;}}_x}{2}\right)}\right)---\left(4\right)$

Wherein, D _{x, x+1}for the distance between adjacent lines two data points in the database that obtains after deleting non-rice habitats data line, V _maxfor the maximum speed limit of urban road, r is the radius of the earth, T _x, φ _x, λ _xfor xth row gps time, GPS latitude and GPS longitude in the database that obtains after deleting non-rice habitats data line, T _x+1, φ _x+1, λ _x+1for (x+1)th row gps time, GPS latitude and GPS longitude in the database that obtains after deleting non-rice habitats data line;

S4: the division of observation scope: the observation scope determined by step S2 is cut into the rectangle plot of the capable n row of m, and wherein m and n is obtained by formula (5) and formula (6), the latitude lat of cut-boundary _{limit i}with longitude lon _{limit j}obtained by formula (7) and formula (8) respectively:

m＝INT((LAT _max-LAT _min)/0.0027)(5)

n＝INT((LON _max-LON _min)/0.003528)(6)

Lat _{limit i}=LAT _min+ (i-0.5*iteration) * (LAT _max-LAT _min)/m, (i=2,3,4 ..., m) (7)

Lon _{limit j}=LON _min+ (j-0.5*iteration) * (LON _max-LON _min)/n, (j=2,3,4 ..., n) (8)

S5: the differentiation of crossing: crossing identification is carried out to the rectangle plot obtained in step S4, also all latitude and longitude coordinates composition subdata bases dropping on the effective GPS data from taxi storehouse in a plot are about to, to all subdata bases, according to the span of longitude and latitude, evenly be divided into 10 × 10 detecting units, then judge according to crossing Rule of judgment the position that whether there is crossing and crossing in each detecting unit, the position of the crossing of all identification is recorded in the detection list of crossing, finally assignment and differentiation are carried out to global iterative variable i teration, wherein, crossing Rule of judgment is: the data line comprised in detecting unit is more than 40 row, and the number of unduplicated " GPS orientation " value that these data lines comprise is greater than 40,

S6: the section based on crossing identifies: the crossing detection list obtained by step S5 is converted into intersection stomion his-and-hers watches to be detected, each intersection stomion to be detected forms by initial crossing and termination crossing; Intersection stomion pair to be detected is organized to each, judges initial crossing and stop whether there is road between crossing, if exist, then by initial crossing with stop the coordinate record of crossing in the detection list of section, and delete and repeat to judge section by accident;

S7: export road attribute list: the crossing detection list that output step S5 obtains and the section detection list that step S6 obtains.

Further, the position of the crossing in described step S5 comprises crossing longitude and crossing latitude, and crossing longitude is the mean value of the longitude of all data lines in detecting unit, and crossing latitude is the mean value of the latitude of all data lines in detecting unit.

Further, carrying out assignment and differentiation to global iterative variable i teration in described step S5 is add 1 to global iterative variable i teration and assignment again, if judge that global iterative variable i teration is not equal to 3, then get back to step S4, if global iterative variable i teration equals 3, then enter step S6.

Further, described step S6 comprises following sub-step:

S6.1: crossing detection list is converted into intersection stomion his-and-hers watches to be detected: according to the distance between the arbitrary crossing in the detection list of formula (4) calculating crossing and all the other crossings, filter out this distance and be less than the intersection stomion of 2km to composition intersection stomion his-and-hers watches to be detected, these his-and-hers watches comprise 5 variablees, i.e. initial crossing longitude lon _start, initial crossing latitude lat _start, stop crossing longitude lon _end, stop crossing latitude lat _end, two crossing distance D _s-e; Wherein, calculate according to formula (9), " s " represents initial crossing, and " e " representative stops crossing:

$D_{s-e}=2r\&CenterDot;\arcsin \left(\sqrt{{\sin }^2\left(\frac{{\mathrm{\&phi;}}_e-{\mathrm{\&phi;}}_s}{2}\right)+cos\left({\mathrm{\&phi;}}_s\right)cos\left({\mathrm{\&phi;}}_e\right){\sin }^2\left(\frac{{\mathrm{\&lambda;}}_e-{\mathrm{\&lambda;}}_s}{2}\right)}\right)---\left(9\right)$

S6.2: intersection stomion pair to be detected is organized to each, judge initial crossing and stop whether there is road between crossing: first, limit Road Detection region, this region by initial crossing with stop crossing coordinate line and respectively to offset two straight lines of 20m to both sides and enclose form respectively through initial crossing, the vertical line that stops these two straight lines of crossing; Then surveyed area is divided into along crossing line orientation average individual rectangular sub blocks, does not contain any effective GPS data from taxi storehouse latitude and longitude coordinates point if be less than in 20 Continuous Rectangular sub-blocks, then shows to there is road between two crossings, be recorded in Road Detection list by 5 right for this intersection stomion variablees;

S6.3: delete and repeat to judge section by accident: according to the angle θ stopping crossing and initial crossing line in the Road Detection list that formula (10) calculation procedure S6.2 obtains and formed, wherein quadrant is included in the arctan function considered by atan2; By Road Detection list successively according to angle θ and initial crossing longitude lon _startsequence, makes the data line of identical initial crossing connect together and sorts according to angle value size; Contrast the angle θ that the adjacent data of identical initial crossing is capable, if difference is less than 40 °, then delete larger distance D _s-edata line; Then the value of the angle θ in all data lines between 270 ° and 360 ° is replaced with θ-360 °, again by Road Detection list successively according to angle θ and initial crossing longitude lon _startsequence, contrasts the angle θ that the adjacent data of identical initial crossing is capable, if difference is less than 40 °, then deletes larger distance D _s-edata line;

θ＝atan2(lat _end-lat _start,lon _end-lon _start)(10)。

Beneficial effect: compared with prior art, the present invention has following beneficial effect:

1) road network information that extracts upgrades, more accurate: tradition extracts the method for road information by the present situation city CAD road network that planning department provides, and relies on the update time of planning department road network to a great extent; China's most cities planning department road network upgrades often synchronous with city planning, within several years, carries out once, and the road network information of extraction certainly will be caused delayed; And the inventive method can be extracted according to recent GPS data from taxi, as long as there is the road of taxi process all can be identified, substantially increase the ageing of road network information and accuracy;

2) fast operation: traditional road network extraction method, on delayed CAD basis, present situation city, make an on-the-spot survey on the spot in conjunction with a large amount of manpowers and could obtain up-to-date road network information, whole process takes time and effort; Fact-finding process can continue the even longer time in a few week; And the inventive method utilizes recent GPS data from taxi, detailed up-to-date road network information can be extracted in several hours, greatly save man-hour, improve the efficiency of traffic Correlative plan;

3) the automatic computing of computing machine, decrease human input, reduce error: this inventive method can be obtained by Computing completely, has both saved the input of manpower survey, prevent again manpower to participate in and the personal error of generation, also greatly save the cost of Data Update.

Accompanying drawing explanation

Fig. 1 is method flow diagram of the present invention;

Fig. 2 is the schematic diagram that observation scope of the present invention divides;

Fig. 3 is the schematic diagram in Road Detection region of the present invention.

Embodiment

Below in conjunction with embodiment, technical scheme of the present invention is further introduced.

A kind of road network extraction method based on taxi GPS provided by the invention, as shown in Figure 1, comprises following step:

S1: extract a workaday GPS data from taxi storehouse: database comprises 8 row, i.e. car number, trigger event, operation state, gps time, GPS longitude, GPS latitude, GPS speed and GPS direction, the extraction call format of variable and implication are such as formula shown in (1):

S2: limit and identify urban border: extract westernmost end, border, municipal administration region, most the east, northernmost with node coordinate southernmost, node coordinate comprises longitude and latitude; The longitude coordinate arranging westernmost end is the longitude minimum value LON of urban border _min, the longitude coordinate that arranges most the east is the longitude maximal value LON of urban border _max, the latitude maximal value LAT that latitude coordinate is northernmost urban border is set _max, the latitude coordinate arranged southernmost is the latitude minimum value LAT of urban border _min, the region that four boundary values enclose formation is observation scope; Reject the data line not meeting formula (2) in the GPS data from taxi storehouse obtained in step S1;

LON _min<GPS longitude <LON _max, and LAT _min<GPS latitude <LAT _max(2)

S3: data cleansing: by the GPS data from taxi storehouse after step S2 process successively according to gps time, the arrangement of car number Two Variables ascending order, make the data line of same vehicle in the new database formed connect together and in chronological sequence arrange; Delete non-rice habitats data line and malposition data line, obtain effective GPS data from taxi storehouse; Setting global iterative variable i teration is 1; Wherein, the data line that non-rice habitats data line is cancellation mark box-like (3) is deleted:

Rule of judgment 1: operation state=2 or 3, and GPS speed=0 (3) delete position abnormal data capable be delete D _{x, x+1}meet (x+1)th row data time formula (4):

Rule of judgment 2:D _{x, x+1}>V _max(T _x+1-T _x) (4) D _{x, x+1}shown in (5):

$D_{x,x+1}=2r\&CenterDot;\arcsin \left(\sqrt{{\sin }^2\left(\frac{{\mathrm{\&phi;}}_{x+1}-{\mathrm{\&phi;}}_x}{2}\right)+cos\left({\mathrm{\&phi;}}_x\right)cos\left({\mathrm{\&phi;}}_{x+1}\right){\sin }^2\left(\frac{{\mathrm{\&lambda;}}_{x+1}-{\mathrm{\&lambda;}}_x}{2}\right)}\right)---\left(5\right)$

Wherein, D _{x, x+1}for the distance between adjacent lines two data points in the database that obtains after deleting non-rice habitats data line, V _maxfor the maximum speed limit of urban road, r is the radius of the earth, T _x, φ _x, λ _xfor xth row gps time, GPS latitude and GPS longitude in the database that obtains after deleting non-rice habitats data line, T _x+1, φ _x+1, λ _x+1for (x+1)th row gps time, GPS latitude and GPS longitude in the database that obtains after deleting non-rice habitats data line;

S4: the division of observation scope: as shown in Figure 2, the observation scope determined by step S2 is cut into the rectangle plot of the capable n row of m, and wherein m and n is obtained by formula (6) and formula (7), the latitude lat of cut-boundary _{limit i}with longitude lon _{limit j}obtained by formula (8) and formula (9) respectively:

m＝INT((LAT _max-LAT _min)/0.0027)(6)

n＝INT((LON _max-LON _min)/0.003528)(7)

Lat _{limit i}=LAT _min+ (i-0.5*iteration) * (LAT _max-LAT _min)/m, (i=2,3,4 ..., m) (8)

Lon _{limit j}=LON _min+ (j-0.5*iteration) * (LON _max-LON _min)/n, (j=2,3,4 ..., n) (9)

S5: the differentiation of crossing: crossing identification is carried out to the rectangle plot obtained in step S4, also all latitude and longitude coordinates composition subdata bases dropping on the effective GPS data from taxi storehouse in a plot are about to, to all subdata bases, according to the span of longitude and latitude, evenly be divided into 10 × 10 detecting units, then judge according to crossing Rule of judgment the position that whether there is crossing and crossing in each detecting unit, the position of the crossing of all identification is recorded in the detection list of crossing, finally assignment and differentiation are carried out to global iterative variable i teration, wherein, crossing Rule of judgment is: the data line comprised in detecting unit is more than 40 row, and the number of unduplicated " GPS orientation " value that these data lines comprise is greater than 40, the position of crossing comprises crossing longitude and crossing latitude, and crossing longitude is the mean value of the longitude of all data lines in detecting unit, and crossing latitude is the mean value of the latitude of all data lines in detecting unit, carrying out assignment and differentiation to global iterative variable i teration is add 1 to global iterative variable i teration and assignment again, if judge that global iterative variable i teration is not equal to 3, then get back to step S4, if global iterative variable i teration equals 3, then enter step S6,

S6: the section based on crossing identifies: the crossing detection list obtained by step S5 is converted into intersection stomion his-and-hers watches to be detected, each intersection stomion to be detected forms by initial crossing and termination crossing; Intersection stomion pair to be detected is organized to each, judges initial crossing and stop whether there is road between crossing, if exist, then by initial crossing with stop the coordinate record of crossing in the detection list of section, and delete and repeat to judge section by accident;

Step S6 comprises following sub-step:

S6.1: crossing detection list is converted into intersection stomion his-and-hers watches to be detected: according to the distance between the arbitrary crossing in the detection list of formula (5) calculating crossing and all the other crossings, filter out this distance and be less than the intersection stomion of 2km to composition intersection stomion his-and-hers watches to be detected, these his-and-hers watches comprise 5 variablees, i.e. initial crossing longitude lon _start, initial crossing latitude lat _start, stop crossing longitude lon _end, stop crossing latitude lat _end, two crossing distance D _s-e; Wherein, calculate according to formula (10), " s " represents initial crossing, and " e " representative stops crossing:

$D_{s-e}=2r\&CenterDot;\arcsin \left(\sqrt{{\sin }^2\left(\frac{{\mathrm{\&phi;}}_e-{\mathrm{\&phi;}}_s}{2}\right)+cos\left({\mathrm{\&phi;}}_s\right)cos\left({\mathrm{\&phi;}}_e\right){\sin }^2\left(\frac{{\mathrm{\&lambda;}}_e-{\mathrm{\&lambda;}}_s}{2}\right)}\right)---\left(10\right)$

S6.2: intersection stomion pair to be detected is organized to each, judge initial crossing and stop whether there is road between crossing: first, limit Road Detection region, as shown in Figure 3, this region by initial crossing with stop crossing coordinate line and respectively to offset two straight lines of 20m to both sides and enclose form respectively through initial crossing, the vertical line that stops these two straight lines of crossing; Then surveyed area is divided into along crossing line orientation average individual rectangular sub blocks, does not contain any effective GPS data from taxi storehouse latitude and longitude coordinates point if be less than in 20 Continuous Rectangular sub-blocks, then shows to there is road between two crossings, be recorded in Road Detection list by 5 right for this intersection stomion variablees;

S6.3: delete and repeat to judge section by accident: according to the angle θ stopping crossing and initial crossing line in the Road Detection list that formula (11) calculation procedure S6.2 obtains and formed, wherein quadrant is included in the arctan function considered by atan2; By Road Detection list successively according to angle θ and initial crossing longitude lon _startsequence, makes the data line of identical initial crossing connect together and sorts according to angle value size; Contrast the angle θ that the adjacent data of identical initial crossing is capable, if difference is less than 40 °, then delete larger distance D _s-edata line; Then the value of the angle θ in all data lines between 270 ° and 360 ° is replaced with θ-360 °, again by Road Detection list successively according to angle θ and initial crossing longitude lon _startsequence, contrasts the angle θ that the adjacent data of identical initial crossing is capable, if difference is less than 40 °, then deletes larger distance D _s-edata line;

θ＝atan2(lat _end-lat _start,lon _end-lon _start)(11)。

S7: export road attribute list: the crossing detection list that output step S5 obtains and the section detection list that step S6 obtains.

Claims (4)

1. based on a road network extraction method of taxi GPS, it is characterized in that: comprise following step:

S1: extract a workaday GPS data from taxi storehouse: database comprises 8 row, i.e. car number, trigger event, operation state, gps time, GPS longitude, GPS latitude, GPS speed and GPS direction;

S2: limit and identify urban border: extract westernmost end, border, municipal administration region, most the east, northernmost with node coordinate southernmost, node coordinate comprises longitude and latitude; The longitude coordinate arranging westernmost end is the longitude minimum value LON of urban border _min, the longitude coordinate that arranges most the east is the longitude maximal value LON of urban border _max, the latitude maximal value LAT that latitude coordinate is northernmost urban border is set _max, the latitude coordinate arranged southernmost is the latitude minimum value LAT of urban border _min, the region that four boundary values enclose formation is observation scope; Reject the data line not meeting formula (1) in the GPS data from taxi storehouse obtained in step S1;

LON _min<GPS longitude <LON _max, and LAT _min<GPS latitude <LAT _max(1)

S3: data cleansing: by the GPS data from taxi storehouse after step S2 process successively according to gps time, the arrangement of car number Two Variables ascending order, make the data line of same vehicle in the new database formed connect together and in chronological sequence arrange; Delete non-rice habitats data line and malposition data line, obtain effective GPS data from taxi storehouse; Setting global iterative variable i teration is 1; Wherein, the data line that non-rice habitats data line is cancellation mark box-like (2) is deleted:

Rule of judgment 1: operation state=2 or 3, and GPS speed=0 (2) delete position abnormal data capable be delete D _{x, x+1}meet (x+1)th row data time formula (3):

Rule of judgment 2:D _{x, x+1}> V _max(T _x+1-T _x) (3)

D _{x, x+1}shown in (4):

$D_{x,x+1}=2r\&CenterDot;\arcsin \left(\sqrt{{\sin }^2\left(\frac{{\mathrm{\&phi;}}_{x+1}-{\mathrm{\&phi;}}_x}{2}\right)+\cos \left({\mathrm{\&phi;}}_x\right)\cos \left({\mathrm{\&phi;}}_{x+1}\right){\sin }^2\left(\frac{{\mathrm{\&lambda;}}_{x+1}-{\mathrm{\&lambda;}}_x}{2}\right)}\right)---\left(4\right)$

Wherein, D _{x, x+1}for the distance between adjacent lines two data points in the database that obtains after deleting non-rice habitats data line, V _maxfor the maximum speed limit of urban road, r is the radius of the earth, T _x, φ _x, λ _xfor xth row gps time, GPS latitude and GPS longitude in the database that obtains after deleting non-rice habitats data line, T _x+1, φ _x+1, λ _x+1for (x+1)th row gps time, GPS latitude and GPS longitude in the database that obtains after deleting non-rice habitats data line;

S4: the division of observation scope: the observation scope determined by step S2 is cut into the rectangle plot of the capable n row of m, and wherein m and n is obtained by formula (5) and formula (6), the latitude lat of cut-boundary _{limit i}with longitude lon _{limit j}obtained by formula (7) and formula (8) respectively:

m＝INT((LAT _max-LAT _min)/0.0027)(5)

n＝INT((LON _max-LON _min)/0.003528)(6)

Lat _{limit i}=LAT _min+ (i-0.5*iteration) * (LAT _max-LAT _min)/m, (i=2,3,4 ..., m) (7)

Lon _{limit j}=LON _min+ (j-0.5*iteration) * (LON _max-LON _min)/n, (j=2,3,4 ..., n) (8)

S5: the differentiation of crossing: crossing identification is carried out to the rectangle plot obtained in step S4, also all latitude and longitude coordinates composition subdata bases dropping on the effective GPS data from taxi storehouse in a plot are about to, to all subdata bases, according to the span of longitude and latitude, evenly be divided into 10 × 10 detecting units, then judge according to crossing Rule of judgment the position that whether there is crossing and crossing in each detecting unit, the position of the crossing of all identification is recorded in the detection list of crossing, finally assignment and differentiation are carried out to global iterative variable i teration, wherein, crossing Rule of judgment is: the data line comprised in detecting unit is more than 40 row, and the number of unduplicated " GPS orientation " value that these data lines comprise is greater than 40,

S6: the section based on crossing identifies: the crossing detection list obtained by step S5 is converted into intersection stomion his-and-hers watches to be detected, each intersection stomion to be detected forms by initial crossing and termination crossing; Intersection stomion pair to be detected is organized to each, judges initial crossing and stop whether there is road between crossing, if exist, then by initial crossing with stop the coordinate record of crossing in the detection list of section, and delete and repeat to judge section by accident;

S7: export road attribute list: the crossing detection list that output step S5 obtains and the section detection list that step S6 obtains.

2. the road network extraction method based on taxi GPS according to claim 1, it is characterized in that: the position of the crossing in described step S5 comprises crossing longitude and crossing latitude, crossing longitude is the mean value of the longitude of all data lines in detecting unit, and crossing latitude is the mean value of the latitude of all data lines in detecting unit.

3. the road network extraction method based on taxi GPS according to claim 1, it is characterized in that: carrying out assignment and differentiation to global iterative variable i teration in described step S5 is add 1 to global iterative variable i teration and assignment again, if judge that global iterative variable i teration is not equal to 3, then get back to step S4, if global iterative variable i teration equals 3, then enter step S6.

4., according to the road network extraction method based on taxi GPS according to claim 1, it is characterized in that: described step S6 comprises following sub-step:

S6.1: crossing detection list is converted into intersection stomion his-and-hers watches to be detected: according to the distance between the arbitrary crossing in the detection list of formula (4) calculating crossing and all the other crossings, filter out this distance and be less than the intersection stomion of 2km to composition intersection stomion his-and-hers watches to be detected, these his-and-hers watches comprise 5 variablees, i.e. initial crossing longitude lon _start, initial crossing latitude lat _start, stop crossing longitude lon _end, stop crossing latitude lat _end, two crossing distance D _s-e; Wherein, calculate according to formula (9), " s " represents initial crossing, and " e " representative stops crossing:

$D_{s-e}=2r\&CenterDot;\arcsin \left(\sqrt{{\sin }^2\left(\frac{{\mathrm{\&phi;}}_e-{\mathrm{\&phi;}}_s}{2}\right)+\cos \left({\mathrm{\&phi;}}_s\right)\cos \left({\mathrm{\&phi;}}_e\right){\sin }^2\left(\frac{{\mathrm{\&lambda;}}_e-{\mathrm{\&lambda;}}_s}{2}\right)}\right)---\left(9\right)$

S6.2: intersection stomion pair to be detected is organized to each, judge initial crossing and stop whether there is road between crossing: first, limit Road Detection region, this region by initial crossing with stop crossing coordinate line and respectively to offset two straight lines of 20m to both sides and enclose form respectively through initial crossing, the vertical line that stops these two straight lines of crossing; Then surveyed area is divided into along crossing line orientation average individual rectangular sub blocks, does not contain any effective GPS data from taxi storehouse latitude and longitude coordinates point if be less than in 20 Continuous Rectangular sub-blocks, then shows to there is road between two crossings, be recorded in Road Detection list by 5 right for this intersection stomion variablees;

S6.3: delete and repeat to judge section by accident: according to the angle θ stopping crossing and initial crossing line in the Road Detection list that formula (10) calculation procedure S6.2 obtains and formed, wherein quadrant is included in the arctan function considered by atan2; By Road Detection list successively according to angle θ and initial crossing longitude lon _startsequence, makes the data line of identical initial crossing connect together and sorts according to angle value size; Contrast the angle θ that the adjacent data of identical initial crossing is capable, if difference is less than 40 °, then delete larger distance D _s-edata line; Then the value of the angle θ in all data lines between 270 ° and 360 ° is replaced with θ-360 °, again by Road Detection list successively according to angle θ and initial crossing longitude lon _startsequence, contrasts the angle θ that the adjacent data of identical initial crossing is capable, if difference is less than 40 °, then deletes larger distance D _s-edata line;

θ＝atan2(lat _end-lat _start,lon _end-lon _start)(10)。