IDL之GPS坐标转高德地图坐标

15195775117

一、需求缘起


气象上制图，经常是大开大合，地图上有个省市轮廓就行了，分辨率比较粗糙；
但是大气环境监测中，很多项目只在一个城市、一个区县中开展，坐标分辨率大约得精确到10m，
这种对地图的高要求，得使用国内的商用地图，其中据说高德最好。


二、代码与要点


;GPS转高德,可批量，lon经度，lat纬度，可以是向量。（现在经度我已经改用lng）
function GPS_amap,lon,lat


key='XXXXXXXXXXXXXXXXXXX';高德地图web开发k码，可免费申请

;坐标转换的url：
url='https://restapi.amap.com/v3/assistant/'+$
  'coordinate/convert?locations='
for i=0,n_elements(lon)-1 do begin
  url=url+strtrim(string(lon),2)+','+strtrim(string(lat),2)+';'
endfor
  url=url+'&coordsys=gps&output=json&key='+key


;将网页数据按json格式下载到本地
  netObject = IDLnetURL()
  jsonfile='.\geocode.json'
  !null = netObject.Get(URL=url,FILENAME=jsonfile)
  netObject = 0


;按字符串读取该json文件
  openr,lun,jsonfile,/get_lun
  xx=''
  readf,lun,xx
  free_lun,lun


;json解析，字符串转结构体(/TOSTRUCT表示转为结构体)
  json=JSON_PARSE(xx,/TOSTRUCT)
  cut=strsplit(json.locations,';',/extract)
  if n_elements(cut) ne n_elements(lon) then begin
    print,'AMAP坐标数量与GPS不一致!'
    stop
  endif
  
;不转为双精度浮点，则后面追加元素时，小数点位数不足，即精度不够
;如果只是float，则有效数字6位，经纬度精度不够，例如120.115,35.2462
  amaplon=double(1)
  amaplat=double(1)


  for i=0,n_elements(lon)-1 do begin
    cut2=strsplit(cut,',',/extract)
    amaplon=[amaplon,double(cut2[0])]
    amaplat=[amaplat,double(cut2[1])]
  endfor


  amaplon=amaplon[1:-1];高德经度
  amaplat=amaplat[1:-1];高德纬度
  
  re={amaplon:amaplon,amaplat:amaplat}
  return,re


end


调用示例：
IDL> lon=116.481499+dindgen(3)/1000
IDL> lat=39.990475+dindgen(3)/1000
IDL> re=GPS_amap(lon,lat)
IDL> re
{
    "AMAPLON": [116.48758626302100,
116.48858289930600, 116.48957926432300],
    "AMAPLAT": [39.991754014756999,
39.992751193577000, 39.993748101129000]
}