矢量地图坐标转换方法及工具.doc

矢量地图坐标转换方法及工具在MapInfo中，通常采用大地坐标的北京54标准。与GPS采用的WGS84在同一点上相差经纬度相差有1020分（仅仅从地图在MAPINFO的显示数据来说没涉及到NumericCoordSys的坐标系统）。虽然现有GIS平台中都预定义有上百个基准面供用户选用，但均没有我们国家的基准面定义。假如精度要求不高，可利用前苏联的Pulkovo 1942基准面(Mapinfo中代号为1001)代替北京54坐标系；假如精度要求较高，如土地利用、海域使用、城市基建等GIS系统，则需要自定义基准面。实际工作中一般都根据工作区内已知的北京54坐标控制点计算转换参数，如果工作区内有足够多的已知北京54与WGS84坐标控制点，可直接计算坐标转换的7参数或3参数；当工作区内有3个已知北京54与WGS84坐标控制点时，可用下式计算WGS84到北京54坐标的转换参数(A、B、C、D、E、F)：x54 = AX84 + BY84 + C，y54 = DX84 + EY84 + F，多余一点用作检验；在只有一个已知控制点的情况下(往往如此)，用已知点的北京54与WGS84坐标之差作为平移参数，当工作区范围不大时精度也足够了。一般使用3参数，将3个控制点的数据代入方程，求解出参数，可以利用Excel求解，如图在黄色区域分别输入3个点的不同坐标系的数据，通过消元法，在红色区域会得出a、b、c、d、e、f参数的值，可以在下图的黄色区域输入数据进行验证，检验误差为了实现地图数据的转换，我用MapBasic编写了一个程序。该程序在MapInfo程序界面中添加了一个转换的菜单，主要包括2部分：一是参数设置，用于设置并保存用于转换的3参数；二是执行转换，选择需要进行转换的表并转换。如图：参数设置内容如图：按钮可以将参数恢复到缺省值，点击按钮会将输入的参数设置进系统，并保存为config.ini文件。转换内容如图：可以选择需要转换的表，如果点击全选则自动全部选中所有表，再次点击则取消全部选中。该程序具有记录操作的功能，将主要的操作记录到文件中。附件：源程序* 坐标转换* Created by 孙国文,江苏电信徐州分公司无线中心* Jan 28, 2005include mapbasic.defdeclare sub maindeclare sub RCVdeclare sub CBdeclare sub HBdeclare sub editlayerdeclare sub minmapdeclare sub restoremapdeclare sub checktabdeclare sub objmovedeclare sub deltabdeclare sub YesButtondeclare sub NoButtondeclare sub EndProdeclare sub Aboutdeclare sub provisiondeclare sub objmove_maindeclare sub getconffromfiledeclare sub setdefaultconfigdeclare Sub getresdeclare Sub changestatemultideclare Sub changestatecheckdim r,v,i,j,h,poly,pnode,cnt as integerdim z as Smallintdim numtab,tabtemp,tabnum,numlayer,objoff,numwin,mx,my,arccnt as integerdim ptcnt,plcnt,lcnt,regcnt,filecnt,rctcnt,rrtcnt,elpcnt,txcnt,ewindow as integerdim sobj as objectdim nx,ny,xoff,yoff,go1,go2,go3,go4,tsp,ta,tlx,tly,gx1,gy1,gx2,gy2 as floatdim tabstr(1) as stringdim tabname,runstr,ts,elayer,ct as stringdim tabobj as aliasdim layerfound,cancel as logicaldim filenum(1) as integerdim tf as fontdim tj,tar as smallintdim tjs as string存储参数dim da,db,dc,dd,de,df as floatdim a,b,c,d,e,f as float配置文件dim configfile,logfile as string需要转换的表数dim num2conv as integersub mainconfigfile=config.inilogfile=convert.logopen file logfile for append as #2print curdate()+ +time(24)+ 开始print #2,curdate()+ +time(24)+ 开始x84=a*x54+b*y54+cy84=d*x54+e*y54+e da = 1.00088323 db = 0.002918798 dc = -0.202027599 dd = -0.000160498 de = 1.001210357 df = -0.022304691 Create Menu 坐标转换 As 参数设置 calling provision, 转换 calling objmove_main, Help calling HB, 关于. calling About, 退出 calling EndPro Alter Menu Bar Add 坐标转换end subsub About dialog title 关于 坐标转换 control statictext position 10,10 title 坐标转换 control statictext position 10,25 title Copyright 01/28/2006 - 孙国文 control statictext position 10,35 title All rights reserved. control statictext position 10,50 title 本工具用于实现矢量图坐标系的转换或地图偏差的修正。 control statictext position 10,60 title 使用方法：在参数设置中填入转换参数(缺省参数可用于徐州地区54转84)， control statictext position 10,70 title 打开需要进行转换的表(直接打开，不要设置缩放范围或通过.wor打开) control statictext position 10,80 title 执行转换菜单，而后耐心等待程序执行 control statictext position 10,90 title 注意：本工具目前仅适用于矢量图的转换，不支持栅格图的转换。 control statictext position 10,100 title 本工具直接修改原数据，因此转换前请“务必手工备份”地图数据 control statictext position 10,110 title 作者：孙国文 control statictext position 10,120 title 联系电话E-mail： control statictext position 10,130 title 欢迎对GIS有兴趣的朋友和我交流 control OKButtonend subsub objmove_mainprint chr$(12) call checktab if numtab=0 then note 请先打开需要转换的table,再进行转换操作! print curdate()+ +time(24)+ 请先打开需要转换的table,再进行转换操作! print #2,curdate()+ +time(24)+ 请先打开需要转换的table,再进行转换操作! exit sub end ifDialog Title 转换 Width 250 Height 80 Control StaticText Position 8,7 Title 选择需要转换的表: Control MultiListBoxPosition 20,18Title from Variable tabstrWidth 80Height 100ID 1Calling changestatecheck Value 1 Control CheckBoxPosition 20,120Title 全选Width 30Height 16ID 2Calling changestatemultiValue False Control OKButtonPosition 10,140Width 50Height 16ID 10Title OK Calling getres Control CancelButtonPosition 70,140Width 50Height 16ID 11Title Cancelend subsub objmove call checktab cancel = FALSE lcnt = 0 plcnt = 0 regcnt = 0 ptcnt = 0 filecnt = filecnt + 1print before ewindow +ewindow set map redraw off call editlayerprint after ewindow +ewindow call minmap Set ProgressBars Off Randomize If filecnt 1 then redim filenum(filecnt) End If filenum(filecnt) = int(100000 * rnd(1) commit table tabname as ProgramDirectory$()+MM+filenum(filecnt) Note Table +tabname+ saved as +ProgramDirectory$()+MM+filenum(filecnt) cnt = 1 runstr = select * from +tabname+ where obj into gotobj run command runstr fetch first from gotobj do while not eot(gotobj) tabobj = gotobj+.obj sobj = tabobj If str$(sobj) = Region thenregcnt = regcnt + 1 End If If str$(sobj) = Polyline then plcnt = plcnt + 1 End If If str$(sobj) = Rectangle or str$(sobj) = Rounded Rectangle or str$(sobj) = Ellipse thensobj = ConvertToRegion(sobj)If str$(sobj) = Rectangle then rctcnt = rctcnt + 1End IfIf str$(sobj)= Rounded Rectangle then rrtcnt = rrtcnt + 1End IfIf str$(sobj) = Ellipse then elpcnt = elpcnt + 1End If End If If str$(sobj) = Arc thensobj = ConvertToPline(sobj)arccnt = arccnt + 1 End If If str$(sobj) = Region or str$(sobj) = Polyline then If str$(sobj) = Region then objoff = 1 Else objoff = 0End If poly = ObjectInfo(sobj,obj_info_npolygons) for i = 1 to poly pnode = ObjectInfo(sobj,obj_info_npolygons+i) for j = 1 to pnode-objoff nx = ObjectNodeX(sobj,i,j) ny = ObjectNodeY(sobj,i,j) Alter Object sobj Node Set Position i,j (a*nx+b*ny+c,d*nx+e*ny+f) Update gotobj set obj = sobj where RowID = cnt next next End If If str$(sobj) = Line then Alter Object sobj Geography obj_geo_linebegx,a*ObjectGeography(sobj,obj_geo_linebegx)+b*ObjectGeography(sobj,obj_geo_linebegy)+c Alter Object sobj Geography obj_geo_linebegy,d*ObjectGeography(sobj,obj_geo_linebegx)+e*ObjectGeography(sobj,obj_geo_linebegy)+f Alter Object sobj Geography obj_geo_lineendx,a*ObjectGeography(sobj,obj_geo_lineendx)+b*ObjectGeography(sobj,obj_geo_lineendy)+c Alter Object sobj Geography obj_geo_lineendy,d*ObjectGeography(sobj,obj_geo_lineendx)+e*ObjectGeography(sobj,obj_geo_lineendy)+f Update gotobj set obj = sobj where RowID = cnt lcnt = lcnt + 1 End If If str$(sobj) = Point then Alter Object sobj Geography obj_geo_pointx,a*ObjectGeography(sobj,obj_geo_pointx)+b*ObjectGeography(sobj,obj_geo_pointy)+c Alter Object sobj Geography obj_geo_pointy,d*ObjectGeography(sobj,obj_geo_pointx)+e*ObjectGeography(sobj,obj_geo_pointy)+f Update gotobj set obj = sobj where RowID = cnt ptcnt = ptcnt + 1 End if If str$(sobj) = Text thenset map window ewindow layer elayer editable ongo1 = ObjectGeography(sobj,obj_geo_minx)go2 = ObjectGeography(sobj,obj_geo_miny)go3 = ObjectGeography(sobj,obj_geo_maxx)go4 = ObjectGeography(sobj,obj_geo_maxy)ts = ObjectInfo(sobj,obj_info_textstring)tf = ObjectInfo(sobj,obj_info_textfont)print tf is +tftsp = ObjectInfo(sobj,obj_info_textspacing)tj = ObjectInfo(sobj,obj_info_textjustify) Do case tj case 0 tjs=LEFT case 1 tjs=CENTER case 2 tjs=RIGHT End case ta = ObjectGeography(sobj,obj_geo_textangle)gx1 = a*go1+b*go2+cgy1 = d*go1+e*go2+fgx2 = a*go3+b*go4+cgy2 = d*go3+e*go4+fif ObjectInfo(sobj,obj_info_textarrow) 0 then tar = ObjectInfo(sobj,obj_info_arrow) tlx = ObjectGeography(sobj,obj_geo_textlinex) tly = ObjectGeography(sobj,obj_geo_textliney) ct = Create Text Into Window +ewindow+ +ts+ (+gx1+,+gy1+) +(+gx2+,+gy2+) +tf+ label line +tar+ (+tlx+,+tly+) spacing +tsp+ justify +tjs+ angle +taprint ct1 = +ct run command ctelse ct = Create Text Into Window +ewindow+ +ts+ (+gx1+,+gy1+) +(+gx2+,+gy2+) +tf+ spacing +tsp+ justify +tjs+ angle +taprint ct2 = +ct run command ctend iftxcnt = txcnt + 1 End If commit table tabname fetch next from gotobj cnt = cnt + 1 loop call restoremap set map redraw off note regcnt+ region(s) moved+chr$(13)+plcnt+ polyline(s) moved+chr$(13)+lcnt+ line(s) moved+chr$(13)+ptcnt+ point(s) moved print curdate()+ +time(24)+ +regcnt+ region(s) moved+chr$(13)+plcnt+ polyline(s) moved+chr$(13)+lcnt+ line(s) moved+chr$(13)+ptcnt+ point(s) moved print #2,curdate()+ +time(24)+ +regcnt+ region(s) moved+chr$(13)+plcnt+ polyline(s) moved+chr$(13)+lcnt+ line(s) moved+chr$(13)+ptcnt+ point(s) moved drop table gotobjend subsub checktab numtab = NumTables() redim tabstr(numtab) for h = 1 to numtab if TableInfo(h, tab_info_mappable) = TRUE then tabstr(h) = TableInfo(h, tab_info_name) end if nextend subsub editlayer print numwin= +numwin for z = 1 to numwin if layerfound = TRUE thenexit for end if numlayer = MapperInfo(WindowID(z),mapper_info_layers) print numlayer=+numlayer for v = 1 to numlayer if LayerInfo(WindowID(z),v,layer_info_name) = tabname then elayer = LayerInfo(WindowID(z),v,layer_info_name) ewindow = WindowID(z) layerfound = TRUE exit for end if next nextend subsub minmap numwin = NumWindows() layerfound = FALSE for z = 1 to numwin* Boolean value - if the inside for.next loop was abnormally terminated. if layerfound = TRUE thenexit for end if if WindowInfo(z, win_info_type) = 1 thenset window windowid(z) min* The following section takes the chosen layer out of view.* It was removed in favor of minimizing all map windows to reduce* redraw time. numlayer = MapperInfo(WindowID(z),mapper_info_layers)for v = 1 to numlayer if LayerInfo(WindowID(z),v,layer_info_name) = tabname then remove map window windowid(z) layer v layerfound = TRUE exit for end if next end if nextend subsub restoremap for z = 1 to numwin if WindowInfo(z, win_info_type) = 1 thenset window windowid(z) restore end if nextend subsub EndPro print curdate()+ +time(24)+ 退出程序 print #2,curdate()+ +time(24)+ 退出程序 close file #2 end programend subsub provisionprint chr$(12)if fileexists(configfile) then call getconffromfileelse print curdate()+ +time(24)+ 不存在配置文件，使用缺省值 print #2,curdate()+ +time(24)+ 不存在配置文件，使用缺省值 call setdefaultconfigend if Dialog Title 设置参数 Width 250 Height 80 Control StaticText Position 8,7 Title x84=a: Control EditTextPosition 32,7Width 50Height 12ID 1Into aValue a Control StaticText Position 85,7 Title *x54+b: Control EditTextPosition 110,7Width 50Height 12ID 2Into bValue b Control StaticText Position 160,7 Title *y54+c: Control EditTextPosition 187,7Width 50Height 12ID 3Into cValue c Control StaticText Position 8,25 Title y84=d: Control EditTextPosition 32,25Width 50Height 12ID 4Into dValue d Control StaticText Position 85,25 Title *x54+e: Control EditTextPosition 110,25Width 50Height 12ID 5Into eValue e Control StaticText Position 160,25 Title *y54+f: Control EditTextPosition 187,25Width 50Height 12ID 6Into fValue f Control OKButtonPosition 85,57Width 60Height 16ID 10Title OK calling RCV Control CancelButtonPosition 150,57Width 60Height 16Title Cancel Control ButtonPosition 18,57Width 60Height 16Title Default calling setdefaultconfig print x84=+a+*x54+b+*y54+c print y84=+d+*x54+e+*y54+fcall YesButton end subsub YesButtonnote OK Button is pressedprint curdate()+ +time(24)+ x84=+a+*x54+b+*y54+cprint curdate()+ +time(24)+ y84=+d+*x54+e+*y54+fprint #2,curdate()+ +time(24)+ x84=+a+*x54+b+*y54+cprint #2,curdate()+ +time(24)+ y84=+d+*x54+e+*y54+fend subsub getconffromfiledim str as stringopen file configfile for input as #1line input #1,strif eof(1) then print curdate()+ +time(24)+ 配置文件错误,使用缺省值 print #2,curdate()+ +time(24)+ 配置文件错误,使用缺省值 call setdefaultconfigelse input #1,a,b,c,d,e,f print curdate()+ +time(24)+ 使用配制文件中参数设置 print #2,curdate()+ +time(24)+ 使用配制文件中参数设置end ifclose file #1end subsub setdefaultconfiga=dab=dbc=dcd=dde=def=dfprint curdate()+ +time(24)+ 使用缺省值print #2,curdate()+ +time(24)+ 使用缺省值end subsub RCVa = ReadControlValue(1)b = ReadControlValue(2)c = ReadControlValue(3)d = ReadControlValue(4)e = ReadControlValue(5)f = ReadControlValue(6)open file configfile for output as #1write #1,#a,b,c,d,e,fwrite #1,a,b,c,d,e,fclose file #1print curdate()+ +time(24)+ 保存配置文件成功print #2,curdate()+ +time(24)+ 保存配置文件成功call YesButtonend subSub changestatecheckdim i