ARTICLE DETAIL

资讯详情

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

武威市区县边界shp文件使用指南:坐标系、投影转换与空间分析实战

武威市区县边界shp文件使用指南:坐标系、投影转换与空间分析实战 简介这份甘肃武威市区县级别行政区划shp文件面向GIS从业者、城市规划人员及地理信息相关专业师生用于区域分析、人口统计、交通规划与自然资源管理等空间分析场景。资源包共21个文件约338KB以shp、shx、dbf、prj等Shapefile核心格式为主分别承载几何边界、索引、属性表与坐标投影信息另含sbn、sbx、cpg、xml等辅助文件可直接在ArcGIS、QGIS等软件中加载使用。数据覆盖武威市下辖凉州区、民勤县、古浪县和天祝藏族自治县包含区县边界线、中心点坐标及名称代码等属性属于GIS应用中的基础数据层。目前已有454人学习下载适合需要精确行政边界底图、开展矢量地图制作与空间查询分析的读者参考使用。1. 武威市区县边界 shp 文件从拿到手到落进 GIS 工程里做河西走廊相关项目的人多半绕不开武威。不管是做农业种植区划、光伏选址、还是做行政区划底图第一步往往就是找一份靠谱的武威市区县级别行政区划 shp 文件。我手上这份就是武威市下辖一区两县一自治县的矢量边界凉州区、民勤县、古浪县、天祝藏族自治县四个县级单元面要素WGS84 地理坐标系。它解决的核心问题很直接——让你在 QGIS、ArcGIS 或者 PostGIS 里有一个能直接参与空间运算、能做裁剪、能做叠加分析的行政边界底图而不是从网页地图上截图描边。适合谁做区域统计的、做专题制图的、做选址分析的以及需要把统计年鉴数据挂到空间单元上的从业者。下面我按拿到文件之后实际会走的流程拆一遍包括坐标系怎么核、属性表怎么改、和统计数据处理时怎么对齐以及几个我踩过的坑。2. 先搞清 shp 文件到底给了你什么坐标系、字段与几何类型2.1 四个文件缺一不可少一个都打不开shp 不是单个文件是一组。拿到手先看目录里有没有这四个文件后缀作用缺失后果.shp几何坐标数据没有几何图层为空.shx几何索引能打开但无法查询、无法编辑.dbf属性表能显示图形但无字段、无法分类.prj坐标系定义坐标系未知面积量算全错常见做法是拿到压缩包先解压到一个纯英文路径下别放在桌面或者带中文和空格的目录里。ArcGIS 对中文路径的容忍度比 QGIS 低尤其是做投影转换的时候路径里带中文偶尔会报 000732 错误。我一般会建一个D:\gis_data\wuwei_admin这样的目录把四个文件放进去再在 QGIS 里拖 .shp 进去看。2.2 坐标系确认别默认它就是 CGCS2000这份文件我打开后 .prj 里写的是 GCS_WGS_1984也就是 EPSG:4326。但这里有个血泪经验很多网上流传的行政区划 shpprj 文件写的是 WGS84实际坐标却是 GCJ02 偏移过的或者干脆是某个地方坐标系。验证方法很简单在 QGIS 里叠加一个在线底图比如 OSM 或者天地图影像看边界和影像上的行政界线是否重合。如果整体偏移几百米那就是坐标系标错了。# 用 ogrinfo 快速查看 shp 的坐标系和字段信息 ogrinfo -al -so wuwei_admin.shp # 输出里重点看这两行 # Geometry: Polygon # SRS: projlonglat datumWGS84 no_defsogrinfo是 GDAL 自带的命令行工具-al表示列出所有要素-so表示只输出摘要不输出具体几何。如果你看到 SRS 那一行是datumWGS84说明文件自认是 WGS84。但自认不等于正确还是要叠底图验证。如果要做面积统计WGS84 地理坐标系下算出来的面积是平方度没有意义必须转成投影坐标系。2.3 属性表里该有什么区县名称和编码打开属性表至少应该有区县名称字段。我手上这份的字段结构大致是字段名类型示例值说明NAMEString凉州区区县中文名ADCODEString620602行政区划代码CITYString武威市所属地市如果只有 NAME 没有 ADCODE做数据关联的时候会很难受。因为统计年鉴里的数据往往用行政区划代码做键中文名会有「天祝藏族自治县」和「天祝县」这种不一致。常见做法是手动补一列 ADCODE武威市下辖的代码是凉州区 620602、民勤县 620621、古浪县 620622、天祝藏族自治县 620623。补完之后用 ADCODE 做连接比用名称稳得多。# 用 geopandas 检查属性表并补全行政区划代码 import geopandas as gpd gdf gpd.read_file(wuwei_admin.shp) print(gdf.columns.tolist()) print(gdf[[NAME]].to_string()) # 如果缺 ADCODE按名称映射补上 code_map { 凉州区: 620602, 民勤县: 620621, 古浪县: 620622, 天祝藏族自治县: 620623 } gdf[ADCODE] gdf[NAME].map(code_map) gdf.to_file(wuwei_admin_fixed.shp, encodingutf-8)这段代码先读入 shp打印字段名和名称列确认四个区县都在。然后用字典映射补 ADCODE。to_file的时候指定encodingutf-8否则中文在 dbf 里可能变成乱码。注意 geopandas 写 shp 时默认会用 UTF-8 编码 dbf但 ArcGIS 老版本读 UTF-8 的 dbf 可能显示乱码如果下游是 ArcGIS可以改成encodinggbk试试。3. 把 shp 用起来投影转换、裁剪与统计挂接3.1 投影转换算面积之前必须做的一步WGS84 下直接算面积得到的是平方度数值没有实际意义。要做面积统计或者缓冲区分析先转成投影坐标系。武威市位于甘肃中部经度大约在东经 101° 到 104° 之间横跨 3 度带和 6 度带的多个带。常见做法是用 CGCS2000 的 3 度带投影中央经线选 102°E 或者 105°E。import geopandas as gpd gdf gpd.read_file(wuwei_admin_fixed.shp) # 转成 CGCS2000 3度带中央经线102°EEPSG:4544 gdf_proj gdf.to_crs(epsg4544) # 算每个区县的面积单位平方米 gdf_proj[area_km2] gdf_proj.geometry.area / 1e6 print(gdf_proj[[NAME, area_km2]].to_string()) gdf_proj.to_file(wuwei_admin_proj.shp, encodingutf-8)to_crs(epsg4544)是 CGCS2000 / 3-degree Gauss-Kruger CM 102E适用于武威大部分区域。如果你做的是全市尺度的分析跨了多个投影带可以考虑用 Albers 等积投影自定义中央经线和标准纬线。area / 1e6是把平方米转成平方公里。算完之后可以对照公开数据核一下凉州区面积大约 4900 多平方公里民勤县约 1.58 万平方公里古浪县约 5100 多平方公里天祝县约 7100 多平方公里。如果算出来差一个数量级大概率是投影没转对。3.2 用区县边界裁剪其他数据这是最常见的用法你有一份武威市的土地利用栅格、或者一份 POI 点数据想按区县拆分统计。用 geopandas 的clip或者overlay都能做。import geopandas as gpd import pandas as pd # 读入区县边界和待裁剪的点数据 admin gpd.read_file(wuwei_admin_proj.shp) poi gpd.read_file(wuwei_poi.shp) # 确保两者坐标系一致 poi poi.to_crs(admin.crs) # 按区县做空间连接给每个点打上所属区县标签 poi_joined gpd.sjoin(poi, admin[[NAME, geometry]], howleft, predicatewithin) # 统计每个区县的点数量 count_by_county poi_joined.groupby(NAME).size().reset_index(namepoi_count) print(count_by_county) # 如果要按区县导出独立的 shp for name in admin[NAME].unique(): county admin[admin[NAME] name] clipped gpd.clip(poi, county) clipped.to_file(fpoi_{name}.shp, encodingutf-8)sjoin的predicatewithin表示点落在面内howleft保留所有点落在边界外的点 NAME 会是空。groupby(NAME).size()就是按区县计数。最后那个循环是按区县导出注意文件名里如果带中文在某些系统上会有问题可以换成 ADCODE 做文件名。裁剪的时候如果数据量大gpd.clip会比sjoin慢但结果更干净。3.3 和统计年鉴数据挂接统计年鉴里通常是 Excel有行政区划代码和一堆指标。用 pandas 读进来和 shp 的属性表按 ADCODE 合并再写回 shp就能在 QGIS 里按指标做分级设色。import pandas as pd import geopandas as gpd # 读统计年鉴数据 stats pd.read_excel(wuwei_stats_2023.xlsx, dtype{ADCODE: str}) # 读 shp gdf gpd.read_file(wuwei_admin_fixed.shp) gdf[ADCODE] gdf[ADCODE].astype(str) # 合并 merged gdf.merge(stats, onADCODE, howleft) # 检查有没有没匹配上的 print(merged[merged[指标列名].isna()][[NAME, ADCODE]]) merged.to_file(wuwei_admin_stats.shp, encodingutf-8)关键点是dtype{ADCODE: str}因为 Excel 里的代码如果被当成数字读前导零会丢620602 变成 620602 还好但如果是 0620602 这种就会出问题。合并之后一定要检查有没有 NaN如果有说明代码对不上可能是统计年鉴里用了老的代码或者名称写法不同。4. 避坑与排查shp 文件处理中最容易翻车的五个点4.1 中文乱码dbf 编码不对现象在 QGIS 里打开属性表区县名称显示成乱码或者问号。原因dbf 文件的编码和软件读取时用的编码不一致。老版本的 shp 常用 GBK 或者 GB2312新写的常用 UTF-8。解决在 QGIS 里右键图层 → 属性 → 源看编码设置手动切成 GBK 或者 UTF-8 试。如果用 geopandas 写文件明确指定encodingutf-8或encodinggbk。ArcGIS 用户注意ArcGIS Pro 对 UTF-8 支持较好ArcMap 可能需要用 GBK。4.2 坐标系标错边界和底图对不上现象把 shp 叠到在线底图上整体偏移几百米到几公里。原因prj 文件写的坐标系和实际坐标不符常见的是 WGS84 和 GCJ02 混用。解决先叠底图目视检查如果偏移是系统性的用 QGIS 的「矢量 → 数据管理工具 → 重新投影图层」试几个候选坐标系看哪个能对上。如果偏移量不固定可能是数据本身被做过非线性变换那就只能找原始来源。4.3 几何无效做叠加分析时报错现象做overlay或者intersection的时候报 TopologyException 或者 Geometry invalid。原因面要素有自相交、悬挂节点或者重复点。解决在 QGIS 里用「矢量 → 几何工具 → 检查有效性」或者用 geopandas 的buffer(0)做一次修复。gdf[geometry] gdf[geometry].buffer(0)buffer(0)是常见的几何修复技巧能把大部分自相交问题消掉。如果还不行用make_valid。4.4 面积算出来是负数或者零现象geometry.area返回负值或者极小值。原因要么是坐标系没转要么是几何方向反了外环顺时针。解决先确认转到了投影坐标系然后用shapely的orient把几何方向统一。from shapely.geometry import polygon from shapely.ops import orient gdf[geometry] gdf[geometry].apply(lambda g: orient(g, sign1.0))sign1.0表示外环逆时针这是 OGC 标准方向大部分空间运算库都认这个。4.5 字段名截断dbf 字段名最长 10 个字符现象写 shp 的时候字段名被截断比如area_km2变成area_km2还好但population_density会变成populati。原因dbf 格式限制字段名最多 10 个字符。解决写文件之前把字段名改短或者用 GeoPackage 格式代替 shpgpkg 没有这个限制。gdf gdf.rename(columns{population_density: pop_dens}) gdf.to_file(output.gpkg, driverGPKG)如果下游必须用 shp那就接受截断或者在属性表里加一个说明文档。5. 进阶技巧用区县边界做空间插值和分区统计5.1 把气象站点数据插值到区县假设你有一份武威市及周边气象站点的气温数据想算每个区县的平均气温。思路是先用站点做空间插值生成栅格再用区县边界做分区统计。import geopandas as gpd import rasterio from rasterstats import zonal_stats # 读区县边界 admin gpd.read_file(wuwei_admin_proj.shp) # 假设已经用站点数据插值生成了气温栅格 temperature.tif # 做分区统计 stats zonal_stats( admin, temperature.tif, stats[mean, min, max], geojson_outTrue ) # 把结果转成 GeoDataFrame result gpd.GeoDataFrame.from_features(stats) print(result[[NAME, mean, min, max]])zonal_stats的第一个参数是面矢量第二个是栅格路径stats指定要算的统计量。geojson_outTrue让结果带几何方便后续导出。注意栅格和矢量的坐标系必须一致不一致的话先用rasterio的warp或者 QGIS 的「对齐栅格」工具处理。5.2 用区县边界做选址分析比如你要找武威市适合建光伏电站的区域条件包括坡度小于 15 度、不在基本农田内、距离道路不超过 2 公里、属于某个区县。这种多条件叠加用 QGIS 的「按位置选择」或者 PostGIS 的ST_Intersects都能做。我一般会在 PostGIS 里做因为数据量大时 SQL 比桌面软件稳。-- 假设有区县边界表 wuwei_admin、坡度栅格转成的矢量表 slope、道路表 roads -- 找出凉州区内坡度小于15度且距离道路2公里内的区域 SELECT a.name, ST_Union(ST_Intersection(a.geom, s.geom)) AS suitable_geom FROM wuwei_admin a JOIN slope s ON ST_Intersects(a.geom, s.geom) JOIN roads r ON ST_DWithin(a.geom, r.geom, 0.02) -- 0.02度约2公里 WHERE a.name 凉州区 AND s.slope 15 GROUP BY a.name;ST_DWithin的距离单位取决于坐标系如果是 WGS84 地理坐标系单位是度0.02 度大约 2 公里但随纬度变化。更稳妥的做法是转成投影坐标系再用米做单位。ST_Union把符合条件的碎片合并成一个大多边形。5.3 一个我常犯的错误早些年我做分区统计的时候直接拿 WGS84 的 shp 去和投影坐标系下的栅格做zonal_stats结果出来一堆 NaN。排查了半天才发现是坐标系不一致。从那以后我每次拿到 shp第一件事就是ogrinfo -al -so看 SRS第二件事是叠底图目视检查第三件事是确认目标分析用的坐标系三步走完再动手。这个习惯帮我省了很多后悔药。希望帮到你。本文还有配套的精品资源点击获取
返回列表