简介：《中国九大农业区划shp》数据集面向农业地理、国土规划与GIS分析方向的研究者、教师及学生，它把全国划分为黄淮海平原、东北平原、长江中下游、北方干旱半干旱、云贵高原、青藏高原、四川盆地、华南、黄土高原九个区，为种植结构对比、土地利用效率评估与区域农业政策研究提供统一的空间边界底图。压缩包共11个文件，约17.21MB，由shp主文件承载区划面要素，配套dbf属性表、prj投影信息与shx索引，另有sbx、sbn等空间索引及xml元数据，可在ArcGIS、QGIS等平台直接加载、裁剪与叠加分析。目前已有4340人学习下载。利用这些数据，读者可完成分区统计、作物产量比较、气候与政策影响的叠加分析，乃至灾害风险评估与农业发展情景模拟，省去自行描边与配准的工作量，适合作为农业区划类课程作业、论文制图与项目前期的底图使用。

1. 九个面、五个文件：这类农业区划 shp 到底能省掉多少活

做农业遥感、乡村规划或者涉农数据看板的人，早晚会撞上同一个需求：手上有一堆点位、一张 NDVI 栅格、一份县级统计表，老板要你「按大区归并成九行」。这时候最省事的路径不是自己拿省级边界去 dissolve，而是直接找一份现成的中国九大农业区划 shp。压缩包解开就是一套面数据：九个多边形、一套区划代码与名称字段、配套的投影与属性文件。它真正解决的是「口径统一」——九个大区的边界是固定的，你算出来的面积、单产、覆盖率才有可比性，不用每次跟别人对数据时先吵半小时「你这个华北是怎么划的」。适合三类人：做分区统计的分析岗、做专题图制图的规划岗、以及要给下游系统塞一张稳定底图的开发岗。不适合想要省、市、县精细边界的人，那是另一套数据。

2. 拆包先体检：shp 五件套、字段编码与 .prj 的三项检查

拿到压缩包，第一步永远不是打开软件画图，而是把文件摊开看清楚。Shapefile 不是「一个文件」，是一组同名前缀、不同后缀的文件捆在一起工作，少一个就可能读不全或者读出来是错的。

2.1 五件套各自管什么，少哪个会翻车

我一般会先把它们放成一个独立目录，前缀统一改成英文短名（比如nyqh_9.shp），因为有些老版本工具对中文路径和中文文件名处理得很糟，路径里带空格也会让命令行脚本莫名失败。改名时务必五个文件一起改，只改.shp会导致找不到同名 dbf——这是新手最常见的第一次翻车。

提示：如果压缩包里还带了.sbn、.sbx、.xml、.qix，不用管，那是空间索引和元数据，删掉不影响读取。

2.2 用 GeoPandas 做第一次体检

不装 QGIS 也能体检，Python 里几行就够。下面这段是我每次拿到新 shp 都会跑一遍的固定动作：

import fiona import geopandas as gpd path = "nyqh_9/nyqh_9.shp" # 1) 用 fiona 先看元信息，不加载几何，快且不会因为几何错误中断 with fiona.open(path) as src: print("driver:", src.driver) # 正常是 ESRI Shapefile print("crs :", src.crs) # 空 = 没有 .prj，后面必须自己补 print("schema:", src.schema) # 字段名/类型，注意 dbf 字段名上限 10 字符 print("count :", len(src)) # 九大区正常应等于 9 # 2) 再真正读进来。中文乱码时把 encoding 换成 gbk 或 gb18030 gdf = gpd.read_file(path, encoding="utf-8") print(gdf.dtypes) print("bounds:", gdf.total_bounds) # 一眼判断坐标系对不对 print(gdf.geom_type.value_counts()) # 期望全是 Polygon 或 MultiPolygon print("invalid:", (~gdf.is_valid).sum())

这段代码的价值集中在三个输出上。crs为空说明.prj丢了，后面所有叠加操作都得先手工指定坐标系，否则 GeoPandas 会把它当无名 CRS 处理，一to_crs就报错。total_bounds是最便宜的坐标系验证手段：如果四个值都在 70～140 和 3～55 这个区间，说明单位是度，是地理坐标系；如果出现六七位数甚至八位数，说明单位是米，是投影坐标系——两者混用就是后面所有错位的根源。geom_type那一行要留意有没有MultiPolygon，九个大区里带岛屿或飞地的单元往往是多部件，这会影响你后面按要素计数和分区统计的写法。

2.3 字段探查：找出区划代码和名称列

