简介：「中国九大农业区划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 五个同名文件的职责边界与缺失后果

每个后缀承担一件事，缺一个就有一类操作直接断掉：

最容易被忽略的是.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.prj

ogrinfo -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", encoding="gbk") # 字段定义：[字段名, 类型, 长度, 小数位] for field in sf.fields[1:]: # 第一个元素是 DeletionFlag，跳过 print(field) # 看第一条记录的属性值，判断中文是否正常 print(sf.record(0)) print(sf.shapeType) # 5 = POLYGON, 15 = POLYGONZ

pyshp的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°E：

+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

lat_1和lat_2是两条没有长度变形的标准纬线，取 25 和 47 是因为它们大致覆盖了全国纬度范围的中段，能把整体面积变形压到很小。lon_0=105是中央经线，让投影带居中。units=m决定输出单位是米，这一项忘了写，算出来的面积数值会离谱。

常见的误用有两个。一是直接用EPSG:3857（Web 墨卡托）算面积，这个投影在纬度 40° 以上的面积放大能到两倍以上，做全国统计会严重失真，它在 Web 地图里好用，在统计分析里是错的。二是用 UTM 分带投影拼全国，不同带之间拼接处会有明显错位，跨带的要素还会被拉变形，全国尺度不适用。

3.2 用 geopandas 做一次可验证的重投影

geopandas的to_crs一行就能转，但转完必须验证：

import geopandas as gpd # 读入时显式指定编码，避免 dbf 乱码 gdf = gpd.read_file("raw/nyqh.shp", encoding="gbk") print(gdf.crs) # 确认源坐标系不是 None print(gdf.shape) # (要素数, 字段数) albers = ("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 " "+x_0=0 +y_0=0 +datum=WGS84 +units=m +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 "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" \ -lco ENCODING=UTF-8 \ -nlt MULTIPOLYGON \ out/nyqh_aea.shp raw/nyqh.shp

-t_srs是目标坐标系，-lco ENCODING=UTF-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 个字节的 LDID（Language Driver ID），0x57对应 ANSI/GBK，0x00表示未指定，未指定时软件会退回系统默认代码页，在中文 Windows 上是 GBK，在 Linux 容器里往往是 Latin-1，于是同一份数据在不同机器上表现不一样，这就是典型的"我这边好好的，那边全是问号"。

第三条最可靠：直接看字节。用十六进制看一眼第一个中文字段的值：

# 取 dbf 前 400 字节看结构，重点看字段值区域的字节 xxd -l 400 nyqh/nyqh.dbf | head -30

GBK 的中文是两个字节，首字节在0x81-0xFE，次字节在0x40-0xFE；UTF-8 的中文是三字节，形如0xE4-0xE9开头。看一眼就能区分。确认之后统一转成 UTF-8 并写出 cpg，后面所有工具都不再需要指定编码：

import geopandas as gpd # 读入时按实际编码解码 gdf = gpd.read_file("raw/nyqh.shp", encoding="gbk") # 写出时统一 UTF-8，geopandas 会自动生成 cpg gdf.to_file("work/nyqh_utf8.shp", encoding="utf-8")

