九大农业区划SHP处理指南：坐标系、属性表与分区统计

简介「中国九大农业区划shp」数据集面向农业地理、区域规划与GIS空间分析的学习者和研究者将全国九大农业区划以可直接加载的矢量边界呈现便于开展种植结构、土地利用效率与区域差异比较等研究。压缩包共11个文件、约17.21MB以shp矢量主文件为核心配合prj投影信息、dbf属性表、shx索引等GIS配套文件可直接在常见GIS软件中打开并参与制图、叠加分析与空间建模。现有4340人学习下载说明其在教学与科研场景中具有一定参考价值。读者可据此快速搭建九大农业区划底图结合气候、产量或社会经济数据分析各区域作物布局、生产力差异并用于灾害风险评估、农业规划与决策支持等延伸研究。1. 九大农业区划shp 拿到手先判断它能不能直接进你的分析链路下载页上只有一个压缩包名解压出来是一堆同名文件加一个 dbf。第一次碰 shp 中国九大农业区划shp 的人几乎都会在同一个地方卡住双击打不开拖进 QGIS 能看到一堆色块点开属性表全是乱码再一看坐标系是 geographic量出来的面积单位是平方度。这套数据的价值恰恰不在几何而在属性表——九个一级农业区的名称、编码、边界线几何只是壳子。它适合三类人做农业面源统计需要分区骨架的、做栅格分区统计需要掩膜的、做专题图需要底图分区的。但它不能拿来就用坐标系、字段编码、几何有效性这三件事没确认过后面每一次统计都是在错误底图上做乘法。下面按我拆包的实际顺序来从文件骨架讲到分区统计中间该踩的坑一个不落。2. 解开 rar 看骨架shp、shx、dbf、prj、cpg 谁管什么一个 ESRI Shapefile 从来不是单个文件而是一组同名、不同后缀文件的集合。rar 解压后如果只看到.shp先别急着分析多半是解压工具把隐藏文件过滤掉了或者压缩包里本来就没带全。判断这份数据能不能用第一步不是加载而是清点文件。2.1 五个同名文件的职责边界与缺失后果每个后缀承担一件事缺一个就有一类操作直接断掉文件内容缺失后的典型症状.shp几何主体记录每个要素的坐标串完全打不开.shx几何索引记录每个要素在 shp 中的偏移部分软件能打开但翻页极慢或报错.dbf属性表dBASE 格式区划名称和编码都在这能看见图形属性全空.prj坐标系定义WKT 文本图层能显示但坐标被当成无投影的裸坐标.cpg属性表编码声明通常一行UTF-8或GBK中文变问号或方块且不一定报错最容易被忽略的是.cpg。它是一个纯文本小文件内容就一行编码名但很多打包者在压缩时把它漏了于是软件只能靠 dBASE 头里的语言驱动 ID 去猜猜错就是满屏乱码。第二容易被忽略的是.prj——有它和没它geopandas读进来的crs属性一个是正常定义、一个是None后面所有重投影都会因为源坐标系未知而失败。拿到包先做一次文件清点比事后debug省事得多。2.2 用 ogrinfo 和 pyshp 做第一次体检命令行体检优先用 GDAL 自带的ogrinfo它不需要写代码输出信息也够全# -so 只输出摘要不打印每个要素-al 针对所有图层 ogrinfo -so -al nyqh/nyqh.shp # 单独看某个图层的坐标系 WKT gdalsrsinfo nyqh/nyqh.prjogrinfo -so会打印图层名、几何类型、要素数量、字段列表和范围框Extent。重点看三处Geometry是不是Polygon如果是Polygon3D或MultiPolygon混排后面要做统一、Feature Count是否和预期一致、Extent的数值量级——如果 x 在 70 到 140 之间、y 在 15 到 55 之间说明是经纬度如果是六位以上的大数说明已经是投影坐标别重复投影。坐标系确认后再用 Python 做一次字段级检查import shapefile # pyshp # 注意 encoding 参数缺 cpg 时必须手动指定 sf shapefile.Reader(nyqh/nyqh.shp, encodinggbk) # 字段定义[字段名, 类型, 长度, 小数位] for field in sf.fields[1:]: # 第一个元素是 DeletionFlag跳过 print(field) # 看第一条记录的属性值判断中文是否正常 print(sf.record(0)) print(sf.shapeType) # 5 POLYGON, 15 POLYGONZpyshp的encoding参数是关键。不传的话它按 Latin-1 解码中文必定乱码传gbk还是utf-8取决于数据实际编码不确定时两个都试一次看哪个出来的是正常汉字。sf.fields输出的字段名如果出现截断比如QYMC_1、QYMC_2说明原始字段名在导出时超了长度被自动改写这时候必须靠字段顺序而不是字段名去认列这一步后面第 4 章还会细说。2.3 目录组织与只读挂载的习惯我把原始包解压到一个raw/目录并打上只读标记所有重投影、字段清洗、格式转换的输出都写到work/和out/。原因很实际shp 是五六个文件联动的格式用 QGIS 或ogr2ogr做一次就地覆盖写入.shp更新了但.prj没跟上、或.dbf写到一半中断留下一个能打开但属性错位的半残文件这种问题排查起来极其费时。# 原始数据只读防止误写 chmod -R a-w raw/ # 输出目录分开命名带上投影和版本避免覆盖 mkdir -p work out目录分层还有一个好处做重投影时源和目标文件名相同会产生读到自己的输出的风险尤其是批处理脚本里遍历目录时。源在raw/、目标在out/路径不同名基本不会出这类问题。这一条看着像洁癖但只要经历过一次跑完脚本发现数据被自己覆盖了、原文件也没了的情况就会自觉加上。3. 坐标系与投影选型为什么经纬度直接算面积一定会翻车属性表里通常有面积字段但那个字段的来源和算法你并不知道——是投影后算的、还是椭球面算的、单位是平方米还是公顷除非数据说明里写清楚否则不能直接用。更常见的做法是自己重投影后重算。这一步的选型如果错了后面所有统计量都跟着错而且错得很隐蔽地图看起来完全正常只有数字偏了百分之几到十几。3.1 WGS84、CGCS2000 与 Albers 等积投影的选型逻辑先说地理坐标系。EPSG:4326是 WGS84 经纬度EPSG:4490是 CGCS2000 经纬度两者在国内的差异在厘米级做农业区划这种千米级尺度的分析可以认为等价但成果交付时最好按项目要求统一别混着用。判断原数据是哪种看.prj里的DATUM字段WGS_1984还是China_2000。再说投影。全国尺度的面状数据选投影只看一个原则要算面积就用等积投影要算距离就用等距投影两个都要就算两遍。做农业区划统计绝大多数场景是算面积和占比所以选 Albers 等积圆锥投影。中国常用的参数是双标准纬线 25°N 和 47°N中央经线 105°Eprojaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defslat_1和lat_2是两条没有长度变形的标准纬线取 25 和 47 是因为它们大致覆盖了全国纬度范围的中段能把整体面积变形压到很小。lon_0105是中央经线让投影带居中。unitsm决定输出单位是米这一项忘了写算出来的面积数值会离谱。常见的误用有两个。一是直接用EPSG:3857Web 墨卡托算面积这个投影在纬度 40° 以上的面积放大能到两倍以上做全国统计会严重失真它在 Web 地图里好用在统计分析里是错的。二是用 UTM 分带投影拼全国不同带之间拼接处会有明显错位跨带的要素还会被拉变形全国尺度不适用。3.2 用 geopandas 做一次可验证的重投影geopandas的to_crs一行就能转但转完必须验证import geopandas as gpd # 读入时显式指定编码避免 dbf 乱码 gdf gpd.read_file(raw/nyqh.shp, encodinggbk) print(gdf.crs) # 确认源坐标系不是 None print(gdf.shape) # (要素数, 字段数) albers (projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs) # 先做几何有效性修复再重投影顺序不能反 gdf[geometry] gdf[geometry].buffer(0) gdf_aea gdf.to_crs(albers) # 面积用平方米转成平方公里便于阅读 gdf_aea[area_km2] gdf_aea.geometry.area / 1e6 print(gdf_aea[area_km2].sum())逻辑上有三处值得说明。第一buffer(0)是修复自相交多边形的常用手法它把自相交的环重组成合法环比make_valid兼容性好但会丢失一部分原始顶点如果对边界精度敏感就用shapely.make_valid。第二修复放重投影之前因为畸形的几何在投影后可能产生更大的坐标跳变。第三geometry.area只有在投影坐标系下才有物理意义如果gdf_aea.crs还是经纬度这行算出来的是平方度数值小到没法看。验证方法很简单把九个区的面积加总和公开的全国陆地面积做个量级对比。差个百分之几属于边界口径不同是否含岛屿、是否含水域但如果差了百分之几十基本可以确定投影选错或者源坐标系认错了。3.3 ogr2ogr 命令行批处理与投影参数落盘不写代码时用ogr2ogr更省事尤其是要批量转多个文件ogr2ogr -f ESRI Shapefile \ -t_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs \ -lco ENCODINGUTF-8 \ -nlt MULTIPOLYGON \ out/nyqh_aea.shp raw/nyqh.shp-t_srs是目标坐标系-lco ENCODINGUTF-8让输出的 dbf 带上 UTF-8 的 cpg 声明这一步能一次性解决后续在 QGIS 和 Python 里的乱码问题。-nlt MULTIPOLYGON强制统一几何类型因为 shp 一个图层只能存一种几何类型原始数据里如果 Polygon 和 MultiPolygon 混着不指定会在写入时报错或者静默丢要素。输出后去看out/目录应该出现.shp/.shx/.dbf/.prj/.cpg五个文件.prj里是完整的 WKT.cpg里是一行UTF-8缺哪个就回去查对应参数。4. 属性表对齐dbf 中文乱码、字段截断与区划编码归一几何处理完真正耗时间的是属性表。shp 的 dbf 用的是 dBASE 格式这套格式诞生在上世纪八十年代限制很多字段名最长 10 个字节一条记录最长 4000 字节左右字符串字段没有原生 Unicode 支持。中文数据在这些限制下会出现三类问题乱码、字段名被截断重命名、中文字段名导致部分软件读不出来。4.1 dbf 编码的三条线索cpg、LDID 与实测字节判断 dbf 编码按优先级看三条线索。第一条是.cpg文件有它就以它为准内容通常是UTF-8、GBK、936之类。第二条是 dbf 文件头第 29 个字节的 LDIDLanguage Driver ID0x57对应 ANSI/GBK0x00表示未指定未指定时软件会退回系统默认代码页在中文 Windows 上是 GBK在 Linux 容器里往往是 Latin-1于是同一份数据在不同机器上表现不一样这就是典型的我这边好好的那边全是问号。第三条最可靠直接看字节。用十六进制看一眼第一个中文字段的值# 取 dbf 前 400 字节看结构重点看字段值区域的字节 xxd -l 400 nyqh/nyqh.dbf | head -30GBK 的中文是两个字节首字节在0x81-0xFE次字节在0x40-0xFEUTF-8 的中文是三字节形如0xE4-0xE9开头。看一眼就能区分。确认之后统一转成 UTF-8 并写出 cpg后面所有工具都不再需要指定编码import geopandas as gpd # 读入时按实际编码解码 gdf gpd.read_file(raw/nyqh.shp, encodinggbk) # 写出时统一 UTF-8geopandas 会自动生成 cpg gdf.to_file(work/nyqh_utf8.shp, encodingutf-8)这里有个细节read_file的encoding只管 dbf 解码不管.cpgto_file的encoding会同时写入 cpg 文件。所以读的时候指定、写的时候指定这两步都要做只做一半的结果是中间过程正常、落地又乱。4.2 字段名 10 字符上限带来的字段重名dBASE 的字段名上限是 10 个字节。英文名正好 10 个字符中文名在 GBK 下只能放 5 个汉字UTF-8 下只能放 3 个。超了会被截断而且如果截断后重名GDAL 会自动加_1、_2后缀于是QYMC、QYMC_1、QYMC_2同时出现语义完全丢失。遇到这种情况不要靠猜靠字段顺序和样例值反推。先导出前几条记录看值import geopandas as gpd gdf gpd.read_file(work/nyqh_utf8.shp) # 只打印字段名和第一行的值快速对照 for col in gdf.columns: if col ! geometry: print(f{col!r:20} - {gdf[col].iloc[0]!r})字段顺序在 shp 里是稳定的所以你完全可以在原数据说明如果有和样例值的组合下确定哪一列是区划名称、哪一列是区划编码、哪一列是上级区。确定之后立刻改名成稳定的英文短名别在后续脚本里继续用截断名——截断名在不同工具、不同转换轮次下会变脚本今天能跑明天就找不到列。# 按位置重命名比按名字匹配稳 rename_map { gdf.columns[0]: zone_id, # 区划编码 gdf.columns[1]: zone_name, # 区划名称 gdf.columns[2]: zone_level, # 层级 } gdf gdf.rename(columnsrename_map) gdf.to_file(work/nyqh_clean.shp, encodingutf-8)改名前后的对照关系建议单独存一份 CSV写进项目 README。半年后回来接手的人可能就是你自己不用再猜一遍。4.3 区划名称到编码的映射表与合并检查区划名称在不同年份、不同来源的数据里写法会有细微差异——多一个区字、少一个地区、平原写成平原区。做空间连接或属性关联时这些差异会导致部分要素连不上而merge默认是内连接连不上的直接消失你看到的统计结果少了一块却不报错。稳妥做法是先做名称归一再关联最后核对数量import pandas as pd # 归一去空格、去括号内容、统一后缀 def norm(s): if pd.isna(s): return s s str(s).strip().replace( , ) s s.split()[0].split(()[0] return s gdf[join_key] gdf[zone_name].map(norm) # 映射表左连接保留全部区划要素 mapping pd.read_csv(work/zone_code_map.csv) merged gdf.merge(mapping, onjoin_key, howleft, indicatorTrue) # 关键检查有没有没匹配上的 print(merged[_merge].value_counts()) print(merged.loc[merged[_merge] left_only, zone_name].tolist())howleft加indicatorTrue是这套流程里最值钱的两个参数前者保证区划要素不丢后者直接告诉你哪几条没匹配上。打印出来的left_only列表就是需要人工补映射的名单通常不会有几条但每一条都可能导致某个区的统计值缺失。这一步不做最后出图时某个区是空白你还以为是数据里本来就没有。5. 避坑与排查五条血泪记录这一章的五条都是我自己踩过或者帮别人排过的症状都不明显共同点是看起来正常数字不对。每条按现象、原因、解决写能对号入座就直接跳过排查过程。5.1 五条具体记录第一条QGIS 里中文全是问号Python 里正常。现象是同一份 shpgeopandas读出来是汉字拖进 QGIS 显示????。原因是 dbf 里没有.cpgQGIS 按自身设置的系统编码去解而读数据的机器编码设置不同。解决是重写一遍带上 cpggdf.to_file(..., encodingutf-8)或者手工建一个同名.cpg文本文件写入UTF-8。注意别用编辑器保存成带 BOM 的 UTF-8某些 GDAL 版本会把 BOM 当成编码名的一部分。第二条面积加起来比实际大了一截。现象是九个区面积合计明显偏大比例在百分之十到几十之间。原因是用 Web 墨卡托或经纬度直接算的。解决是先重投影到 Albers 等积再算并且确认crs的单位是米不是度。判断方法很直接gdf.crs.axis_info[0].unit_name打印出来是metre才算对。第三条某个区在做空间连接后消失了。现象是连接完要素数量比原来少。原因是属性表里的区划名称有不可见字符尾部空格、全角空格、换行或者用了howinner。解决是连接前统一strip并去掉全角空格\u3000连接后用indicatorTrue检查left_only。全角空格这一条尤其隐蔽肉眼完全看不出来len()长度却对不上。第四条裁剪后边界出现细缝或重叠。现象是用区划做掩膜裁剪栅格或者用区划做叠加分析相邻区之间出现一两像素的缝。原因是相邻多边形共享边界但顶点不重合这是不同批次数据拼接的常见结果。解决是裁剪前对整体做一次节点对齐shapely里可以用unary_union后重新切分或者在 GIS 里跑一次拓扑修复。对精度要求不高的场景栅格统计时用all_touchedTrue能掩盖大部分缝但矢量的叠加分析掩盖不了。第五条换台机器后脚本报字段不存在。现象是本地跑通的脚本放到服务器上读 shp 时报列名不存在。原因是字段名被截断成了QYMC_1这种带后缀的形式而不同 GDAL 版本对重名后缀的生成规则不完全一致。解决是读取后立刻按位置重命名成固定英文名脚本里永远不出现原始中文字段名。这条是从根上断掉兼容性问题比在每台机器上对齐 GDAL 版本省事。5.2 一套最小自检脚本每次拿到新的 shp不管是这份区划还是别的我都会先跑一遍这个脚本几十秒的事能提前发现八成问题import geopandas as gpd def health_check(path, encodingutf-8): gdf gpd.read_file(path, encodingencoding) print(要素数:, len(gdf)) print(几何类型:, gdf.geom_type.value_counts().to_dict()) print(坐标系:, gdf.crs.to_string() if gdf.crs else None) print(无效几何:, (~gdf.geometry.is_valid).sum()) print(空几何:, gdf.geometry.is_empty.sum()) print(重复几何:, gdf.geometry.duplicated().sum()) # 字段名长度检查超过 10 字节的会在下次导出时被截断 for c in gdf.columns: if c ! geometry and len(c.encode(utf-8)) 10: print(超长字段名:, c) health_check(work/nyqh_clean.shp)六个指标里无效几何和重复几何是要处理的坐标系为 None是要补的超长字段名是要改的。空几何一般说明源数据里有占位要素视需求删掉。5.3 出问题时的排查顺序顺序定死能省时间先看文件齐不齐再看坐标系有无再看几何有效性最后才看属性。反过来查最常见的结果是在属性表里翻半天最后发现是.prj丢了导致整体偏移。另一个习惯是把中间结果落盘——每一步转换后的 shp 都保留出问题时能快速定位是哪一步引入的而不是从头重跑一遍猜。磁盘空间远比排查时间便宜。6. 进阶用法把区划当分区骨架做栅格统计与成果核验矢量本身只是骨架真正的产出通常是把栅格数据按区聚合。比如手上有一份 NDVI 或降水栅格想算九个农业区各自的均值用rasterstats的zonal_stats是常规做法。这一节写清楚参数怎么选、结果怎么验。6.1 zonal_stats 的 all_touched 与 nodata 选择from rasterstats import zonal_stats stats zonal_stats( work/nyqh_clean.shp, # 矢量必须是投影坐标系或与栅格一致 work/ndvi.tif, # 栅格 stats[mean, count, std], # 需要的统计量 all_touchedFalse, # 只统计像元中心落在面内的 nodata-9999, # 显式指定别依赖 tif 自带的 geojson_outFalse, ) # 结果是一个列表顺序与矢量要素顺序严格对应 import pandas as pd df pd.DataFrame(stats) print(df.head())all_touched是这套流程里最需要想清楚的参数。False时只有像元中心落在多边形内的才计入面积统计偏保守边界处会丢半个像元True时只要像元与多边形有任何接触就计入相邻两个区会重复统计边界像元导致各区之和大于全域总量。做区际对比用False做全域覆盖估算用True两个都跑一遍看差异有多大差异超过百分之五就说明你的区划边界和栅格分辨率不匹配该考虑重采样了。nodata必须显式指定。tif 自带的 nodata 值有时是-3.4e38这种极端值如果不声明它会被当成真实数值参与均值计算结果直接崩掉。这一点用打印核对最快统计结果里count应该接近该区面积除以像元面积差太多就是 nodata 没处理好。结果回填到矢量上再做图顺序不能反gdf gpd.read_file(work/nyqh_clean.shp) gdf[ndvi_mean] [s[mean] for s in stats] gdf.to_file(out/nyqh_stats.shp, encodingutf-8)6.2 用面积占比反查区划边界是否可靠统计做完别急着出图先做一次反向核验把每个区算出来的像元数乘以像元面积和该区的几何面积比一比比值应该接近 1。偏差超过一成通常是三个原因之一——坐标系不一致矢量和栅格不在同一投影下zonal_stats不会报错但结果全错、区划边界和栅格范围没对齐、或者有要素的几何是无效的。还有一个更粗但有效的检查九个区的面积占比加总必须是 100%。如果加出来是 103%说明区划之间存在重叠这在多来源拼接的数据里不罕见。处理方式是在做统计前先做一次unary_union再切分或者接受重叠但在报告里注明口径。我一般会把这项检查写成脚本里的断言跑不通就中止不让自己带着错误结果往下走。从那以后我每次拿到新的区划 shp都强制走一遍文件清点、坐标系确认、几何修复、字段重命名、面积占比核验这五步一步不跳。跳过哪一步后面就得花几倍时间回头找。希望帮到你。本文还有配套的精品资源点击获取