ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

粤港澳大湾区shp数据实战:坐标系校验、格式转换与落地避坑指南

粤港澳大湾区shp数据实战:坐标系校验、格式转换与落地避坑指南 简介这份粤港澳大湾区shp数据面向从事GIS制图、空间分析与区域规划的研究人员、学生及从业者提供开箱即用的矢量底图素材可用于城市群边界绘制、专题地图制作与空间可视化练习。资源包共20个文件约272KB以shp、shx、dbf、prj等Shapefile核心文件为主分别承载几何图形、索引、属性表与坐标投影信息另含sbx、sbn空间索引及xml元数据覆盖粤港澳大湾区整体范围与hk_line、Macon_line等线要素图层便于直接加载到ArcGIS、QGIS等平台进行叠加分析与出图。目前已有291人学习下载适合需要快速获取区域边界数据、搭建制图底图或开展空间统计的读者参考使用能有效减少数据搜集与格式转换的时间成本。1. 粤港澳大湾区 shp 数据从哪拿、怎么用、坐标系怎么不翻车做城市分析、选址评估、物流路径规划的人迟早会撞上同一个需求要一份能直接扔进 GIS 或代码里跑的粤港澳大湾区 shp 数据。它通常包含大湾区 11 个城市广州、深圳、珠海、佛山、惠州、东莞、中山、江门、肇庆、香港、澳门的行政区划边界有时还带区县一级、路网、水系、建成区。拿到手之后你能干的事很具体按城市裁剪遥感影像、统计各区县面积、做渔网分割 shp 后抽样、把边界转成 3dtiles 做三维底图、导出 wkt 塞进数据库做空间查询。但真正卡住人的从来不是「有没有数据」而是「这份数据能不能用」。坐标系是 4490 还是 4326边界有没有拓扑错误属性表里城市名是中文还是拼音直接决定你后面会不会返工。这篇就按一线做法把获取、校验、转换、落地整条链路讲清楚新手能照着跑熟手能对着参数挑刺。2. 拿到 shp 之后先别急着画图坐标系与属性校验2.1 为什么 4490 和 4326 混用会让你白干一天粤港澳大湾区 shp 数据最常见的两种坐标系是 CGCS2000EPSG:4490和 WGS84EPSG:4326。两者在米级精度上差异很小很多人觉得「差不多随便用」结果在叠加分析时发现边界错位几十米尤其是跨珠江口做缓冲区分析时错位会被放大。判断方法很简单用 Python 读一下.prj文件或者直接看 crsimport geopandas as gpd # 读取 shp注意 encoding 参数中文属性表常用 gbk 或 utf-8 gdf gpd.read_file(greater_bay_area.shp, encodingutf-8) # 打印坐标系和范围 print(CRS:, gdf.crs) print(Bounds:, gdf.total_bounds) print(要素数量:, len(gdf)) print(字段列表:, list(gdf.columns))逻辑说明gpd.read_file会自动解析.prj如果返回None说明这份 shp 缺投影文件必须手动指定。total_bounds返回[minx, miny, maxx, maxy]大湾区大致落在东经 111°~115°、北纬 21°~24° 之间如果数值是几百万甚至上千万那多半是投影坐标如 UTM 或高斯克吕格不是经纬度。参数说明encoding是血泪经验点。很多从 ArcGIS 导出的 shp 属性表用 GBK直接读会报UnicodeDecodeError换成gbk即可。如果字段名出现乱码优先怀疑编码而不是数据损坏。提示如果 crs 为 None 且坐标范围在 111~115 之间直接gdf.set_crs(EPSG:4490, inplaceTrue)不要用 4326 硬套后续要转再转。2.2 属性表里藏着的三个坑空值、重复、命名不一致拿到一份县域行政区划边界 shp第一件事不是画图是看属性表。常见问题有三个第一城市名字段有空值或空白字符串导致按城市分组时丢数据。第二同一城市出现多条记录比如飞地或重复录入做面积统计时翻倍。第三命名不统一有的写「广州市」有的写「广州」有的写拼音「Guangzhou」。# 检查空值和重复 print(gdf[city_name].isna().sum()) print(gdf[city_name].value_counts()) # 统一命名去掉「市」字后缀便于后续匹配 gdf[city_std] gdf[city_name].str.replace(市, , regexFalse).str.strip() # 检查几何有效性 invalid gdf[~gdf.geometry.is_valid] print(无效几何数量:, len(invalid))逻辑说明is_valid能查出自我相交、环方向错误等拓扑问题。无效几何在做intersection、union时会直接抛异常或返回空结果。修复用gdf.geometry gdf.geometry.buffer(0)这是常用急救手段但会轻微改变边界精度要求高时慎用。参数说明buffer(0)对大多数自相交多边形有效但对完全退化的线状多边形无效那种情况只能手动删记录。value_counts()输出里如果某个城市数量明显偏多基本就是重复录入。3. 从 DWG、DEM、渔网到 shp几种常见转换的实操路径3.1 dwg 转 shp图层筛选和闭合检查是成败关键规划院给的基础图经常是 DWG要转成 shp 才能做空间分析。常见做法是用 ArcGIS 的「CAD 至地理数据库」工具或者用 ODA File Converter 先转 DXF 再用 GDAL 读。我一般走 GDAL 路线因为可脚本化。# 先用 ogrinfo 看 DWG/DXF 里有哪些图层 ogrinfo -so input.dxf # 只导出需要的图层到 shp注意 -nlt POLYGON 强制多边形 ogr2ogr -f ESRI Shapefile output.shp input.dxf \ -sql SELECT * FROM boundary_layer \ -nlt POLYGON \ -a_srs EPSG:4490逻辑说明DWG 里线和面混在一起-nlt POLYGON会尝试把闭合线转成面但前提是线必须真正闭合。如果导出后要素数量对不上多半是线没闭合。-a_srs指定输出坐标系DWG 本身通常不带坐标系信息必须手动给。参数说明-sql里的图层名要和ogrinfo输出一致大小写敏感。如果 DWG 版本太新 GDAL 读不了先用 ODA 转成 DXF R12 或 2000 格式。注意DWG 转出来的面经常有自相交转完必须跑一遍is_valid检查别直接进分析流程。3.2 从 DEM 提取 shp等值线和流域边界的两种玩法从 DEM 提 shp 有两个典型需求一是提取等高线做地形图二是提取流域或汇水区边界。等高线用 GDAL 的gdal_contour一条命令搞定# 从 DEM 提取每 10 米一条的等高线 gdal_contour -a elev -i 10 input_dem.tif contour.shp逻辑说明-a elev把高程值写进属性表字段elev-i 10是等高距。生成的 shp 是线要素如果要面还得再做一次闭合处理通常没必要。流域边界复杂一些需要先用gdaldem算流向和流量累积再阈值提取河网最后用r.water.outletGRASS或 WhiteboxTools 做汇水区。这条链路参数多新手建议先用 QGIS 的「流域分割」插件跑通一次再考虑脚本化。参数说明等高距-i要根据 DEM 分辨率选30 米 DEM 用 10 米等高距会非常密建议 20~50 米。DEM 如果有 NoData 空洞gdal_contour会在空洞边缘生成异常线先用gdal_fillnodata补洞。3.3 渔网分割 shp做空间抽样的标准前置步骤渔网分割 shp 是空间抽样、格网统计的常用手段。比如你要在大湾区范围内均匀抽 500 个点做实地调查直接随机撒点可能落到海里用渔网先切再筛选就稳得多。import geopandas as gpd from shapely.geometry import box import numpy as np gdf gpd.read_file(gba_boundary.shp) minx, miny, maxx, maxy gdf.total_bounds # 创建 0.1 度间隔的渔网 cell_size 0.1 cols np.arange(minx, maxx, cell_size) rows np.arange(miny, maxy, cell_size) cells [box(x, y, x cell_size, y cell_size) for x in cols for y in rows] grid gpd.GeoDataFrame(geometrycells, crsgdf.crs) # 只保留落在边界内的格子 grid_clipped gpd.overlay(grid, gdf, howintersection) print(有效格子数:, len(grid_clipped)) grid_clipped.to_file(gba_grid.shp, encodingutf-8)逻辑说明box生成矩形overlay做相交裁剪只保留与大湾区边界有交集的格子。cell_size用度为单位0.1 度约 11 公里做城市级抽样够用做区县级要降到 0.01 度。参数说明overlay的how参数选intersection会保留裁剪后的几何如果只想要完整格子改用grid.sjoin(gdf, predicateintersects)再筛。数据量大时overlay很慢可以先做 bounding box 粗筛。4. shp 数据落地避坑编码、拓扑、导出格式的常见翻车4.1 中文乱码从 shp 到 GeoJSON 到数据库的编码链路现象shp 在 QGIS 里显示正常导出 GeoJSON 后中文变问号再入库变成乱码。原因shp 的.dbf文件编码没有统一标准ArcGIS 默认用系统编码中文 Windows 是 GBK而 GeoJSON 规范要求 UTF-8。中间任何一步没转编码就会断链。解决读取时显式指定encodinggbk写出时强制encodingutf-8。如果已经乱码用gdf[name] gdf[name].str.encode(latin1).str.decode(gbk)尝试修复但只对特定乱码模式有效。# 读 GBK写 UTF-8一步到位 gdf gpd.read_file(input.shp, encodinggbk) gdf.to_file(output.geojson, driverGeoJSON, encodingutf-8)4.2 拓扑错误自相交和重叠面让叠加分析直接崩现象做intersection时程序卡死或报TopologyException。原因多边形自相交、相邻面共享边但节点不重合、面与面重叠。这些在手工数字化数据里非常常见。解决先is_valid定位再buffer(0)修复修复后重新检查。如果buffer(0)后仍有无效几何用make_validShapely 1.8。from shapely.validation import make_valid gdf[geometry] gdf[geometry].apply( lambda g: make_valid(g) if not g.is_valid else g )4.3 导出格式选错shp 的字段名长度和类型限制现象导出 shp 后字段名被截断成 10 个字符或者数值字段变成文本。原因ESRI Shapefile 格式本身限制字段名最长 10 字符且不支持某些数据类型如 64 位整数、日期时间带时区。解决如果字段名重要改用 GeoPackage.gpkg或 GeoJSON。如果必须用 shp提前把字段名缩到 10 字符以内数值统一用 float。4.4 坐标系丢失转了一圈回来发现位置偏了现象shp 经过几次转换后叠加到在线底图上位置偏移。原因中间某一步用了set_crs而不是to_crs或者导出时没带.prj文件。解决set_crs只改标签不改坐标to_crs才做实际转换。转换前确认源坐标系转换后立刻验证total_bounds是否合理。# 正确从 4490 转到 4326 gdf_4326 gdf.to_crs(EPSG:4326) # 错误把 4490 的数据硬标成 4326坐标值没变 # gdf_wrong gdf.set_crs(EPSG:4326)4.5 大数据量卡顿4 万条以上要素的处理策略现象读取或写出几万条要素的 shp 时内存爆掉或耗时过长。原因shp 格式本身不适合大数据量且 GeoPandas 默认全量加载。解决用pyogrio引擎替代默认的 Fiona读取速度能快几倍或者分块处理用bbox参数只读当前范围。# 用 pyogrio 引擎加速读取 gdf gpd.read_file(large.shp, enginepyogrio, use_arrowTrue) # 只读指定范围 gdf_sub gpd.read_file(large.shp, bbox(113.0, 22.0, 114.0, 23.0))5. 把 shp 用起来转 3dtiles、导 WKT、接数据库的进阶技巧数据校验完、转换完最后一步是让它真正进入你的系统。这里挑三个高频场景讲。转 3dtiles 做三维底图shp 本身是二维的要转 3dtiles 需要先给高度。常见做法是用gdal_rasterize把 shp 烧录成栅格再用gdal_translate转成带高程的 GeoTIFF最后用 CesiumLab 或py3dtiles切片。如果只是要平面贴地直接用ogr2ogr转成 GeoJSON 再喂给 Cesium 的GeoJsonDataSource更省事。导出 WKT 入库PostGIS 和很多空间数据库都支持 WKT。用 GeoPandas 一行就能导出# 导出 WKT 和属性到 CSV方便入库 gdf[wkt] gdf.geometry.apply(lambda g: g.wkt) gdf.drop(columnsgeometry).to_csv(gba_wkt.csv, indexFalse, encodingutf-8)入库时用ST_GeomFromText(wkt, 4490)构造几何注意 SRID 要和数据一致。如果 CSV 里 WKT 太长检查是不是用了g.wkt而不是g.simplify(0.001).wkt后者能在不影响可视化的前提下大幅缩短字符串。接 Python 做自动化把上面所有步骤串成一个脚本输入一个 shp输出校验报告 转换后的 GeoPackage WKT CSV。我一般会加一个--city参数只处理指定城市避免每次全量跑。import argparse parser argparse.ArgumentParser() parser.add_argument(--input, requiredTrue) parser.add_argument(--city, defaultNone) args parser.parse_args() gdf gpd.read_file(args.input, encodinggbk) if args.city: gdf gdf[gdf[city_std] args.city] # 后续校验、转换、导出...最后一个习惯每次拿到新的 shp先跑一遍is_valid、crs、total_bounds、value_counts四件套再动手做任何分析。这个习惯帮我省掉的返工时间比任何优化技巧都多。希望帮到你。本文还有配套的精品资源点击获取
返回列表