资讯详情

资讯详情

长江经济带GIS矢量图全流程处理:坐标系检查、拼接转换与缓冲区分析

简介长江经济带GIS矢量图是一份面向GIS科研、城乡规划与区域经济分析人员的2015年长江经济带地级市专题底图数据。数据覆盖沿线全部地级行政区包含城市行政边界、地理位置等基础要素可在ArcMap中直接导入使用适用于区域差异比较、产业布局研究、生态环保评估和政策规划等场景。rar压缩包共7个文件以shp矢量主文件为核心配有dbf属性表、prj投影信息、sbn/sbx空间索引及shp.xml元数据整包大小约294KB文件结构清晰便于制图与空间分析。目前已有1660人学习下载。借助这份数据使用者可快速获取2015年长江经济带地级市城市格局并结合人口、交通、经济等属性字段深入探究区域协同、协调发展与生态保护问题为学术研究、政府决策与公众科普提供可靠的地理底图支持。1. 长江经济带GIS矢量图不是一张画而是一套能计算的空间数据底座第一次拿到的“长江经济带GIS矢量图”可能只是一个画了边界的图片也可能是几十个省级面片散落在一起属性表里连“经济带”三个字都没有。真正能用于区域分析的长江经济带矢量图应该是带坐标系、闭合边界、属性标记的空间数据集它能在GIS软件里叠加规划图层、做缓冲区分析、做数值统计。下面按我自己的落地流程展开从验收、清洗、拼合、格式转换到缓冲分析。适合GIS开发、数据分析、区域规划方向的一线从业者也适合刚接触空间数据的同学照着走一遍。如果你只想要一张静态示意图这篇未必合适如果你要对长江经济带做任何形式的量化分析后面每章都值得读。2. 验收一套长江经济带矢量图先查坐标系、属性表和几何完整性原始数据拿到手别急着写可视化脚本先花10分钟做验收能省掉后面一天的排障。下面这三个检查项是我每次处理区域边界数据的第一课。2.1 坐标系决定了这张图钉在地球上的位置坐标系是第一个检查项。长江经济带覆盖范围很广跨越多个投影分带。如果数据是地理坐标系单位是度适合底图展示如果数据是投影坐标系单位是米适合计算面积和距离。两者混用最常见的现象是“图层叠加位置不对”专业上叫套合偏差。你以为图形画错了实际上多半是坐标系没统一。用Python的geopandas读入矢量文件后直接看crs属性。这是最直观的检查方式。import geopandas as gpd gdf gpd.read_file(cjsjjj_raw.shp) print(坐标系:, gdf.crs) print(要素数量:, len(gdf)) print(范围:, gdf.total_bounds)这段代码完成三件事读取文件、打印坐标系、打印要素数量与数据范围。crs如果是None意味着数据文件缺少投影定义对于shapefile通常是同名的.prj文件没有一起拷贝。total_bounds返回四个数值分别是xmin、ymin、xmax、ymax即数据覆盖范围的左下角与右上角。如果这组数值在73到135之间且输入坐标是经纬度说明范围正常如果出现类似400000这样的七位数值说明数据是投影坐标需要先做一次to_crs转换再决定是否显示。常见误区是看到底图偏移就手动平移整个图层。在这个尺度下手动平移很容易“看着对上分析全偏”。我一般不会移动坐标系而是要求原始数据必须带清楚的crs定义。如果确实没有再结合已知控制点做空间校正。还有一个细节shapefile的投影信息并不总是正确。比如文件声称是WGS84实际是用GCJ02偏移过的坐标肉眼在沿江城市边界处能看出道路偏移这种只能通过比对已知地物坐标来验证不能轻信文件头。2.2 属性表里到底有没有“经济带”这个边界长江经济带不是自然界线而是一个经济区域概念。数据提供方通常会按省级行政区或流域面来组织要素属性表里需要有明确的字段说明一个要素是否属于经济带范围内。例如“区域名称”“所属流域”“经济带标记”这类字段。如果完全没有将来做统计分析时无法直接过滤出经济带内要素。用value_counts看字段的分布情况。print(gdf.columns.tolist()) if region_name in gdf.columns: print(gdf[region_name].value_counts(dropnaFalse)) else: print(没有找到 region_name 字段需要自行标记经济带范围)这段代码会先输出所有字段名再检查region_name字段是否存在。如果存在按分组统计要素数量如果不存在后续拼接时需要自己加统一标记。常见情况是一个省的面在原始数据里被切成很多小面片比如按县界发布省界特征不明显或者反过来整个经济带被提前合并成一个MultiPolygon属性表只有一行内部行政层级全部丢失。这些都要根据业务决定保留哪种粒度。我习惯把“经济带”理解为一个图层组而不是单个图层边界轮廓、行政点、水系、交通线各放一层属性里都带一个zone_code字段。这样后续任何图层都能快速过滤出经济带内的对象。如果数据来自公开基础要素通常会有省级行政区划代码字段这是后续连接统计年鉴、人口数据的关键。如果只有中文名称要先对齐编码表避免后面做空间连接时对不上。2.3 跑一遍完整性检查空几何、无效几何和边界范围一份看起来正常的矢量图内部可能藏有大量空几何、自相交多边形或孤立碎片。直接拿来做空间连接结果会被污染统计面积时出现负值生成缓冲区时出现异常自相交Web端渲染时出现白线。所以验收第三步是几何质量检查。gdf[is_valid] gdf.geometry.is_valid gdf[is_empty] gdf.geometry.is_empty gdf[geom_type] gdf.geometry.geom_type print(gdf.groupby(is_empty).size()) print(gdf[~gdf[is_valid]][geometry].head(10))代码中is_valid是shapely的拓扑有效性判断返回False表示该几何不符合拓扑规则。is_empty判断是否为空几何。geom_type返回Polygon、MultiPolygon、LineString等类型用来检查数据里是否混入了不该出现的线或点。如果is_valid为False的要素很多可以用第3章里的buffer(0)修复但如果修复后边界仍然异常就需要回到原始数据去筛查。结合total_bounds还可以判断数据是否漏了一块。比如总范围只覆盖中下游那说明打包时截断了上游区域。这种问题靠代码不一定能发现但结合业务常识很容易察觉。长江经济带覆盖11个省级行政区如果你的数据要素数量只有1说明做过一次合并边界轮廓有但省界信息丢如果要素数量是几十万个说明粒度太细直接做全量分析容易卡死。验收做完你已经能判断这包数据能不能用、用在哪一层。下一章处理最常见的坏数据。3. 从省级面片拼成长江经济带dissolve、buffer与属性标记如果拿到的是散碎的省级边界就需要自己拼成经济带轮廓。这一步最容易翻车有人用union合并后丢属性有人用clip裁出锯齿还有人不管拓扑直接在浏览器里画结果边界有缝。这里说说我的做法。3.1 用dissolve按字段合并而不是用union硬拼把多个面合并成一个整体在geopandas里有两种思路unary_union合并后返回shapely几何对象不带属性表dissolve则保留GeoDataFrame结构可以继续做属性操作。对于要做成图层的经济带首选dissolve。gdf[zone_code] CJJJ gdf_dissolved gdf.dissolve(byzone_code, aggfuncfirst) print(gdf_dissolved)代码先给每个要素打上同一个分组值然后按这个字段聚合。dissolve的aggfunc支持first、sum、mean等聚合方式。如果只是把多个省面拼成一个整体几何部分会融合属性部分用first取第一个要素的值即可。如果有面积字段想同时得到总面积可以写成aggfunc{area_field: sum}。为什么要用by指定字段而不是直接dissolve因为如果以后数据里还有“城市群”和“经济带”两类面一个zone_code可以运行一次生成多个类别。相比手动选择再merge这种方式更可复现。注意输入数据里不能有重复的行政代码否则可能把不同县也合并掉。还有一个容易误解的地方dissolve只是消除了面之间的内部边界不会自动修复拓扑错误。如果输入数据面片之间有交叉dissolve的结果仍然可能自相交。因此下一步要做几何清理。3.2 用buffer(0)处理细缝与自相交最常用的后悔药buffer(0)在GIS社区里被当作拓扑修复工具。它能消除零距离缓冲造成的不稳定让重叠部分合并、自相交处拆分同时不会显著移动原边界。但它的适用条件是有微小缝隙或纯拓扑冲突不是空间上相距很远的碎片。from shapely.ops import unary_union cleaned gdf_dissolved.geometry.buffer(0).unary_union gdf_cleaned gpd.GeoDataFrame(geometry[cleaned], crsgdf_dissolved.crs) gdf_cleaned gdf_cleaned.explode(index_partsFalse).reset_index(dropTrue)buffer(0)的作用是把每个多边形稍微膨大一丁点再缩回一丁点数值上半径虽然是0但shapely在计算时会做一次拓扑重构。如果两块区域相距数公里buffer(0)不会缝合。遇到这种大缝隙要先检查原始数据是否漏了中间要素而不是盲目补洞。explode是把合并后的MultiPolygon拆成多个Polygon要素保留每个岛屿或飞地作为独立记录。这一步在实际业务中经常需要因为后续要对每个面片单独命名。如果只关心整体轮廓可以省略explode。经过buffer(0)后几何顶点数量会变多直接保存会增大文件体积建议在下一章做simplify。为什么叫后悔药因为很多人在合并前没意识到边界有缝直到出图时发现统一填充的颜色里透出底图才回头处理。加一个buffer(0)能解决80%的案例但不要误以为它能修复所有拓扑损坏。比如自相交多边形buffer(0)之后可能变成两个或多个新多边形这会改变统计口径使用时要留意。3.3 补上区域标识字段让数据可被过滤、可被连接合并后的数据如果只有一个面且属性表空荡荡用起来还是不方便。你需要给它补充业务字段比如区域名称、代码、版本号。有的项目还需要保留合并前覆盖的省级行政区数量作为元数据记录。gdf_cleaned[zone_name] 长江经济带 gdf_cleaned[zone_code] CJJJ gdf_cleaned[province_count] 11 gdf_cleaned[update_year] 2024 gdf_cleaned.to_file(cjsjjj_merged.shp, encodingutf-8)多出来的几个字段中zone_code建议全英文否则在旧版GIS软件里容易因字段名问题报错。province_count是手工写入的说明合并前按业务口径应该覆盖11个省级行政区。如果后续发现实际是10个要通过这个字段回溯过程而不是只看几何。保存shapefile时geopandas默认写UTF-8但很多桌面GIS软件默认按GBK读所以最好在读取端显式设置encoding或者直接保存GeoJSON避免编码问题。我通常保存两个版本一个原始投影的shapefile给团队内部做分析一个WGS84的GeoJSON给Web端显示。这一步做完经济带矢量图的“底座”就干净了。4. 数据格式与转换Shapefile、GeoJSON、TopoJSON怎么选怎么转一份矢量图最终给谁用决定了它的存储格式。本地跑分析、Web渲染、第三方交接三种场景的选型不一样。这部分从实际项目视角讲清楚取舍给出可直接复制的转换代码。4.1 三种常见格式各在什么场景下用格式适用场景优点注意点Shapefile传统桌面GIS、批量制图、历史项目兼容性极强基本所有GIS软件都能打开由多个文件组成字段名限10个字符属性编码容易乱GeoJSONWeb地图、前后端开发、开源GIS单文件、JSON结构直观浏览器可解析文件偏大没有内置拓扑压缩坐标必须是经纬度TopoJSON大范围区域概览、Web端轻量展示拓扑压缩明显边界共享只存一份需要前处理部分软件读取需要额外插件按照我的个人习惯长江经济带这种大范围边界给Web端优先用TopoJSON给数据分析同事用shapefile给需要跨部门协作的用GeoJSON。三者之间用geopandas一行代码即可转换。4.2 用geopandas完成转换驱动和编码是最大的两个坑转换代码看起来简单细节都在参数里。以转GeoJSON为例gdf_4326 gpd.read_file(cjsjjj_merged.shp) gdf_4326 gdf_4326.to_crs(epsg4326) gdf_4326.to_file(cjsjjj_merged.geojson, driverGeoJSON, encodingutf-8)read_file读入后先转到EPSG:4326。GeoJSON标准要求使用WGS84经纬度如果把带投影坐标的数据直接写成GeoJSON浏览器上会错位。driver参数显式指定GeoJSON虽然从扩展名也能推断但写出来能避免某些环境下默认驱动不一致。encoding指定UTF-8防止中文字段名或属性值乱码。Shapefile的“多文件”问题也值得注意。shapefile不是单文件.shp、.shx、.dbf、.prj一起才算完整。很多人只拷贝.shp结果打开报错。我在交接时都会把整个目录打包成zip并在zip里放一个README说明坐标系和编码。这里没有额外命令但很值得养成习惯。4.3 简化边界让Web端加载不再卡顿长江经济带边界如果来自高精度省界节点数可能达到几十万。原样转成GeoJSON一张图几MB浏览器渲染会卡。简化是必须的。shapely的simplify基于道格拉斯-普克算法减少顶点gdf_simple gdf_4326.copy() gdf_simple[geometry] gdf_simple.geometry.simplify( tolerance0.01, preserve_topologyTrue ) gdf_simple.to_file(cjsjjj_simple.geojson, driverGeoJSON)tolerance的单位与当前坐标系一致。这里是经纬度0.01度大约对应1公里。如果想保留更多细节可以设为0.005如果只是做全国概览0.05即可。不要为了精简把tolerance设到0.5否则边界会变成锯齿严重失真。preserve_topologyTrue能防止简化后产生自相交面但也会让算法更保守速度与文件大小的平衡需要试算几轮。保存后对比文件大小并抽查经济带轮廓是否平滑。我会在GIS软件里把简化前后两个图层重叠对比肉眼差别不大才算合格。tolerance选0.01只是起点具体值取决于原始边界精度和用途。5. 落地避坑指南长江经济带矢量图使用时的5个典型翻车现场数据整理完不代表项目能顺利跑通。下面这5个问题是我在实际项目中反复见过的按出现频率排列。每条都按现象、原因、解决来写。5.1 底图偏移几十米坐标系基准被悄悄替换现象经济带边界叠加到在线卫星底图上沿江边界明显偏到河一侧东北方向位移几十到几百米。原因原始数据用的是GCJ02或BD09坐标系但文件头标注为WGS84。这两种坐标系与标准经纬度存在非线性偏移且不同区域偏移量不同。解决先检查图层属性里的坐标信息如果确认坐标偏移用geopandas转换到EPSG:4326并不能自动纠偏。更可靠的做法是找到原始数据说明确认坐标基准。没有任何说明时只能基于已知控制点评估偏移量并明确告诉最终用户该数据只能用于相对位置分析不能和标准影像严格套合。关键是不要手动平移硬凑因为不同区域偏移量不是常数。5.2 合并后的细尖角dissolve没有修复拓扑冲突现象dissolve出的经济带轮廓上有一堆细长三角尖刺放大看是很多极窄的条带。原因原始省际边界之间有很多微小重叠dissolve合并了属性但几何上这些重叠区的边界没有真正对齐。解决在dissolve后执行buffer(0)然后再做unary_union。这一步能合并重叠边段消除尖角残片。如果仍然有小刺可以再做一次simplify它会将摆动幅度小于tolerance的凸起抹平。gdf_fixed gpd.GeoDataFrame( geometry[gdf_dissolved.geometry.buffer(0).unary_union], crsgdf_dissolved.crs )这里buffer(0)后直接构造新的GeoDataFramegeometry列表里只有一个unary_union结果。注意要先buffer(0)再simplify顺序不能反否则简化反而会把未处理的接缝重新放大出来。5.3 属性表中文乱码Shapefile的dbf编码黑匣子现象用GIS软件打开属性表中文全变成问号或乱码字段名也显示异常。原因Shapefile的dbf编码需要读取端猜测常见有GBK和UTF-8两种。写入时用UTF-8读取端按GBK解就会乱。解决统一在写文件时指定encodingutf-8读文件时显式指定encodingutf-8或gbk。如果数据已经乱了用Python读出来重新写一遍。最稳妥的方案是字段名全部用英文属性值里的中文尽量转成拼音或代码中文说明单独放到说明文档中。这不是能力问题是省心。5.4 Web端加载GeoJSON白屏文件太大顶点太多现象在浏览器里拖拽缩放地图渲染明显掉帧移动端直接白屏。原因原始边界顶点数量巨大GeoJSON没有拓扑压缩浏览器需要解析长数组并重复计算边界关系。解决用simplify做减法再考虑转为TopoJSON。如果要保留多级行政边界可以分层加载底图用低精度轮廓放大到一定级别再请求高精度边界。经验值是浏览器单图层控制在2MB以内如果目标是移动端最好低于500KB。简化tolerance要根据缩放设置展示经济带整体轮廓时0.02度已经够用。5.5 空间统计少算两个城市点面关系与线面关系混用现象用城市中心点做空间连接统计经济带内的城市结果明显漏了边界附近的区域。原因中心点不在多边形内但城市行政区边界与经济带有大量重叠。用点面相交判定会把这类城市整体排除。解决改用面面相交判断。先让经济带面与城市行政边界求交得到交集面积再定义纳入规则。例如“交集面积占比超过50%”或“只要相交就纳入”。用geopandas overlay可以拿到每个城市的交集面再groupby汇总。这比点面判断更符合“流域范围覆盖城市”的业务含义。6. 进阶用缓冲区分析快速定位“沿江50km覆盖了哪些区域”到了这一步你已经有了干净的长江经济带矢量图。最后一个具体技巧以长江干流线为基准生成50km缓冲区统计这个缓冲区内涉及的县级行政区。这是我们经常用到的“沿江可达域”分析基础。先准备两个数据源长江干流线图层和县级行政区面图层。第一步把坐标系统一到米制投影我这里选用EPSG:32650适合中下游局部区域如果项目跨整个经济带可以按流域分段选择不同投影。第二步调buffer。第三步用overlay求交按县汇总面积。river gpd.read_file(yy_river_line.shp).to_crs(epsg32650) counties gpd.read_file(counties.shp).to_crs(epsg32650) buf gpd.GeoDataFrame( geometryriver.geometry.buffer(50000), crsriver.crs ) buf[buffer_type] 50km hit gpd.overlay(counties, buf, howintersection) hit[area_km2] hit.geometry.area / 1e6 result hit.groupby(county_name)[area_km2].sum().round(2) print(result.sort_values(ascendingFalse).head(15))buffer(50000)生成的缓冲区半径是5万米即50公里。overlay howintersection返回两个图层相交后的面片也就是县域落在缓冲区内的那部分。area_km2是对相交面片面积换算成平方千米groupby后得到每个县被覆盖的面积。这个方法可以无缝推广到“离高铁站3km”“离产业园区5km”等可达性场景。我第一次做这个分析时忘记先统一投影直接在经纬度上buffer结果缓冲区在东西方向被拉伸统计面积偏差很大。后来才意识到投影坐标系下的距离才可信。如果分析涉及全国尺度可以用等积投影如果只是中下游局部区域UTM分带投影足够。这个习惯我一直保持到现在。希望帮到你。本文还有配套的精品资源点击获取
觉得有用,分享给同行:

为您的企业打造数字门面

稳重轻奢商务风格,端正雅致视觉,长效耐看不易过时。

立即咨询 →