不同来源的区划数据字段命名五花八门，QHDM、QH_CODE、AREA_ID、NAME、区划名称都见过。与其猜，不如把所有文本型字段打出来看：

# 属性表摸底：文本型字段的取值数、样例、长度 for col in gdf.columns: if col == "geometry": continue if gdf[col].dtype == object: vals = gdf[col].astype(str) print(f"{col}: nunique={vals.nunique()} len_max={vals.str.len().max()}") print(" sample:", vals.head(3).tolist())

判断标准很朴素：唯一值正好是 9、长度在 2～10 之间、内容是中文的那一列，就是区划名称；唯一值是 9 且形如01…09或1…9的，就是区划代码。找到之后建议立刻做两件事——把代码列转成两位字符串（zfill(2)），避免1和01在后续 join 时对不上；把名称列重命名成统一字段（我习惯用zone_name），这样换一份区划数据时下游脚本不用改。

顺便看一眼 dbf 字段名有没有被截断成QHMC_1、DESCRIP_2这种。Shapefile 的 dbf 字段名硬上限是 10 个字符，超过就被截，两个长字段前 10 位相同时还会自动加后缀，这是 Shapefile 格式的天生缺陷，不是数据做错了。

3. 面积口径定生死：从经纬度到亩，投影怎么选、参数怎么设

九大区划数据最常见的用途就是算面积、算占比。而面积这件事，方法选错会让结果差出一个量级，且错得很隐蔽——图看着是对的，数字是错的。

3.1 为什么不能在经纬度坐标里直接算面积

gdf.geometry.area在 EPSG:4326 下返回的是「平方度」，这是个没有物理意义的量。更麻烦的是，一平方度在不同纬度代表的实际面积完全不同：低纬度一平方度接近 1.2 万平方公里，高纬度只有它的一半左右。如果你直接拿平方度去排序或者算占比，等于默认全国各地的「一度」一样大，北方的大区会被系统性低估。这不是精度问题，是逻辑错误。

正确做法是先把数据转换到一个等积投影，再算面积。等积投影保证投影前后面积比例不变，代价是形状会轻微变形——用来算面积、做统计没问题，用来做印刷级制图就得换回等角投影或者直接上经纬度加标注。

3.2 Albers 等积投影的参数怎么填

国内做全国尺度面状统计，惯用 Albers 圆锥等积投影。下面这套参数是我在多个项目里反复用的：

from pyproj import CRS # 双标准纬线 Albers：适合东西跨度大、南北跨度大的全国尺度 albers = CRS.from_proj4( "+proj=aea " "+lat_1=25 +lat_2=47 " # 双标准纬线，取数据南北跨度的 1/6 与 5/6 附近 "+lat_0=0 " # 原点纬度，全国尺度取 0 即可 "+lon_0=105 " # 中央经线，取数据东西范围的中点 "+x_0=0 +y_0=0 " "+ellps=GRS80 " # 椭球，配 CGCS2000；若原始数据是 WGS84 也用 GRS80 差别可忽略 "+units=m +no_defs" ) gdf_m = gdf.to_crs(albers) gdf_m["area_km2"] = gdf_m.geometry.area / 1e6 gdf_m["area_mu"] = gdf_m["area_km2"] * 1500.0 # 1 km² = 1500 亩 print(gdf_m[["zone_name", "area_km2", "area_mu"]])

参数里真正需要你根据数据调整的是lat_1、lat_2和lon_0。双标准纬线的经验取法是：设数据纬度范围为 φmin 到 φmax，则lat_1 = φmin + (φmax-φmin)/6、lat_2 = φmax - (φmax-φmin)/6。这两条纬线上投影无变形，中间和外侧的变形被压到最小。lon_0取经度范围中点，能让中央经线两侧的变形对称。如果这两个参数随便填，面积误差不会大到离谱（通常千分之几），但当你需要跟官方公布的播种面积对账时，这点误差就足够让你怀疑人生。

3.3 三套单位互转与一次交叉校验

农业口用亩，规划口用公顷，遥感口用平方公里，一份数据要伺候三种说法。换算是死数：1 平方公里 = 100 公顷 = 1500 亩。真正值得花时间的是校验，我固定做两步：

第一步，把九个区的面积加总，看是否落在全国陆地面积的合理区间。九个大区理论上应该基本铺满，加总明显偏小说明区划之间有缝或者被裁掉了海疆部分，明显偏大说明区与区之间有重叠。

total_km2 = gdf_m["area_km2"].sum() print(f"九区合计: {total_km2/1e4:.1f} 万 km²") print(gdf_m["area_km2"].describe())

