IDL从GPS计算车速

15195775117

附件GPS数据内容为日期-时间-经度-纬度，是GPS装在车上测的，可以计算车速一开始求车速的时候，是把距离和时间一一对应，发现车速震荡很大，时速经常超过200，显然不合理
我试着隔一段求个车速，发现随着间隔的增大，车速曲线的质量先升高后降低，大概间隔8-15位，曲线就比较好看了
下图是在不同位数（Y轴左的数字）滑动平均下车速的变化曲线(km/h)

                               
登录/注册后可看大图

以下几楼我会给出代码和分析

15195775117

Part1.读取GPS数据这里有个精度的问题，IDL计算式经常有8位有效数字的限制，而GPS是精确到小数点后第6位的，大概是米级别，我们发现经度数据共9位有效数字，最后一位会被略去，而纬度数据正好是8位有效数字，不受影响，我把经度数据拆成整数部分和小数部分返回，其实我想着，只用经度的小数部分计算距离没有问题，等于沿着经度线转，这对计算距离没有影响

function GPS,no
;GPS数据的路径：
  gps_file='C:\Users\Administrator\Desktop\car_speed\GPS.txt'
  ;把GPS数据重新整理后存放的文件路径，后面用：
  gps_file2='C:\Users\Administrator\Desktop\car_speed\gps2.txt'
  
  a=read_ascii(gps_file,DELIMITER=';',data_start=0)
  x=a.field1
  sizex=size(x)
  ;如果数据不是二维的就报错，这样情况有可能是数据只有一行或者空文件
  if(sizex[0] ne 2)then begin
    print,'error!'
  endif
  ;行数是raw
  raw=sizex[2]
  openr,lun,gps_file,/get_lun
  x=strarr(1,raw)
  readf,lun,x
  ;把间隔符都统一为分号，写出到GPS2，再重新读取
  x=x.replace('/',';')
  x=x.replace(' ',';')
  x=x.replace(':',';')
  x=x.replace('.',';')
  openw,lun,gps_file2,/get_lun
  printf,lun,x
  close,/all
  a=read_ascii(gps_file2,DELIMITER=';',data_start=0)
  x=a.field01
  x=long(x)
  year=x[0,*]
  month=x[1,*]
  day=x[2,*]
  hour=x[3,*]
  minute=x[4,*]
  sec=x[5,*]
  ;经度的整数部位lon1：
  lon1=double(x[6,*])
  ;经度小数部分lon2：
  lon2=(x[7,*]+0.0)/(10.0^ceil(alog10(x[7,*])))
  ;纬度整数部分lat1：
  lat1=double(x[8,*])
  ;纬度小数部分lat2：
  lat2=(x[9,*]+0.0)/(10.0^ceil(alog10(x[9,*])))
  ;因为纬度整数部分2位+小数部分6位，正好是8位，没有超过8位有效数字，精度可以保障
  ;而经度整数位有3位，多了一位，精度丢失，所以不能直接加
  lat=lat1+lat2
  time=julday(month,day,year,hour,minute,sec)
  ;按结构体输出返回值
  re={year:year,month:month,day:day,hour:hour,minute:minute,sec:sec,$
    lat:lat,time:time,lon1:lon1,lon2:lon2}
  return,re
end

15195775117

Part2.计算车速
;获取经纬度和时间
x=gps(0)
;时间：
t=x.time
;纬度：
lat=x.lat
;经纬度的距离大小不随经度变化，鉴于8位有效数字的限制，用经度小数点后数字代表经度不会影响算出的距离
lon=x.lon2
;以下部分都在该for循环内，查看不同位数滑动平均的效果，每nn个坐标点计算下速度
for nn=1,25 do begin
n=n_elements(lat)
;book用来存放，第一列时间间隔，第二列距离，第三列车速=距离/时间间隔
book=[-99,-99,-99]
for i=0,n-nn-1 do begin
;距离：
d=map_2points(lon,lat,lon[i+nn],lat[i+nn],/meters)
;时间间隔：
time=t[i+nn]-t
;把时间分解成年月日时分秒
caldat,time,month,day,year,hour,minute,sec
;根据对GPS数据的观察，其只精确到秒，所以秒取整
delta_time=round(sec);精确到秒即可，与原始数据对照十分准！
book=[[book],[delta_time,d,d/delta_time]]
endfor
book=book[*,1:-1]
car_spd=transpose(book[2,*])

;画车速曲线：
fig=plot(car_spd*3.6,color='r',yrange=[0,100],dimensions=[4000,400],$
/buffer,xrange=[200,2100],position=[0.05,0.05,0.95,0.95])
;注明滑动平均位数：
t2 = TEXT(0.01,0.5,strtrim(string(nn),2),font_size=50)

endfor

15195775117

Part3.车速水平分布西北和西南路段连续红色部分车速较快，其他部分都是走走停停


代码：
;该部分没变：
x=gps(0)
t=x.time
lat=x.lat
lon=x.lon2
;间隔14个点算个速度
nn=14
n=n_elements(lat)
book=[-99,-99,-99]
for i=0,n-nn-1 do begin
d=map_2points(lon,lat,lon[i+nn],lat[i+nn],/meters)
time=t[i+nn]-t
caldat,time,month,day,year,hour,minute,sec
delta_time=round(sec)
book=[[book],[delta_time,d,d/delta_time]]
endfor
book=book[*,1:-1]
car_spd=transpose(book[2,*])
;以下为新加代码
;匹配颜色，不过自己做的颜色序列IDL不会自动生成色条，可能还要自己做个色条
;不过我是用33号色条，可以直接拷过来，再自己加数字
vert_colors=[-99,-99,-99]
minv=min(car_spd)
maxv=max(car_spd)
for j=0,n-nn-2 do begin
  color=rgb_table_33(car_spd[j],minV,maxv)
  vert_colors=[[vert_colors],[color]]
endfor
vert_colors=vert_colors[*,1:-1]
fig=plot3d(lon[0:-2],lat[0:-2],car_spd*3.6,thick=5,vert_colors=vert_colors)
fig.Rotate, /RESET
fig.Rotate, 0, /ZAXIS
fig.Rotate, 0, /XAXIS
fig.Rotate, 0, /YAXIS
endfor;出图结束

;色条辅助程序，该程序之前没有处理x超过上下限的情况，现在按两端颜色匹配
function rgb_table_33,x,minV,maxV
  MINV=MINV+0.0
  MAXV=MAXV+0.0
  L=(MAXV-MINV)/8.0
  ;如果x超过上下限的情况则按最边沿颜色配色
  if(x lt minv)then begin
    return,[0,0,255]
  endif
  if(x gt maxv)then begin
    return,[255,0,0]
  endif
  IF(X GE MINV AND X LT MINV+L)THEN BEGIN
    R=0.0
    G=0.0
    B=127.5/L*X+127.5
  ENDIF
  IF(X GE MINV+L AND X LT MINV+3*L)THEN BEGIN
    R=0.0
    G=127.5/L*(X-L)
    B=255.0
  ENDIF
  IF(X GE MINV+3*L AND X LT MINV+5*L)THEN BEGIN
    R=127.5/L*(X-3*L)
    G=255.0
    B=-127.5/L*(X-5*L)
  ENDIF
  IF(X GE MINV+5*L AND X LT MINV+7*L)THEN BEGIN
    R=255.0
    G=-127.5/L*(X-7*L)
    B=0.0
  ENDIF
  IF(X GE MINV+7*L AND X LE MAXV)THEN BEGIN
    R=-127.5/L*(X-9*L)
    G=0.0
    B=0.0
  ENDIF
  RETURN,[R,G,B]
end