ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

成渝城市群矢量数据shp坐标系转换与清洗避坑指南

成渝城市群矢量数据shp坐标系转换与清洗避坑指南 简介成渝城市群矢量数据包面向城市规划、地理信息分析与科研项目人员提供学校、医院等兴趣点以及道路、铁路、建筑轮廓等空间要素可作为空间分析与制图的基础底图。包内共一百五十五个文件压缩包约四十点三兆主体为二十二套shp矢量文件每套均配齐shx索引文件、dbf属性表、prj投影信息等配套格式方便直接加载到常用地理信息系统软件。数据年限为二零一九年其中兴趣点数据通过网络地图抓取并裁剪覆盖幼儿园、中小学、高校以及综合医院、专科医院、动物医院等类别道路、铁路与建筑轮廓则来自开放街道地图并按成渝城市群边界完成裁剪整体采用WGS84椭球投影。该数据包免去自行爬取和坐标转换流程适宜开展设施分布、可达性分析、路网规划等研究。已有约一千一百九十一人学习下载适合需要现成矢量底图的科研人员。1. 成渝城市群矢量数据不是一堆shp文件而是一套可落地的空间底座做区域经济分析、商业选址、交通可达性研究时成渝城市群矢量数据几乎是绕不开的底图。很多人以为拿到几个shp文件就能开工结果导入ArcGIS后要么边界跑到贵州要么POI点叠不上道路要么建筑轮廓一缩放就变形。这套数据实际包含行政区边界、道路网络、建筑轮廓、POI四类核心图层分别解决“范围在哪”“怎么过去”“哪里有人活动”“具体设施是什么”的问题。对规划师、GIS开发者和数据分析师来说关键是把它整理成坐标系一致、字段规整、拓扑干净的数据集。下面我按自己处理成渝数据的流程把获取、清洗、转换、避坑和验证讲清楚。2. 成渝城市群shp数据从哪来公开数据源、坐标系差异与选型清单2.1 四大类数据的常见来源与精度取舍成渝城市群覆盖重庆27个区县和四川15个市总面积约18.5万平方公里。不同图层的数据来源和精度差异很大先搞清楚每个图层的“最优解”再动手能省掉后面一半的返工。行政区边界shp最权威的来源是国家基础地理信息中心发布的1:100万公众版数据以及各省民政厅公布的行政区划图。这些数据的边界精度在百米级用于城市群尺度展示和统计汇总足够但如果要做某个街道的用地分析就需要更高精度的县级或乡镇级边界。另一个常用来源是OpenStreetMapOSM导出的国界和省界优点是免费、更新快缺点是边界在局部存在锯齿或未闭合问题。道路数据我一般优先用OSM的公路图层分类字段齐全motorway、trunk、primary等成渝两地的城市快速路和高速基本都有。注意OSM道路是中心线不是路面边界宽度需要按车道数估算。建筑轮廓数据有两个来源一是OSM的building多边形在重庆主城区和成都绕城内覆盖较好郊区缺失严重二是微软和高德公开的建筑足迹数据但成渝地区的开放建筑数据集不如北美和欧洲完整。实际项目中我会把OSM和本地天地图对比着看缺的地方用遥感影像手动补绘。POI数据是这四类里最“活”的也是最需要清洗的。常见来源包括高德开放平台、百度地图开放平台和OSM的amenity图层。高德POI经过GCJ-02加密不能和WGS-84坐标的边界、道路直接叠加OSM的POI是WGS-84坐标但类别少、覆盖稀疏。我的建议是如果分析精度要求不高直接用OSM的POI如果要商业选址得从地图API拉取高德或百度POI并做好坐标纠偏。提示公开数据源的坐标系统一是个大问题。成渝城市群范围内天地图、国家基础地理信息中心产品多为CGCS2000或WGS84高德、百度是加密坐标系OSM是WGS84。第2.2节会讲怎么统一。2.2 拿到手先统一坐标系WGS84还是GCJ02成渝城市群横跨105°E到108°E按高斯-克吕格投影涉及3度带的35、36带按6度带涉及18、19带。如果你只想做可视化直接用经纬度WGS84即可但要做面积量算、缓冲区分析、网络分析必须投影到平面坐标系否则一个缓冲区会被拉成椭圆。常见做法是先统一到WGS84地理坐标系再根据分析范围投影到合适的UTM或高斯投影。对于成渝城市群这种大范围区域我一般用Albers等积投影或兰伯特投影避免面积变形。在Python中可以用geopandas统一坐标。下面是一段从GCJ02转WGS84再转Albers的代码成渝POI基本都会用到。import geopandas as gpd import pandas as pd from pyproj import Transformer import numpy as np # 读取高德导出的POI CSV包含lng/latGCJ02 df pd.read_csv(chengyu_poi_gcj02.csv) # 定义GCJ02 - WGS84 的转换器 trans_gcj_to_wgs Transformer.from_crs(EPSG:4490, EPSG:4326, always_xyTrue) def gcj02_to_wgs84(lng, lat): # 这里用pyproj无法直接处理GCJ02需要先转成EPSG:4490CGCS2000 # 实际做法是先用算法解算偏移然后再投影 return trans_gcj_to_wgs.transform(lng, lat) # 简化的火星坐标系转WGS84单点转换示例 def simple_gcj02_to_wgs84(lng, lat): # 如果手头没有加密规则可以用近似网格查表法 # 实际项目中我通常用coord_convert库处理 return lng - 0.0035, lat - 0.0020 # 仅示意不用于生产 df[wgs_lng], df[wgs_lat] zip(*df.apply(lambda r: simple_gcj02_to_wgs84(r[lng], r[lat]), axis1)) # 转成GeoDataFrame并设置原始坐标为WGS84 gdf gpd.GeoDataFrame(df, geometrygpd.points_from_xy(df[wgs_lng], df[wgs_lat]), crsEPSG:4326) # 投影到适合成渝城市群的Albers等积投影 chengyu_albers gdf.to_crs(projaea lat_125 lat_233 lat_00 lon_0105 x_00 y_00 datumWGS84) chengyu_albers.to_file(chengyu_poi_albers.shp, encodingutf-8)这段代码的逻辑是先读CSV再对经纬度做火星坐标转换最后设置CRS并投影。实际中GCJ02转WGS84不是简单地减去固定值成渝地区偏移量在50到500米不等必须用完整的火星坐标算法或调用公共库。我一般用coord_convert或gcj02_to_wgs84的公开实现避免手写公式出错。参数说明projaea定义Albers等积投影lat_1和lat_2是标准纬线lon_0取成渝中间经度105°E。如果你只处理重庆主城可以改用UTM 48N或高斯3度带35带面积精度会更高。2.3 用QGIS快速查看与属性检查拿到shp后先用QGIS做一次快速体检。我通常在QGIS里按图层加载叠加“四川省重庆市”边界然后查看每个图层的属性表、坐标系和几何类型。步骤很简单图层右键 → Properties → Source看CRS和Encoding用Check geometries插件跑一遍拓扑检查。成渝数据最常见的两个毛病就在这一步暴露边界shp的Geometry collection类型和POI的重复点。QGIS里还可以用“Field Calculator”快速统计每个区县的POI数量验证数据完整性。比如要检查成渝城市群21个地市有没有遗漏用县界shp和POI做join后count一下如果某地市count为0大概率是边界或POI坐标系不匹配。3. 把原始数据做成可分析的城市群数据集字段整理、空间连接与POI清洗3.1 POI数据的清洗与去重坐标系、类别映射POI数据是所有图层里最脏的。从高德或百度拉下来的POI经常有重复记录、字段缺失、类别混乱。比如同一家银行网点既出现在“金融”又出现在“商务住宅”或者“餐饮”和“美食”两个类别并存。第一步是统一类别体系。我用高德POI时会把它的大类字段type归一成自己定义的category_level1和category_level2。比如“餐饮服务;中餐厅;川菜”拆成一级“餐饮”二级“川菜”。这个映射表需要人工维护成渝地区特色类别很多“火锅店”“茶馆”要单独列出来。下面是用pandas清洗POI并去重的代码。import pandas as pd import geopandas as gpd df pd.read_csv(chengyu_poi_raw.csv) # 按名称地址类别去重保留首次出现的记录 df[dup_key] df[name].str.strip() | df[address].str.strip() | df[type] df df.drop_duplicates(subsetdup_key, keepfirst).drop(columnsdup_key) # 类别映射 def map_category(t): t str(t) if 餐饮 in t or 美食 in t: return 餐饮 if 购物 in t or 商场 in t: return 购物 if 医疗 in t or 医院 in t: return 医疗 if 教育 in t: return 教育 return 其他 df[category_level1] df[type].apply(map_category) # 校正POI点落在行政区外的情况用边界shp做空间过滤 boundary gpd.read_file(chengyu_boundary.shp) gdf gpd.GeoDataFrame(df, geometrygpd.points_from_xy(df[wgs_lng], df[wgs_lat]), crsEPSG:4326) gdf gpd.sjoin(gdf, boundary[[region, geometry]], howinner, predicatewithin) gdf gdf.drop_duplicates(subset[name, address, category_level1]) gdf.to_file(chengyu_poi_clean.shp, encodingutf-8)这里的去重逻辑不只是按名称因为成渝有很多“分店”同名不同址。加上地址和类别能降低误杀如果两个POI距离很近但地址写法不同“XX路1号”和“XX路01号”还需要先做地址标准化。空间过滤是处理“坐标偏移导致POI飞出去”的兜底手段配合sjoin的predicatewithin把落在边界外的点剔除。但注意如果边界本身是GCJ02WGS84的点会被剔除一半所以顺序一定是先统一坐标系再做sjoin。3.2 道路与建筑轮廓的拓扑修复道路和建筑轮廓shp拿到后经常有拓扑错误。OSM建筑多边形的典型问题是自相交和不闭合。自相交会导致面积计算错误QGIS里显示正常但做Overlay分析时直接报错。GIS软件里都有修复工具ArcGIS的Repair Geometry和QGIS的Fix Geometries。我更推荐用geopandas批量处理尤其在数据集很大的时候。import geopandas as gpd from shapely.validation import make_valid gdf gpd.read_file(chengyu_building_osm.shp) # 记录修复前的几何类型分布 print(gdf.geometry.type.value_counts()) # 批量修复无效几何 def fix_geom(geom): if geom is None: return None if not geom.is_valid: # make_valid会处理自相交、洞等问题但可能返回几何集合 return make_valid(geom) return geom gdf[geometry] gdf[geometry].apply(fix_geom) # 把多几何拆成单个多边形避免后续叠加分析翻车 gdf gdf.explode(index_partsFalse) gdf gdf[gdf.geometry.type Polygon] gdf gdf[gdf.geometry.is_valid] # 面积检查成渝建筑一般小于5000平米过大说明几何有问题 gdf gdf.to_crs(EPSG:4326) gdf[area_km2] gdf.geometry.area * 111.32 * 111.32 # 经纬度下的近似面积仅筛查用 outlier gdf[gdf[area_km2] 5] print(f异常大建筑 {len(outlier)} 个)make_valid不是万能的它可能把自相交多边形拆成MultiPolygon甚至GeometryCollection所以修复后要explode并按Polygon过滤。经纬度坐标直接area算出来的是平方度不能当真实面积我这里的*111.32*111.32只是粗筛异常值精确面积必须用投影坐标系重算。道路数据则重点检查悬挂节点和伪节点用unary_union后的线段是否断裂来判断。建筑轮廓修复后还要做“负缓冲”检查。成渝和全国一样OSM建筑轮廓常包含公共边两个相邻建筑的边完全重合这本身不是错误但做叠加分析时会被算成重叠面积。如果要做建筑密度统计必须先对要素做dissolve再算每个地块内的建筑覆盖而不是直接统计shp里的多边形面积。3.3 用空间连接给POI挂接区县与城市群边界整理好的POI只有经纬度和类别还不够分析时要能回答“哪个区县”“属于哪个城市群片区”。这就需要把区县边界shp的属性挂到POI上。空间连接用geopandas的sjoin是最顺手的。但要注意连接时选择predicatewithin还是intersects——POI点落在边界线上时intersects会返回多个匹配导致重复。import geopandas as gpd poi gpd.read_file(chengyu_poi_clean.shp) county gpd.read_file(chengyu_county.shp) # 统一坐标系再sjoin poi poi.to_crs(EPSG:4326) county county.to_crs(EPSG:4326) # 点在哪条边界里就挂哪个区县 joined gpd.sjoin(poi, county[[county_name, city_name, geometry]], howleft, predicatewithin) # 用drop_left_index避免重复 joined joined.drop_duplicates(subset[name, address, category_level1], keepfirst) # 检查没挂上区县的POI missing joined[joined[county_name].isna()] print(f未匹配到区县的POI数量: {len(missing)})如果未匹配的点很多先别急着删除看看它们落在哪。成渝边界附近存在飞地比如重庆潼南和四川遂宁交界处有些乡镇飞地会脱离主边界更常见的是POI坐标偏移恰好压线。我此时会把predicate换成intersects再试一次## 4. 从dwg到shp再到3D场景格式转换与坐标系偏移排查4.1 dwg转shp的常见工作流很多成渝本地规划项目的历史数据是CAD的dwg格式尤其是建筑轮廓和道路中线。dwg里的实体没有“要素”概念只有线、块、多段线直接转shp会遇到图层混乱的问题。我的习惯是先在CAD里整理图层再转出。dwg转shp有两条路一是用FME或ArcGIS的“CAD到地理数据库”工具二是用QGIS的AutoCAD DXF/DWG导入插件。前者对复杂块和标注的支持更好后者胜在免费。转之前必须做三件事炸开所有块、把文字标注放在独立图层、统一单位毫米还是米。成渝测绘项目常用“成都坐标系”和“重庆独立坐标系”转shp时要知道它们的EPSG代码或转换参数否则坐标会偏几公里。实际项目里我见过最坑的是dwg里“坐标是米但图形画在毫米单位上”转shp后全图缩小1000倍。解决办法是先用CAD的SCALE命令除以1000再用QGIS加载检查。对于带高程的dwg转shp时Z值会保留但建筑轮廓不需要Z值的话建议在CAD里把Z全部置0免得后续分析出现三维几何。转换后需要用第3.2节的拓扑修复流程处理一遍。dwg转出的多边形几乎必然有重复节点、极小碎边和自相交。我在QGIS中处理成渝dwg数据的标准流程是导入 → 检查坐标系 → 炸开 → 用v.clean清理 → 输出shp。4.2 shp转WKT/TXT导出给非GIS场景成渝城市群数据不光在GIS软件里用很多数据分析师和算法工程师需要把shp转成文本格式比如WKT或GeoJSON才能喂给PostgreSQL或机器学习模型。shp转WKT常用ogr2ogr或geopandas的to_wkt。# 用ogr2ogr把shp转成GeoJSON再转WKT ogr2ogr -f GeoJSON chengyu_boundary.geojson chengyu_boundary.shpogr2ogr的-t_srs参数可以在转换时重投影但要注意源shp的.prj文件是否完整。如果srcshp缺少prj文件一定要在-s_srs里手动指定否则转换结果会是错的。用Python转WKT是另一个方案适合把POI导出给后端服务。import geopandas as gpd gdf gpd.read_file(chengyu_poi_clean.shp) gdf gdf.to_crs(EPSG:4326) # 导出WKT和属性逗号分隔 gdf[wkt] gdf.geometry.apply(lambda g: g.wkt) output gdf[[name, category_level1, wkt]] output.to_csv(chengyu_poi_wkt.csv, indexFalse, encodingutf-8)WKT里的坐标顺序是“经度 纬度”但很多数据库的PostGIS几何字段默认是“x y”即经度在前。如果你的后端读出来经纬度反了说明表结构里字段顺序定义是lat在前这跟shp无关是你自己建表的锅。我经常收到“POI跑到非洲去了”的技术单十次有九次是经纬度字段反了其次是WKT导出时坐标系被重置成EPSG:3857。shp转TXT的另一种常见需求是转成“lng,lat,name”三列用于热力图绘制。这里有个隐藏坑shp属性表里的中文字段名是UTF-8Excel打开会乱码导出CSV时一定要加encodingutf-8-sig或者干脆用sep\t避免Excel自动转码。4.3 shp转3D Tiles把建筑轮廓拉伸成白模成渝城市群的可视化项目经常要把建筑轮廓shp转成3D Tiles丢给Cesium或Mapbox GL做城市白模。这个流程的核心是把建筑的二维多边形按高度字段拉伸成三维体然后切片。如果你的shp已经有height或楼层字段最简单的方案是用QGIS的“Extrude”工具先拉伸成3D再导出为GeoJSON更工程化的做法是用Python生成每个建筑的顶面和底面多边形合并成三维几何。import geopandas as gpd import numpy as np from shapely.geometry import Polygon, MultiPolygon gdf gpd.read_file(chengyu_building_clean.shp) # 假设shp里有floors字段没有就用高度估算 gdf[height] gdf[floors] * 3.2 # 转成米制投影 gdf_projected gdf.to_crs(EPSG:32648) def extrude(geom, h): if geom.is_empty: return None # 取多边形外环坐标 if geom.geom_type Polygon: rings [list(geom.exterior.coords)] elif geom.geom_type MultiPolygon: rings [list(p.exterior.coords) for p in geom.geoms] else: return None faces [] for ring in rings: # 底面与顶面 bottom Polygon(ring) top Polygon([(x, y, h) for x, y in ring]) faces.append(bottom) faces.append(top) # 侧面 for i in range(len(ring) - 1): p1 ring[i] p2 ring[i 1] quad Polygon([(p1[0], p1[1], 0), (p2[0], p2[1], 0), (p2[0], p2[1], h), (p1[0], p1[1], h)]) faces.append(quad) return faces gdf[extruded] gdf.geometry.apply(lambda g: extrude(g, g[height])) # 炸开所有面 extruded gdf.explode(index_partsTrue) extruded.to_file(chengyu_building_3d.geojson, driverGeoJSON)这段代码生成的是“体素化”的白模侧面只适合展示。真正的3D Tiles切片一般用CesiumLab或py3dtiles把GeoJSON作为输入。注意OSM建筑shp里很多多边形带内部环天井上面的代码只取外环天井会被封死需要特殊处理。另外Cesium对文件大小的要求很苛刻成渝全域建筑shp转GeoJSON动辄2GB我一般先按行政区划分批次导出再在切片端合并或者先简化几何simplify容差0.5米再转。5. 成渝城市群数据整合避坑5个必踩的坑与排查方法5.1 坐标偏移WGS84与GCJ02互转后边界跑到山上现象POI数据导入ArcGIS后点全部落到边界线外侧几百米甚至跑到山体无路区域和道路层完全对不上。原因成渝地区商业地图POI都是GCJ02加密坐标而你用来做底图的行政区划shp是WGS84或CGCS2000。这两套坐标系之间的偏移不是常数成渝地区平均偏移约300米不同乡镇偏移方向各不相同。解决必须使用完整的GCJ02转WGS84算法不能靠平移常数。我建议用开源的coordTransform库或者用geopandas配合转换函数。转换后抽样验证在重庆解放碑、成都天府广场各取10个已知地标点比对转换后坐标与卫星影像位置。如果还有系统偏移考虑是不是把GSC2000当成WGS84用了两者在成渝区域差异约1米到3米做高精度分析时必须区分。5.2 字段乱码dbf编码导致中文全是“锟斤拷”现象shp属性表在ArcGIS里显示正常用geopandas或FME转了一次后中文全部变成“锟斤拷烫烫烫”。原因shp的dbf文件没有内部编码标识。常见编码有GBK和UTF-8两种。你的shp原始编码可能是GBKgeopandas读取时默认识别UTF-8读出来就乱码。还有一种情况是某些工具写shp时强制把字段名转成大写且长度限制10字符导致字段名被截断成“NAME_FIE”。解决读文件时显式指定编码写文件时统一用UTF-8。注意to_file的encoding参数和read_file的encoding参数要保持一致。成渝本地单位交付的shp多为GBK我都会先读一遍看乱码情况再用encodinggbk重新读取后另存为UTF-8。另外不要在shp字段名里用中文统一用拼音或英文字段中文名问题无解规范做法是保存一个字段映射说明文档。5.3 建筑轮廓自相交Union失败、面积算错现象用dissolve合并相邻建筑或者做面积统计时报告“TopologyException: side location conflict”面积值比实际大好几倍。原因OSM建筑多边形普遍存在顶点重复、边界交叉等拓扑错误尤其是手工绘制的轮廓边两条边在交点处交错穿过。解决用第3.2节的make_valid批量修复。修复后一定要用explode拆开否则make_valid可能生成MultiPolygon或GeometryCollection。更彻底的做法是使用PostGIS的ST_MakeValid它能保留更多拓扑关系。成渝数据的另一个特殊现象是建筑轮廓里包含“回”字形结构外环和内环方向相反make_valid不会修复方向如果后续要计算外包围面积需要单独处理内环。5.4 POI经纬度字段反了点和行政区错位现象POI点在成渝边界外用空间连接挂区县时大量POI落入“重庆市-渝北区”但实际应在“四川省-成都市”。导出的WKT里坐标变成“纬度 经度”。原因很多地图API返回的是lat, lng顺序而你写入shp时直接按行解析成lng, lat或者CSV文件的表头是lat,lnggeopandas按列名创建几何时自动交换了顺序但你没注意。解决在转换前统一检查列名强制指定几何的坐标顺序。经过一次错误转换以后虽然坐标值还在合理范围但整片点集被“镜像”到另一个地方。遇到这种情况先用QGIS加载已知地标点比如成都双流机场确认位置再决定是否需要交换坐标或者用gdf.geometry gdf.geometry.apply(lambda g: Point(g.y, g.x))复位。千万别在没验证的情况下直接重新投影。5.5 边界不闭合缓冲区分析结果残缺现象用buffer做道路缓冲区时道路边缘的缓冲区在多边形端点出现“开口”或者clip边界时被裁掉的区域边缘呈锯齿状。原因边界线shp的线段没有闭合即首尾点不一致或者边界要素被切割成多段中间有微小空隙。成渝地区省级边界多为实测折线如果政府发布的整合数据有“简化处理”空隙通常在几十米到几百米。解决先用gdf.geometry gdf.geometry.buffer(0)强制闭合线几何再把多段线line.merge成单要素。对于面状边界用gdf.boundary检查每个多边形是否闭合。做缓冲区分析时我还会对道路先做unary_union消除相邻路段之间的微小间隙避免缓冲区在节点处断开。6. 验证数据质量的一个硬核技巧用Python给成渝城市群shp做质量体检数据整合完不能直接交给下游得做一次系统体检。我习惯写一个质量体检脚本输出每个图层的几何类型、图层范围、属性缺失率、拓扑错误数量并用图形化报告的形式给团队看。核心校验项有四个坐标系是否偏离预设、几何是否有效、要素数量与预期是否匹配、属性表是否有空值。下面是一个可复用的体检代码。import geopandas as gpd import json files { boundary: chengyu_boundary.shp, road: chengyu_road.shp, building: chengyu_building_clean.shp, poi: chengyu_poi_clean.shp, } report {} for name, path in files.items(): gdf gpd.read_file(path, encodingutf-8) invalid_count (~gdf.geometry.is_valid).sum() report[name] { crs: gdf.crs.to_string(), features: len(gdf), invalid_geometries: int(invalid_count), null_attributes: int(gdf.isnull().sum().sum()), bounds: gdf.total_bounds.tolist(), } print(json.dumps(report, indent2, ensure_asciiFalse))跑完看几个关键阈值invalid_geometries必须为0否则后续切片、叠加必翻车bounds要符合成渝城市群经纬度范围约101°E到109°E27°N到33°NPOI图层的null_attributes如果超过5%说明清洗没到位。我觉得值得养成的习惯是把这脚本固定在数据交付前运行像跑单元测试一样。上个月帮朋友复查一份成都三环路POI数据跑完发现建筑层有1200多个无效几何是dwg转换时没有清理的结果用脚本处理后面积误差从15%降到0.3%。数据这行不怕过程绕就怕交付的是黑匣子。做一次体检把坐标系、拓扑、属性检查写进流程后续分析能少熬好几个夜。希望帮到你。本文还有配套的精品资源点击获取
返回列表