第二步，换一套投影再算一遍，看两次结果的相对偏差。偏差在 0.5% 以内说明数据本身没问题；如果某个区单独偏差超过 2%，八成是这个区的几何有自相交或者多部件没合并，第 6 章会讲怎么查。

注意：网上流传的「1 平方公里 = 1500 亩」是精确值（1 亩 = 666.67 平方米），别用 1000 或者 1495 这类近似值，几十个大区汇总下来误差会被放大。

4. 分区叠加实战：把点位、栅格和县级面挂到九个大区上

数据体检过关、投影定好，接下来才是真正的业务动作。九大区划 shp 本质上是一张「归类底图」，所有活儿都是拿别的数据去贴它。三种贴法，写法完全不同。

4.1 点位落区：sjoin 的三个参数决定结果对不对

把采样点、监测站、企业点位分配到九个大区，用空间连接。这段代码看着短，坑全在参数上：

import geopandas as gpd zones = gpd.read_file("nyqh_9/nyqh_9.shp")[["zone_code", "zone_name", "geometry"]] pts = gpd.read_file("sites.shp") # 关键 1：两边坐标系必须一致，不是「看起来差不多」而是必须完全相同 pts = pts.to_crs(zones.crs) # 关键 2：predicate 用 within 表示「点落在面内部」， # 落在边界线上的点用 within 会匹配失败，需要 intersects 兜底 joined = gpd.sjoin(pts, zones, how="left", predicate="within") # 关键 3：left join 会留下未匹配的点，index_right 为空就是没落进任何大区 miss = joined[joined["index_right"].isna()] print("未落区点数:", len(miss))

how="left"是我强烈建议的默认选择。用inner虽然结果干净，但会把落在边界缝隙里、或者根本在区划范围外的点悄悄丢掉，你看到的统计结果永远是「全部落区」，问题被掩盖了。用 left 保留下来，先看未匹配数量：如果是 0，皆大欢喜；如果是几十个，先查坐标有没有写错（比如把经纬度顺序搞反）；如果是几千个，那多半是点位数据本身覆盖了区划之外的区域，属于正常，按需丢弃并记录。

predicate的选择要看业务定义。「点在大区内」严格说是within，但海岸线、大区交界处的点经常因为浮点精度落在边界上，within判定为 False。这类点如果业务上允许归属任一侧，就改用intersects，然后用drop_duplicates按点 ID 去重——因为边界上的点可能同时命中两个大区。

4.2 栅格分区统计：all_touched 是个统计口径问题

NDVI、降水、积温这类栅格，要按大区求均值，用rasterstats：

