YAOTU INSIGHTS

从shapefile到TIF:羌塘高原区GIS底图制作与海拔提取全流程

从shapefile到TIF:羌塘高原区GIS底图制作与海拔提取全流程
简介这份面向地理信息与科研制图的数据包围绕羌塘高原区地理位置与海拔高度数据组织适合需要开展高原区地图绘制、空间分析或论文配图的GIS用户。压缩包共14个文件涵盖可编辑mxd工程、完整shapefile矢量组、TIF标准成图并附世界国家级与中国省级行政区划shape素材其中prj、dbf等配套文件保证了投影坐标和属性信息完整整体包体约186MB方便不同熟练程度的用户直接调用或二次编辑。已有80人浏览学习内置Word预览图可在下载前确认成图样式解压全部压缩文件夹后可在ArcGIS等平台中一键打开mxd使用。除基础地理数据外还兼顾矢量、栅格与工程文件三种形态mxd已预设常用图层符号和比例尺TIF成图可直接放入文档或报告结合行政区划底图还能灵活裁剪研究区范围为高原区专题制图和成果展示节省整理时间。1. 羌塘高原区底图三件套先解决坐标系再谈制图接到羌塘高原区的底图需求时大多数人第一反应是去下载DEM和边界数据结果往往卡在坐标系不一致、边界精度差、成图风格不统一三件事上。羌塘高原区这份资源把研究区边界、海拔信息、可编辑工程和标准成图打包在了一起mxd文件保留了ArcGIS里调好的图层组织标准shape文件提供可直接读取的位置与属性结构TIF成图则是带坐标的底图能直接放进论文或汇报。适合要做区域底图、野外路线图、海拔分级图的GIS从业者也适合需要拿现成范围去做裁剪、统计和批处理的Python用户。2. Shapefile八件套从shp到prj解析羌塘高原区数据底座2.1 一份shapefile为什么是八件套不少人把.shp当成单个文件拖进ArcGIS没问题换到QGIS、GeoPandas或者自己写读取脚本时出了问题才发现少了文件。shapefile是多文件地理数据格式主文件与附属文件必须在同一目录且文件名前缀一致任意缺一个都会影响要素读取、属性关联或坐标解释。“羌塘高原区shp文件”目录下是一套很典型的组合完整构成如下文件扩展名存储内容缺失时的影响.shp几何要素多边形边界数据完全无法打开.shx几何与属性的对应索引部分软件读不出形状.dbf属性表名称、面积、海拔等字段图形在但属性丢失.prj坐标系WKT定义投影信息丢失叠图错位.cpg属性编码声明中文属性乱码.sbn / .sbxArcGIS空间索引读取无影响编辑变慢.shp.xml元数据影响检索与识别表格里的.sbn和.sbx不算必需文件这份数据里带了说明它经历过ArcGIS的编辑和空间索引重建。验收shp时不能只看能不能打开还要确认.prj存在、.dbf里字段数量符合预期。我一般会把整包文件复制到工程目录保持八个文件在同一层且不改名否则后续用Python读dbf时会出现字段对不上或字符编码报错的问题。2.2 prj坐标系定义先读再动手打开羌塘高原区.prj里面是一段WKT形式的坐标系描述。WKT虽然冗长却是所有GIS软件判断两个数据能否在同一空间对齐的关键。文件里通常会出现PROJCS、GEOGCS、DATUM、PROJECTION、PARAMETER这些关键字分别对应投影坐标系、地理坐标系、基准面、投影方式和投影参数。用记事本或者命令行直接查看内容cat 羌塘高原区.prj查看时重点确认三件事椭球体是WGS84还是CGCS2000投影是UTM还是高斯-克吕格带号是多少。这组参数决定了后续叠加全国高程TIF时能否完全对齐。如果prj定义的是经纬度而高程TIF里是投影坐标叠加前必须先做投影转换大多数高程TIF采用WGS84或CGCS2000经纬网所以我通常先用gdalwarp -t_srs EPSG:4326把TIF统一到shp坐标系再继续裁剪。prj文件本身不要随意编辑它一旦被改动所有依赖该坐标系的显示和计算都会跟着偏移。2.3 命令行核对几何和属性完整性拿到shp后先不要急着拖进ArcGIS看图形用ogrinfo做一次快速体检更高效。ogrinfo是GDAL自带的命令行工具不需要打开界面就能输出图层范围、要素数量和字段定义ogrinfo -so -al 羌塘高原区.shp-so表示只读取概要信息-al表示遍历该文件中所有图层。输出结果里能看到要素数量、Extent范围以及每个字段的名称和类型。这一步能顺带确认海拔数据到底是存放在属性表里还是需要从高程TIF中提取。这份资源的海拔信息来自“全国高程数据”TIFshp本身是研究区边界面要素真正的海拔值要通过栅格提取获得也就是第四章用gdalwarp分割、第五章用rasterio统计的那条路径。提示如果ogrinfo报“Unable to open datasource”先检查.shp、.shx、.dbf三个文件是否在同一目录且文件名完全一致。3. mxd文件路径丢失排查与图层样式还原3.1 先解压所有压缩包再谈打开mxdmxd记录图层数据源时既支持绝对路径也支持相对路径但这处资源里的全国高程和行政区数据各自封装成了zip。只解压最外层压缩包时mxd引用的文件实际还躺在内部zip里打开后所有图层都会挂上红色感叹号。正确顺序是把资源zip解压至工作目录后再把TIF格式全国高程数据.zip和china-shapefiles.zip也分别解压出来让目录结构与mxd记录保持一致。解压完成后还要检查路径层级有些zip解压后会多出一层同名根目录导致mxd按原路径找不到文件。整个流程走完后用ArcMap打开羌塘高原区.mxd图层列表里没有红点才算通过。3.2 用ArcPy排查所有图层的真实路径手动右键Data - Set Data Source可以修但图层一多效率太低。用ArcGIS自带的Python环境跑一段脚本几十个图层的路径一次看完import arcpy mxd_path rD:\qiangtang\羌塘高原区.mxd mxd arcpy.mapping.MapDocument(mxd_path) for lyr in arcpy.mapping.ListLayers(mxd): if lyr.supports(DATASOURCE): print(lyr.name, -, lyr.dataSource)arcpy.mapping.MapDocument用来加载mxd工程文件ListLayers返回工程中所有图层对象supports(DATASOURCE)判断图层是否绑定数据源打印出的路径如果已经不存在说明该图层需要修复。批量修复场景下用replaceDataSourcemxd.replaceDataSource( rD:\qiangtang\merged, RASTER_WORKSPACE, None, True ) mxd.save()第一个参数是数据统一存放的目录第二个是工作空间类型栅格数据用RASTER_WORKSPACE矢量文件用SHAPEFILE_WORKSPACE或FILE_GDB_WORKSPACE第三个参数传None表示不替换图层名True表示验证路径存在性。执行完保存再打开mxd检查红感叹号。打开时现象原因处理图层红感叹号数据源路径失效Data - Set Data Source指回TIF显示为黑色块渲染未拉伸符号系统-拉伸-直方图中文注记变问号cpg丢失或字体缺失补cpg安装中文字体比例尺指北针错位布局视图模板被改动布局-恢复默认模板3.3 海拔分级符号化的可编辑点mxd的价值在于图层顺序和符号是调好的边界线宽、注记字体、海拔色带都在里面。直接使用省事但建议先复制一份工作副本保留原始mxd作为备份。做海拔专题图时常用做法是双击图层打开符号系统选择“分类-分级色彩”分级数量设在5到8级配色用渐变色带分类方法选Natural Breaks (Jenks)最后把图层透明度调到20%左右再与底图叠加。关键参数是分级方法与透明度分级方法决定海拔区间切在哪透明度影响TIF底图与边界线的关系。注记层如果与研究区边界重叠把注记字号调小并开启随比例尺缩放避免出图时字压线。4. TIF成图坐标对齐从投影参数到CAD落位4.1 TIF成图与tfw世界文件的关系标准成图TIF是带地理参考的栅格图片。ArcGIS导出TIF时通常会生成一个同名.tfw世界文件用六个数字描述影像左上角坐标和像元大小。用gdalinfo可以快速确认这张TIF是否具备坐标信息gdalinfo 羌塘高原区.tif输出结果里的Origin代表左上角坐标Pixel Size代表X与Y方向的像元尺寸。拿到数据后先跑这一条命令确认它是含坐标的GeoTIFF。如果输出里只能看到Image Structure Metadata而找不到Origin说明这张图只是普通图片改名必须经过Georeferencing手动配准才能与其他图层对齐。4.2 全国高程TIF与shp范围的精确裁剪全国高程TIF覆盖范围大直接铺进工程会拖慢显示做统计时也会把研究区外的像元算进去。用羌塘高原区shp做cutline来裁剪高程栅格gdalwarp -cutline 羌塘高原区.shp -crop_to_cutline \ -dstnodata -9999 -overwrite \ 全国高程.tif 羌塘裁剪.tif-cutline指定作为裁剪边界的矢量文件-crop_to_cutline让输出栅格范围与shp边界严格一致-dstnodata -9999把边界外像元统一写成无效值-overwrite允许覆盖已有输出文件。裁剪后产生的羌塘裁剪.tif范围、投影、像元尺寸都与shp匹配后续放进mxd做海拔统计或坡度分析不会出现范围错位。有人习惯用gdal_translate -projwin按矩形范围裁那种方式适合规则边界研究区是不规则多边形时-cutline才是正确选项否则会出现大块黑色背景进而污染平均海拔统计。4.3 TIF影像图在CAD类软件自动定位的原理经常有人问TIF影像图怎么在CAD中自动定位这个问题的前提是TIF带地理参考信息。CAD本身不识别投影坐标它靠tfw提供的仿射参数把影像放到世界坐标对应的图面位置。可行做法是确保TIF同一目录下有tfw文件插入影像时选择关联世界文件另一种做法是把TIF在GIS里导出为带坐标的DXF再进CAD。tfw的六行参数含义如下行数含义常见取值特征第1行X方向像元宽度正数值越小分辨率越高第2行X方向旋转项一般取0第3行Y方向旋转项一般取0第4行Y方向像元高度通常为负数第5行左上角像素中心X坐标与投影带有关第6行左上角像素中心Y坐标与纬度有关手动编辑tfw时必须保持六行结构删除任意一行或调整顺序都会让CAD插入后的位置和比例完全失控。如果插入后影像落在世界原点附近而不是目标位置优先检查tfw是否存在、坐标是否被下载过程重命名过。5. 行政区shp裁剪与批量出图的实用脚本5.1 用省级行政区数据搭邻域背景china-shapefiles.zip里是全国省级行政区划shp适合做研究区周边背景。裁剪之前先统一坐标系通常以shp为准再用空间筛选只保留与研究区相交的省份避免全图数据冗余。import geopandas as gpd pref gpd.read_file(china-shapefiles/china_province.shp) aoi gpd.read_file(羌塘高原区.shp).to_crs(pref.crs) near pref[pref.geometry.intersects(aoi.unary_union)] near.to_file(邻域省份.shp, encodingutf-8)intersects判断两个几何是否相交aoi.unary_union把研究区所有要素合并成单个几何对象避免逐个要素重复输出。编码参数指定为utf-8防止后续在QGIS中打开出现中文乱码。5.2 提取研究区平均海拔提取海拔的核心思路是用shp边界裁剪全国高程TIF然后对有效像元做统计import rasterio from rasterio.mask import mask import geopandas as gpd aoi gpd.read_file(羌塘高原区.shp) dem rasterio.open(全国高程.tif) out_img, out_transform mask( dem, aoi.geometry, cropTrue, nodata-9999 ) values out_img[out_img ! dem.nodata] print(平均海拔, values.mean().round(1))mask按矢量几何裁剪栅格cropTrue让返回的数组与变换矩阵贴合研究区范围nodata-9999必须与gdalwarp阶段的无效值保持一致否则无效像元会混进平均值计算。如果TIF与shp坐标系不一致要先转换shp到TIF的坐标系mask才会按正确的空间范围取样。5.3 导出草图验证渲染结果mxd里的版面调好后用ArcPy的ExportToPNG输出一组草图横向DPI设150就能满足预览需求。导出后先看边界是否出现白边再看海拔色带是否连续过渡边界异常优先回查4.2节的cutline参数色带断层则检查无效值设置确认裁剪TIF和原始高程TIF的nodata一致。本文还有配套的精品资源点击获取