本文详细介绍了如何利用百度地图API获取城市边界数据，并通过C++程序将数据从百度坐标系转换为WGS84坐标系，再用ArcGIS进行处理，最终得到城市边界线的shp文件。过程中涉及数据格式转换、坐标系转换及ArcMap的使用。
参考文章：百度地图API+ArcGIS软件—城市出行时空数据可视化_WenWu_Both的博客-CSDN博客
这篇博客在介绍的时候遗漏了很多关键步骤，我对此进行了必要的补充。
一、获取轮廓线的代码（getData.html）
<!DOCTYPE html
PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>获取地区轮廓线</title>
<script type="text/javascript" src="http://api.map.baidu.com/api?v=2.0&ak=lpAE5ZQGsVIyTHdChEo9CwPnCXyTbCFl">
</script>
<style type="text/css">
body,
html {
width: 100%;
height: 100%;
margin: 0;
font-family: '微软雅黑';
}
#container {
height: 500px;
width: 90%;
margin-left: 50px
}
#Div1 {
width: 100%
}
</style>
</head>
<body>
<div id="container"></div>
<br />
输入省、直辖市或县名称：<input type="text" id="districtName" style="width:80px" value="杭州市">
<input type="button" onclick="getBoundary()" value="获取轮廓线">
<textarea id="Div1" style="width:100%;height:200px"></textarea>
<script type="text/javascript">
var map = new BMap.Map('container');
map.centerAndZoom(new BMap.Point(116.403765, 39.914850), 5);
map.addControl(new BMap.MapTypeControl()); //添加地图类型控件
map.setCurrentCity("杭州市"); // 设置地图显示的城市 此项是必须设置的
map.enableScrollWheelZoom(true); //开启鼠标滚轮缩放
function getBoundary() {
var bdary = new BMap.Boundary();
var name = document.getElementById("districtName").value;
bdary.get(name, function (rs) { //获取行政区域
map.clearOverlays(); //清除地图覆盖物
document.getElementById('Div1').innerText = rs.boundaries;
var count = rs.boundaries.length; //行政区域的点有多少个
for (var i = 0; i < count; i++) {
var ply = new BMap.Polygon(rs.boundaries[i], {
strokeWeight: 2,
strokeColor: "#ff0000"
}); //建立多边形覆盖物
map.addOverlay(ply); //添加覆盖物
map.setViewport(ply.getPath()); //调整视野
}
var blob = new Blob(rs.boundaries, {
type: "text/csv,charset=UTF-8"
})
var csvUrl = URL.createObjectURL(blob)
var aEle = document.createElement("a")
aEle.download = "data.csv" //文件名随意
aEle.href = csvUrl
aEle.click()
});
}
</script>
</body>
</html>注意一下这行代码
<script type="text/javascript" src="http://api.map.baidu.com/api?v=2.0&ak=lpAE5ZQGsVIyTHdChEo9CwPnCXyTbCFl">其中ak=lpAE5ZQGsVIyTHdChEo9CwPnCXyTbCFl
ak=后面的是我申请的百度地图的开发者许可密钥（不知道什么时候会过期，如果过期了建议自己去百度地图申请一个）
打开getData.html，这里以杭州为例（如果想获取别的城市的shp文件，可以自行在文本框内更改），点击获取轮廓线后浏览器会弹出下载文件的弹窗，选择将文件另存为到指定的目录，这里我选择将文件保存到C:/Users/Administrator/Desktop/GIS/data文件夹中。
二、数据转换（dataTransform.cpp）
使用excel打开data.csv，格式如下
为了便于arcgis(arcmap)处理，我们希望得到这种格式的数据
于是我编写了一个C语言程序（dataTransform.cpp）负责格式转换，这里的代码需要根据个人情况修改原文件路径和转换后的文件路径，直接运行会出错。
同时，由于百度地图API获取得到的边界坐标是基于百度坐标系BD09，为了方便后续使用，需要将其转换到通用的WGS84坐标系，有关地理坐标系的相关知识可以参考我的这篇博客
主流地理坐标系、投影坐标系和投影方法的区别和联系_ 一只博客-CSDN博客
https://blog.csdn.net/qq_42276781/article/details/122597363有关BD09和WGS84转换的知识可以阅读我的这篇博客
地理坐标系：WGS84和BD09互转_ 一只博客-CSDN博客https://blog.csdn.net/qq_42276781/article/details/122596796dataTransform.cpp
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define PI 3.1415926535897932384626
double* BD09ToWGS84(double, double);
double* WGS84ToBD09(double, double);
double* BD09ToGCJ02(double, double);
double* GCJ02ToBD09(double, double);
double* GCJ02ToWGS84(double, double);
double* WGS84ToGCJ02(double, double);
int main(){
//input file path
FILE *fp = fopen("C://Users//Administrator//Desktop//GIS//data//data.csv", "r");
//output file path
FILE *fq = fopen("C://Users//Administrator//Desktop//GIS//data//data2.csv", "w");
double lon, lat;
int i = 1;
fprintf(fq, "id, lon, lat\n");
while(fscanf(fp, "%lf, %lf;", &lon, &lat)>0){
double *temp = BD09ToWGS84(lon, lat);
fprintf(fq, "%d, %f, %f\n", i++, temp[0], temp[1]);
}
fclose(fp);
fclose(fq);
}
// (BD09 <=> WGS84) base on (BD09 <=> GCJ02) and (WGS84 <=> GCJ02)
double* BD09ToWGS84(double lon_BD, double lat_BD){
double *GCJ = BD09ToGCJ02(lon_BD, lat_BD);
return GCJ02ToWGS84(GCJ[0], GCJ[1]);
}
double* WGS84ToBD09(double lon_WGS, double lat_WGS){
double *GCJ = WGS84ToGCJ02(lon_WGS, lat_WGS);
return GCJ02ToBD09(GCJ[0], GCJ[1]);
}
// BD09 <=> GCJ02
double baiduFactor = (PI * 3000.0) / 180.0;
double* BD09ToGCJ02(double lon_BD, double lat_BD){
double *GCJ = (double *)malloc(2*sizeof(double));
double x = lon_BD - 0.0065;
double y = lat_BD - 0.006;
double z = sqrt(x * x + y * y) - 0.00002 * sin(y * baiduFactor);
double theta = atan2(y, x) - 0.000003 * cos(x * baiduFactor);
GCJ[0] = z * cos(theta);
GCJ[1] = z * sin(theta);
return GCJ;
}
double* GCJ02ToBD09(double lon_GCJ, double lat_GCJ){
double *BD = (double *)malloc(2*sizeof(double));
double z = sqrt(lon_GCJ * lon_GCJ + lat_GCJ * lat_GCJ) + 0.00002 * sin(lat_GCJ * baiduFactor);
double theta = atan2(lat_GCJ, lon_GCJ) + 0.000003 * cos(lon_GCJ * baiduFactor);
BD[0] = z * cos(theta) + 0.0065;
BD[1] = z * sin(theta) + 0.006;
return BD;
}
// WGS84 <=> GCJ02
double a = 6378245;
double ee = 0.006693421622965823;
// roughly check whether coordinates are in China.
int isInChinaBbox(double lon, double lat){
return lon >= 72.004 && lon <= 137.8347 && lat >= 0.8293 && lat <= 55.8271;
}
double transformLat(double x, double y){
double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(fabs(x));
ret += ((20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0) / 3.0;
ret += ((20.0 * sin(y * PI) + 40.0 * sin((y / 3.0) * PI)) * 2.0) / 3.0;
ret += ((160.0 * sin((y / 12.0) * PI) + 320.0 * sin((y * PI) / 30.0)) * 2.0) / 3.0;
return ret;
}
double transformLon(double x, double y){
double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(fabs(x));
ret += ((20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0) / 3.0;
ret += ((20.0 * sin(x * PI) + 40.0 * sin((x / 3.0) * PI)) * 2.0) / 3.0;
ret += ((150.0 * sin((x / 12.0) * PI) + 300.0 * sin((x / 30.0) * PI)) * 2.0) / 3.0;
return ret;
}
double* delta(double lon, double lat){
double *d = (double *)malloc(2*sizeof(double));
double dLon = transformLon(lon - 105.0, lat - 35.0);
double dLat = transformLat(lon - 105.0, lat - 35.0);
double radLat = (lat / 180.0) * PI;
double magic = sin(radLat);
magic = 1.0 - ee * magic * magic;
double sqrtMagic = sqrt(magic);
d[0] = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * PI);
d[1] = (dLat * 180.0) / ((a * (1.0 - ee)) / (magic * sqrtMagic) * PI);
return d;
}
double* WGS84ToGCJ02(double lon_WGS, double lat_WGS){
double *GCJ = (double *)malloc(2*sizeof(double));
if (!isInChinaBbox(lon_WGS, lat_WGS)){
GCJ[0] = lon_WGS;
GCJ[1] = lat_WGS;
return GCJ;
}
double *d = delta(lon_WGS, lat_WGS);
GCJ[0] = lon_WGS + d[0];
GCJ[1] = lat_WGS + d[1];
return GCJ;
}
double* GCJ02ToWGS84(double lon_GCJ, double lat_GCJ) {
double *WGS = (double *)malloc(2*sizeof(double));
if (!isInChinaBbox(lon_GCJ, lat_GCJ)){
WGS[0] = lon_GCJ;
WGS[1] = lat_GCJ;
return WGS;
}
WGS[0] = lon_GCJ; WGS[1] = lat_GCJ;
double *temp = WGS84ToGCJ02(WGS[0], WGS[1]);
double dx = temp[0] - lon_GCJ;
double dy = temp[1] - lat_GCJ;
while (fabs(dx) > 1e-6 || fabs(dy) > 1e-6) {
WGS[0] -= dx;
WGS[1] -= dy;
temp = WGS84ToGCJ02(WGS[0], WGS[1]);
dx = temp[0] - lon_GCJ;
dy = temp[1] - lat_GCJ;
}
return WGS;
}使用devc++运行，转换后的文件被保存到C://Users//Administrator//Desktop//GIS//data//data2.csv，即
C:\Users\Administrator\Desktop\GIS\data\data2.csv
最终在data文件夹下，有一个data.csv原数据文件，一个data2.csv转换后的数据文件
三、Arcmap处理
如果ArcMap的界面被搞乱了不知道如何还原，可以参考我这篇博客
ArcMap10.7界面初始化_ 一只博客-CSDN博客
https://blog.csdn.net/qq_42276781/article/details/122468720
打开ArcMap后使用默认空白地图，点击确定
点击添加数据，点击连接到文件夹，选择数据所在的文件夹，选择data2.csv，点击添加
右键data2.csv，点击显示XY数据，X字段选择lon，Y字段选择lat，输入坐标的坐标系为地理坐标系=>World=>WGS 1984（由于之前我们已经把BD09转换为WGS84，所以这里需要指定WGS作为地理坐标系），点击确定
点击ArcToolbox，点击数据管理工具，点击要素，点击要素转点，输入要素选择data2.csv，输出要素类指定为data目录下的point.shp文件，点击确定
至此，我们得到了杭州市边界的点shp文件，考虑到边界一般都是线shp文件，所以接下来还需要将点shp文件转为线shp文件。
点击点集转线，输入要素为point.shp，输出要素类为指定目录data下的line.shp，线字段选择ORIG_ID，排序字段选择FID，勾选闭合线，点击确定。
通过勾选这三个图层可以发现，point.shp和data2.csv都是点，而line.shp是线
打开data文件夹，可以发现杭州市的线shp文件已经导出了。
扫一扫
3390
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