这里有个细节：read_file的encoding只管 dbf 解码，不管.cpg；to_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(columns=rename_map) gdf.to_file("work/nyqh_clean.shp", encoding="utf-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, on="join_key", how="left", indicator=True) # 关键检查：有没有没匹配上的 print(merged["_merge"].value_counts()) print(merged.loc[merged["_merge"] == "left_only", "zone_name"].tolist())

how="left"加indicator=True是这套流程里最值钱的两个参数：前者保证区划要素不丢，后者直接告诉你哪几条没匹配上。打印出来的left_only列表就是需要人工补映射的名单，通常不会有几条，但每一条都可能导致某个区的统计值缺失。这一步不做，最后出图时某个区是空白，你还以为是数据里本来就没有。

5. 避坑与排查：五条血泪记录

这一章的五条都是我自己踩过或者帮别人排过的，症状都不明显，共同点是"看起来正常，数字不对"。每条按现象、原因、解决写，能对号入座就直接跳过排查过程。

5.1 五条具体记录

第一条：QGIS 里中文全是问号，Python 里正常。现象是同一份 shp，geopandas读出来是汉字，拖进 QGIS 显示????。原因是 dbf 里没有.cpg，QGIS 按自身设置的系统编码去解，而读数据的机器编码设置不同。解决是重写一遍带上 cpg：gdf.to_file(..., encoding="utf-8")，或者手工建一个同名.cpg文本文件写入UTF-8。注意别用编辑器保存成带 BOM 的 UTF-8，某些 GDAL 版本会把 BOM 当成编码名的一部分。

第二条：面积加起来比实际大了一截。现象是九个区面积合计明显偏大，比例在百分之十到几十之间。原因是用 Web 墨卡托或经纬度直接算的。解决是先重投影到 Albers 等积再算，并且确认crs的单位是米不是度。判断方法很直接：gdf.crs.axis_info[0].unit_name，打印出来是metre才算对。

第三条：某个区在做空间连接后消失了。现象是连接完要素数量比原来少。原因是属性表里的区划名称有不可见字符（尾部空格、全角空格、换行），或者用了how="inner"。解决是连接前统一strip并去掉全角空格\u3000，连接后用indicator=True检查left_only。全角空格这一条尤其隐蔽，肉眼完全看不出来，len()长度却对不上。

第四条：裁剪后边界出现细缝或重叠。现象是用区划做掩膜裁剪栅格，或者用区划做叠加分析，相邻区之间出现一两像素的缝。原因是相邻多边形共享边界但顶点不重合，这是不同批次数据拼接的常见结果。解决是裁剪前对整体做一次节点对齐，shapely里可以用unary_union后重新切分，或者在 GIS 里跑一次拓扑修复。对精度要求不高的场景，栅格统计时用all_touched=True能掩盖大部分缝，但矢量的叠加分析掩盖不了。

第五条：换台机器后脚本报"字段不存在"。现象是本地跑通的脚本放到服务器上，读 shp 时报列名不存在。原因是字段名被截断成了QYMC_1这种带后缀的形式，而不同 GDAL 版本对重名后缀的生成规则不完全一致。解决是读取后立刻按位置重命名成固定英文名，脚本里永远不出现原始中文字段名。这条是从根上断掉兼容性问题，比在每台机器上对齐 GDAL 版本省事。

5.2 一套最小自检脚本

每次拿到新的 shp（不管是这份区划还是别的），我都会先跑一遍这个脚本，几十秒的事，能提前发现八成问题：

import geopandas as gpd def health_check(path, encoding="utf-8"): gdf = gpd.read_file(path, encoding=encoding) 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_touched=False, # 只统计像元中心落在面内的 nodata=-9999, # 显式指定，别依赖 tif 自带的 geojson_out=False, ) # 结果是一个列表，顺序与矢量要素顺序严格对应 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", encoding="utf-8")

6.2 用面积占比反查区划边界是否可靠

统计做完别急着出图，先做一次反向核验：把每个区算出来的像元数乘以像元面积，和该区的几何面积比一比，比值应该接近 1。偏差超过一成，通常是三个原因之一——坐标系不一致（矢量和栅格不在同一投影下，zonal_stats不会报错但结果全错）、区划边界和栅格范围没对齐、或者有要素的几何是无效的。

还有一个更粗但有效的检查：九个区的面积占比加总必须是 100%。如果加出来是 103%，说明区划之间存在重叠，这在多来源拼接的数据里不罕见。处理方式是在做统计前先做一次unary_union再切分，或者接受重叠但在报告里注明口径。我一般会把这项检查写成脚本里的断言，跑不通就中止，不让自己带着错误结果往下走。

从那以后我每次拿到新的区划 shp，都强制走一遍"文件清点、坐标系确认、几何修复、字段重命名、面积占比核验"这五步，一步不跳。跳过哪一步，后面就得花几倍时间回头找。希望帮到你。

本文还有配套的精品资源，点击获取