|
1.背景——为什么要转换用的更多的场景是把经纬度转化为xy平面坐标,因为经纬度是方便我们确定地理位置的,我们可以很容易的从地图数据(可利用高德开放平台)上获取某一个地址它的经纬度,但是我们看到的地图是平面的,所以要利用各种投影把经纬度转换为平面坐标便于我们自己分析~xy平面坐标转化为经纬度的使用场景不多,但我们可以作为一个验证来运用。关于如何从一串地址,如“重庆市沙坪坝区沙正街174号”这样的一个文本地址得出它的经纬度,我在之前发在微信的文章里有过介绍,可见:https://mp.weixin.qq.com/s/Z2I2ufb-AB_gQcfTBD6vmg也可点此跳转2.我的操作环境python3.8我习惯在jupyter里运行python,大家的pycharm、vscode都可3.经纬度和xy平面坐标相互转换----实战importmathimportpandasaspd#定义经纬度转换为米勒坐标的方法defmillerToXY(lon,lat):xy_coordinate=[]L=6381372*math.pi*2#地球周长W=L#平面展开,将周长视为X轴H=L/2#Y轴约等于周长一半mill=2.3#米勒投影中的一个常数,范围大约在正负2.3之间#循环,因为要批量转换forx,yinzip(lon,lat):x=x*math.pi/180#将经度从度数转换为弧度y=y*math.pi/180#将纬度从度数转换为弧度y=1.25*math.log(math.tan(0.25*math.pi+0.4*y))##这里是米勒投影的转换x=(W/2)+(W/(2*math.pi))*x#这里将弧度转为实际距离,转换结果的单位是kmy=(H/2)-(H/(2*mill))*y#这里将弧度转为实际距离,转换结果的单位是kmxy_coordinate.append((int(round(x)),int(round(y))))returnxy_coordinate#xy坐标转换成经纬度的方法(该方法未定义循环,仅能单个坐标转换defxy_to_coor(x,y):lonlat_coordinate=[]L=6381372*math.pi*2W=LH=L/2mill=2.3lat=((H/2-y)*2*mill)/(1.25*H)lat=((math.atan(math.exp(lat))-0.25*math.pi)*180)/(0.4*math.pi)lon=(x-W/2)*360/W#TODO最终需要确认经纬度保留小数点后几位lonlat_coordinate.append((round(lon,8),round(lat,8)))returnlonlat_coordinate#读取数据文件df=pd.read_csv('village.DAT',header=None,encoding='utf-8',delimiter='')#delimiter=''代表分隔符是空格df.to_csv('village.DAT',header=None,encoding='utf-8',sep='')#取一下经纬度的数据x_data=df.iloc[:,2]#对应第3列数据,即经度数据y_data=df.iloc[:,3]#对应第4列数据,即纬度数据#调用经纬度转化为xy坐标方法,直接输出转换后的xy坐标xy_data=millerToXY(x_data,y_data)xy_data#直接输出#调用xy坐标转化为经纬度方法,输出经纬度print(xy_to_coor(31327850,8251850))#这里直接随机给了一个xy坐标参数1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950问题:把经纬度坐标转换为xy平面坐标后,我们会发现新坐标的数值变得非常大,比如(31327850,8251850),这非常不利于我们的可视化过程,比如把各点在一个坐标系上标注出来,那怎么解决呢?常见的解决方案:我们可以根据原来的经纬度,找一个和所有点靠的比较近,能在一张图上的点,以这个点作为原点,重新建立一个平面的直角坐标系,这样所有点的坐标都会比较好看,不会那么大,也不至于离原点太远。同时需注意这个新坐标原点的选择,可以尽量选在所有点的左下方,这样可以确保咱的坐标都是正数,方便后面分析~4.在xy平面坐标上,以一个新点为原点,重新建立坐标系,输出各点新坐标importmathimportpandasaspd#方法:以一个新点为原点,重新建立坐标系,输出各点新坐标defxyToNew(lon,lat):xy_coordinate=[]#这里给下坐标原点,我是随机选的,使得所有点的坐标都为正,且尽量大于1x0=31327850y0=8251850forx,yinzip(lon,lat):x=x-x0y=y-y0xy_coordinate.append((int(round(x)),int(round(y))))returnxy_coordinate#读取转换成平面坐标后的坐标文件df=pd.read_csv('xymile.txt',header=None,encoding='utf-8',delimiter='')df.to_csv('xymile.txt',header=None,encoding='utf-8',sep='')#df#可以打印下看看这个df数据有啥#取一下xy平面坐标系上xy的数据x_data=df.iloc[:,0]#对应第1列,即xy_data=df.iloc[:,1]#对应第2列,即y#调用转换坐标的方法xy_data=xyToNew(x_data,y_data)xy_data#直接输出坐标1234567891011121314151617181920212223242526272829
|
|