from rasterstats import zonal_stats stats = zonal_stats( zones, # 面数据，坐标系要与栅格一致 "ndvi_2023.tif", stats=["mean", "count", "std"], # count 是参与计算的像元数，用来判断可信度 all_touched=False, # False：只算像元中心落在区内的 nodata=-9999, # 必须显式指定，否则 NoData 会被当数值参与平均 categorical=False, ) zones["ndvi_mean"] = [s["mean"] for s in stats] zones["px_count"] = [s["count"] for s in stats]

all_touched不是精度开关，是口径开关。设成False时只统计像元中心落在区内的像元，结果是保守估计；设成True时只要像元与区有任意重叠就算进去，边缘区会被系统性放大。大区面积大、形状规整时两者差别很小；但对于形状狭长或者带大量碎岛的区域，True会让有效像元数多出一两成，均值被边缘混合像元拉偏。我的习惯是统一用False，然后在文档里写清楚口径，避免半年后有人拿另一套结果来对数时扯皮。

nodata一定要显式写。很多栅格产品用的是-32768或者nan，不声明的话这些值会直接进平均，结果会变成一堆负数或者nan，而你还以为是数据质量问题。

4.3 面与面叠加：字段爆炸和面积重算

要算「每个大区内各县的耕地面积」，就得做面面相交。用overlay：

county = gpd.read_file("county.gpkg")[["adcode", "adname", "geometry"]] county = county.to_crs(zones.crs) inter = gpd.overlay(zones, county, how="intersection", keep_geom_type=True) # 叠加后几何变了，面积必须在等积投影下重新算，不能沿用原始面积字段 inter_m = inter.to_crs(albers) inter_m["area_km2"] = inter_m.geometry.area / 1e6 # 按大区汇总 summary = inter_m.groupby("zone_name", as_index=False)["area_km2"].sum()

这里有两个必须知道的副作用。第一，overlay会给两边的同名字段自动加后缀（_1、_2），如果你的两个图层都有name字段，结果里就会出现name_1、name_2，后续代码引用错列是高频事故。稳妥做法是叠加前先把字段裁到最少、并统一命名。第二，叠加后的碎多边形面积必须重算——原表里带的Shape_Area字段是叠加前的，直接 sum 会得到「比实际大」的结果，因为一个县跨两个大区时它的面积会被重复计入两次。

keep_geom_type=True也值得说一句。面面相交理论上还是面，但因为边界处的浮点误差，可能会产出极小的线或点几何，设成 True 会把它们丢掉；如果你的业务不能容忍丢要素，就设 False 然后自己检查geom_type分布。

5. 避坑与排查：乱码、无 prj、多部件和缝隙的五个血泪现场

这一章全是真金白银换来的。九大区划 shp 这类「流传已久、来源不明」的数据包，问题往往不是「能不能读」，而是「读出来是不是你以为的那个东西」。

5.1 属性表中文全是问号或乱码

现象：read_file之后区划名称列显示成????或者一堆锟斤拷，但代码列正常。

原因：dbf 是上世纪格式，没有强制的编码声明，全靠.cpg文件或者读取端猜。国内数据大多用 GBK/GB18030 存中文，而 GeoPandas 默认按 UTF-8 读，一读就崩。有些包干脆把.cpg删了，读取端只能瞎猜。

解决：先试encoding="gb18030"（GBK 的超集，兼容性最好），不行再试gbk。如果.cpg文件存在但内容是UTF-8而实际数据是 GBK，把.cpg内容改成GBK再读。实在搞不定就用 QGIS 打开：它会在识别失败时弹出编码选择框，选对了再另存为 GPKG，以后就再也不会有编码问题——GPKG 强制 UTF-8。

5.2 没有 .prj，叠加后整个图飞到地图外面

现象：把区划层和底图一叠，区划缩成一个小点或者跑到屏幕外，坐标数值看着像 100 多。

原因：.prj丢失后，GeoPandas 读进来的 CRS 是None。此时如果对另一个有 CRS 的图层做to_crs反投影，或者直接用某个库自动猜测，坐标系会被当成「米」或者被当成「其他地理坐标系」，结果就是你熟悉的玄学错位。

解决：先看total_bounds。数值在 ±180 和 ±90 内，就是经纬度，用gdf = gdf.set_crs("EPSG:4326")声明一下。注意是set_crs不是to_crs——前者只贴标签不改坐标，后者会真的做数学变换，用错了数据就被二次扭曲，这是最容易被忽略的一步。声明完再to_crs到目标投影，一切正常。

5.3 一个区对应多个面，按要素计数永远是错的

现象：属性表里只有 9 行，但len(gdf.explode())出来是十几行；某些区做count统计时数量翻倍。

原因：带岛屿、飞地或跨河段的大区，几何类型是MultiPolygon，一个要素里塞了多个独立面片。多数情况下这没问题，但一旦做空间连接或者栅格统计，多部件会被按面片分别处理。

解决：想按面片算就显式explode(index_parts=False)，想按大区算就在聚合前dissolve(by="zone_code")合并回单一面。我一般会在体检阶段就判断：gdf.geom_type.value_counts()里只要有MultiPolygon，后面的所有分组逻辑都按zone_code走，绝不按行号走。

5.4 区与区之间有缝隙和重叠，加总对不上

现象：九个区面积加起来跟预期差了几万平方公里；或者overlay之后多出一堆面积接近 0 的碎屑多边形。

原因：这类区划数据往往是多来源拼接或者从别处矢量化来的，相邻区边界不共线，A 的右边界和 B 的左边界有几十米的错位，形成一条细缝；也有数字化时手抖造成的重叠。1:400 万以下的小比例尺数据尤其明显。

解决：先量化再修。用unary_union把所有几何合起来，看合并后的面积与各自面积之和的差值，差值就是净重叠量。修的时候优先用shapely.make_valid()处理自相交，缝隙则用 0.001 度级别的buffer外扩后buffer回缩做一次「闭合」——注意这会轻微改变边界，做面积统计前要记录修正量。如果数据要长期用，值得花半天在 QGIS 里对相邻边界做一次拓扑对齐（snap），一劳永逸。

5.5 同名字段在 join 之后类型变了

现象：zone_code原本是字符串01，跟另一张表 join 完变成整数1，匹配率骤降。

原因：pandas 读 dbf 时会把看起来像数字的字符串列推断成数值型，而另一张表里存的是字符串，1 != "01"，静默匹配失败。

解决：所有用作连接键的区划代码，读入后立刻统一astype(str).str.zfill(2)，两边都做，不要只做一边。

提示：以上五条里，第 2 条和第 4 条最容易造成「结果看起来正常但其实全错」，因为它们不报错。任何一次分区统计出结果后，先看面积加总、再看未匹配数，这两个数不对就别往下走。

6. 进阶用法：拓扑校验与一份能反复用的九大区底图

数据修完之后，真正让它产生长期价值的一步是「固化」。我现在的习惯是：任何一份区划数据，只要在一个项目里用过一次，就一定沉淀成一份带校验的底图文件，后面所有项目直接引用。下面这套流程跑完大约十分钟，能省掉未来几十次重复排查。

先做拓扑体检，三条硬指标：

import geopandas as gpd from shapely.validation import make_valid z = gpd.read_file("nyqh_9/nyqh_9.shp").to_crs(albers) # 1) 自相交检查：无效几何会让面积和叠加全部失真 z["geom_ok"] = z.is_valid if not z["geom_ok"].all(): z["geometry"] = z.geometry.apply(make_valid) # 2) 重叠检查：两两求交，非空且面积超过阈值的才算真重叠 overlaps = [] for i in range(len(z)): for j in range(i + 1, len(z)): inter = z.geometry.iloc[i].intersection(z.geometry.iloc[j]) if not inter.is_empty and inter.area > 1e6: # 1 km² 以上才算问题 overlaps.append((z.index[i], z.index[j], inter.area / 1e6)) print("重叠对数:", len(overlaps), overlaps[:5]) # 3) 唯一性检查：区划代码必须唯一且非空 assert z["zone_code"].is_unique and z["zone_code"].notna().all()

重叠面积阈值我设成 1 平方公里，理由是小于这个量级基本是浮点误差，修它反而引入新问题；超过就说明边界确实有问题，得回去查源数据。这个阈值不是标准，是按全国尺度数据常见的数字化精度调的，换成县级细数据要往下调到 0.01。

校验通过后再做两件事：一是用simplify(tolerance=200, preserve_topology=True)生成一份轻量版，专供 Web 端渲染——保留拓扑意味着简化后相邻区依然是共边的，不会裂开；二是导出成两种格式，GPKG 存原始精度给分析用，GeoJSON（EPSG:4326）给前端用。

# 一次导出两种用途的底图，GDAL 命令行比写脚本更快 ogr2ogr -f GPKG nyqh_9_full.gpkg nyqh_9.shp -nln zones -t_srs EPSG:4547 ogr2ogr -f GeoJSON nyqh_9_web.geojson nyqh_9.shp -simplify 0.002 -t_srs EPSG:4326

-t_srs EPSG:4547是 CGCS2000 的三度带投影，如果数据本来就是这套基准，转过去几乎无损；-simplify 0.002是度单位的容差，约等于 200 米，Web 端在全国尺度下看不出来，文件体积能小一个量级。要落到数据库里做多用户共享，就再补一句 PostGIS 建表，把zone_code建成主键，几何列建 GiST 索引：

CREATE TABLE zones ( zone_code varchar(4) PRIMARY KEY, zone_name varchar(32) NOT NULL, area_km2 numeric(12,2), geom geometry(MultiPolygon, 4547) ); CREATE INDEX idx_zones_geom ON zones USING GIST (geom);

最后说个验证方法，比任何断言都直观：把九区底图和一个已知的省级边界层叠在一起，用「大区面积之和 ÷ 省级面积之和」算个比值。正常情况下这个比值应该非常接近 1（差额来自海疆和边界精度），如果落在 0.9 以下说明区划有大量空缺，超过 1.02 说明重叠严重。这个数字我每次交付前都会算一遍，当成最后一道后悔药。

这套流程不是天生就有的。早些年我直接拿一份来路不明的区划 shp 做统计，图出得漂漂亮亮，报告交上去两周后被人发现九区面积加总比实际少了将近三十万平方公里，原因是相邻区之间有系统性缝隙，当时没人做加总校验。从那以后我每次拿到任何面状区划数据，都强制先跑一遍「面积加总 + 重叠检查 + 未匹配点计数」这三件事，跑完才允许进业务流程。九大农业区划这份数据本身质量算不错的，但版本多、来源杂，落到你手上是哪一版只有天知道，体检这十分钟省不得。希望帮到你。

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