1. 先搞清楚为什么要用PostGIS去取一个像素值做GIS的人都干过这种事手里拿到一份栅格数据比如30米分辨率的DEM、坡度图、NDVI影像然后你需要知道某个点位、某条线路或某个地块上像素值到底是多少。以前最常见的做法是什么用ArcMap的识别工具逐点去点或者用多值提取至点跑一遍再把属性表导成Excel慢慢看。点位少还好说点一多或者你需要在PostgreSQL里跟业务表做关联查询、做统计分析这套流程就非常痛苦了。PostGIS里有个ST_Value函数就是专门干这个的给定一张栅格表和一组几何坐标位置直接返回那个位置对应的像素值。这样就能把栅格取值直接塞进SQL里跟矢量表做JOIN跟业务数据做联动统计还能写成存储过程自动化处理。这篇文章我会把栅格数据载入、ST_Value的参数用法、批量取值、性能优化和常见坑全部拆开讲适合刚接触PostGIS栅格功能或者在ArcMap跟PostGIS之间切换工作流的朋友参考。先说清楚一个核心概念ST_Value本身并不复杂复杂的是它背后涉及到栅格数据组织形式、坐标系、波段、NODATA值、采样方式这些容易被忽略的细节。很多人用这个函数返回NULL或者数值明显不对十有八九不是函数用错而是数据准备阶段出了问题。所以我不会一上来就丢函数文档而是从数据准备开始把整个链路串起来讲。1.1 平时大家都在用什么方式取像素值在我接触过的团队里取栅格像素值的方案大致可以分为三类各有各的问题。第一类是手工点选在ArcMap里加载栅格用识别工具一个个点点完记录在Excel里。这个方法只适合几个点几十个点以上纯粹是折磨人而且完全没有留痕机制操作过程不可追溯。第二类是脚本批量提取。ArcMap里的提取多值至点工具或者Python的arcpy.Sample函数都能批量把栅格值提取到点要素的属性表里。这种方法适合一次性处理但毛病也很明显每次处理都要打开桌面软件数据无法实时联动如果你的点表是天天在变的业务数据总不能天天跑一遍手动导出吧。第三类是我现在推荐的做法——直接在数据库里算。把栅格载入PostGIS点表作为PostGIS的矢量表一句SQL就可以把几百个点的栅格值全部取出来。好处至少有三点一是过程可重复SQL写一次以后数据更新了重跑一遍就行二是可以跟其他表做复杂的JOIN和统计比如按行政区汇总平均坡度、按地块提取像元值后再做回归分析三是原生跑在服务器上处理几十万点也不会累死你的电脑。1.2 ST_Value到底解决什么问题ST_Value在PostGIS栅格功能里属于基础函数它的定位就是从raster对象里在指定位置取出对应像素值。注意我说的是raster对象不是raster文件。PostGIS里一张栅格表存储的是若干条raster记录每条记录是一个瓦片需要ST_Value在这个瓦片上定位然后返回值。为什么单独一个取值函数值得专门写一篇因为取值这个动作一旦放到数据库环境里就涉及三个问题位置用什么表达几何还是坐标、坐标系要不要转换、栅格外的点怎么处理。这三个问题不解决你写出来的SQL看着没问题结果就是不对。ST_Value本身是个很轻的函数但用对它的前提是理解整个栅格数据模型这也是很多人觉得PostGIS栅格功能难上手的原因。1.3 哪些业务场景用得上ST_Value从我的实践经验看这套能力的应用面比想象中广。最典型的是农业和林业的样地调查外业人员带GPS采了一堆样点需要把这些样点对应的高程、坡度、土壤属性栅格值自动提取出来再做后续建模分析。其次是环保和水利的断面监测河流断面点位需要提取多年遥感影像的NDVI或水质反演值做时间序列分析SQL里写个循环就能批量完成。还有城市规划和交通领域的沿线分析一条规划路线做缓冲区后需要统计缓冲区内的平均高程、坡向分布甚至要按像素值做成本计算。这类任务在ArcMap里做要导来导去在PostGIS里一条统计SQL就出来了。2. 准备工作把栅格塞进PostGIS再说所有ST_Value使用的前提是栅格数据已经正确载入到PostGIS数据库里。这一步很多朋友草草用工具点了导入就完事结果切片方式不对、坐标系没指定、NODATA值丢了后面全乱套。所以这个环节值得花点篇幅讲透。2.1 raster2pgsql载入栅格的基本姿势栅格载入PostGIS的主力工具是raster2pgsql它是随PostGIS一起装的可执行文件。我见过有些朋友试图用SQL的方式手动INSERT栅格数据实际上在大多数场景下完全没必要raster2pgsql生成的SQL脚本已经足够灵活。最基础的一条命令长这样raster2pgsql -s 4326 -I -C -M -t 128x128 -F -n name ./dem.tif public.dem | psql -U postgres -d gisdb这条命令里的参数每个都有讲究。-s 4326是指定SRID这里用4326意思是数据本身已经是WGS84经纬度坐标。如果数据是投影坐标系比如UTM或者高斯克吕格强烈建议先搞清楚原始坐标系再指定千万别乱填否则你后面用ST_Value取出来的坐标点全是错位的。-I是创建GiST索引对后面的空间查询性能至关重要。-C是添加栅格约束包括SRID、像素类型、NODATA值、切片范围等。-M表示对整个栅格文件生成vacuum analyze信息让查询计划器心里有数。-t 128x128是按这个瓦片大小切块。至于为什么切块下一节详细说。加一个-F和-n参数是把文件名作为一个属性列存起来当你需要同时管理多张栅格时非常有用。如果你不想在命令行里处理也可以用PostGIS自带的图形化工具shp2pgsql-gui里附带的功能或者用QGIS的DB Manager导入栅格。但命令行是最不容易出错的方式而且脚本化之后可以反复用强烈建议习惯这种方式。2.2 栅格数据组织形式切块、金字塔与元数据这一节是整个PostGIS栅格功能最容易让新手懵圈的地方。很多人下意识觉得栅格表里就是一大块完整的栅格记录其实不是。raster2pgsql默认会把整张栅格按指定的尺寸切成一个个瓦片每片作为一行存储。为什么要切块因为PostgreSQL的数据模型里一个raster对象是作为一个整体存储在一条记录里的如果一整幅全国30米分辨率DEM是几千万个像元数据库根本没法高效读取。切成128x128或者256x256后每次只读取涉及查询范围的那几条记录配合GiST索引精确找到命中的瓦片查询效率能提升几个数量级。你可以把切块想象成把一本书拆成可以单页检索的活页用ST_Value取值时只需要翻到对应那页就行了。切块之后还有一个重要的元数据视图raster_columns。你在创建扩展并载入栅格后可以通过这个视图查看所有栅格列的资料比如SRID、瓦片大小、波段数量、NODATA值。用这样一句就能看核心信息SELECT r_table_name, r_raster_column, srid, pixel_types, nodata_values, tile_size_x, tile_size_y FROM raster_columns;我调数据问题时的第一件事永远是先查这个视图确认SRID和NODATA有没有配对。很多人ST_Value取出来的结果是NODATA值而不是真实数值十有八九就是载入时NODATA信息丢了或者设错了。关于金字塔再多说一句。PostGIS的栅格类型本身不做金字塔所谓概述表的优化手段也是全局性的。如果你的需求是快速做大范围浏览可以在raster2pgsql里用-l参数生成概述表比如-r降低分辨率到原来的四分之一。但ST_Value是精确定位取值走的是原始分辨率数据跟概述表关系不大所以做不做只看你的可视化需要就行。2.3 PostGIS Raster安装失败的常见坑既然热词里有PostGIS安装失败那这个坑确实值得单说。很多朋友装完PostGIS正常建库create extension postgis也成功了但执行和栅格相关的SQL时报错说找不到函数或者创建postgis_raster扩展时说control file not found。这是因为从PostGIS 3.0开始栅格功能从核心扩展里拆分成了独立的postgis_raster扩展。你得确认安装的时候是否勾选了栅格组件。Windows下的安装包有PostGIS Raster选项Linux下一般是postgis package的分包比如postgis-3.x。装完之后要在数据库里单独执行CREATE EXTENSION IF NOT EXISTS postgis_raster;如果报could not open extension control file那你大概率只是装了核心PostGIS没装栅格组件。Debian/Ubuntu系可以用apt install postgis postgis-docCentOS系要看repo里有没有postgis的raster子包。还有种情况是装了旧版PostGIS比如2.x栅格功能是包含在postgis扩展里的不需要单独创建。所以先了解你装的版本再决定怎么做别一上来就蒙头执行命令。3. ST_Value函数用法拆解参数、坐标与采样那些事ST_Value的官方文档其实很短但实际使用中牵扯到的细节很多。这一节我把几种调用方式、坐标转换、波段和采样参数全部展开讲。3.1 最常用的三种调用写法ST_Value最常见的写法有以下几种选哪种取决于你的数据形态。第一种传入栅格对象和像素行列号。这是最原始的定位方式你直接告诉它第几行第几列它返回那个像素的值SELECT ST_Value(rast, 1, 100, 200) AS val FROM dem WHERE rid 1;这里的参数分别是raster对象、波段索引、列号、行号。现实中用得不多因为你通常没有行列号只有经纬度或投影坐标但理解它是理解后面函数的基础。第二种传入栅格对象和点几何SELECT ST_Value(rast, (SELECT geom FROM sample_points WHERE id 1)) AS val FROM dem WHERE ST_Intersects(rast::geometry, (SELECT geom FROM sample_points WHERE id 1));这个写法是把点几何传给ST_Value内部会帮你定位但前提是栅格和点的坐标系一致。如果点列表很大效率上还有优化空间后面专门说。第三种传入栅格对象和坐标数值SELECT ST_Value(rast, 1, ST_SetSRID(ST_MakePoint(120.5, 30.8), 4326)) AS val FROM dem WHERE ST_Intersects(rast::geometry, ST_SetSRID(ST_MakePoint(120.5, 30.8), 4326));这种适合临时查一下某个坐标值不用建表。在写自动化脚本时也很实用比如从业务表里取出经纬度字段动态构建点再去取值。需要特别注意的是ST_Value省略波段参数时默认取波段1。如果你的栅格是多波段的影像想取第二波段的值一定要把波段索引写清楚。还有一种常见的错误是把栅格本身传错位置写成ST_Value(geom, rast)参数顺序搞反函数直接报错或者行为异常。3.2 坐标系不一致是最大的坑ST_Value内部不会做坐标转换。它拿到点几何之后直接把点坐标对标到栅格的地理位置如果点坐标系和栅格坐标系不一样轻则取到错误的像素值重则定位到栅格范围之外直接返回NULL。我举个例子。一张SRID为32650的DEM栅格UTM 50N投影点表是SRID 4326经纬度。如果你不做任何处理直接查ST_Value会把经纬度坐标当成投影坐标去栅格上定位。经纬度的数值范围是三十几度、一百二十几度投影坐标是几百万米量级定位结果自然乱套。正确的做法是先用ST_Transform把点转到栅格的SRID或者反过来把栅格转过去。实际工程中我强烈建议转点不转栅格因为栅格可能很大ST_Transform栅格的重采样开销远大于转换一个点。写法如下WITH point_geom AS ( SELECT ST_Transform(geom, 32650) AS geom FROM sample_points WHERE id 1 ) SELECT ST_Value(r.rast, p.geom) AS val FROM dem r, point_geom p WHERE ST_Intersects(r.rast::geometry, p.geom);这里ST_Transform也需要目标SRID正确无误。很多人以为只要转换了就行但转换时没指定正确的目标SRID照样出错。建议每次写取值SQL前先在raster_columns视图里确认栅格的SRID再动手。还有一个容易忽略的细节如果栅格是4326经纬度坐标系点的X是经度、Y是纬度那么点的坐标顺序必须正确。PostGIS遵循传统的X、Y顺序也就是经度在前纬度在后别跟某些软件里纬度、经度的显示顺序搞混。3.3 NODATA、波段与采样方式怎么选栅格数据里通常会有一个NODATA值表示这个像素没有有效数据。比如DEM的边缘、云遮挡区域、统计范围外的区域。ST_Value取到NODATA像素时不会自动帮你过滤它返回的就是NODATA值本身。所以如果你想统计有效像素的数量或者计算均值必须先排除NODATA。假设栅格表的NODATA值是-9999典型写法SELECT COUNT(*) AS valid_count, AVG(val) AS avg_val FROM ( SELECT ST_Value(r.rast, ST_Transform(p.geom, 4326)) AS val FROM dem r, sample_points p WHERE ST_Intersects(r.rast::geometry, ST_Transform(p.geom, 4326)) ) t WHERE val IS NOT NULL AND val -9999;注意我在这里对点做了ST_Transform因为DEM是4326。如果你已经预先在点表里把坐标转好就不需要每次查询都转了。还有一点ST_Value返回类型是double precision理论上不会返回NULL只会返回NODATA值所以NULL判断往往是多余的但留着没坏处万一你手动修改过栅格数据呢。像素类型也值得留意。GDT_Byte类型的栅格取值是0到255的整数GDT_Float32是浮点。ST_Value统一返回浮点数你可以用ST_BandPixelType来确认栅格类型。在做地形分析时如果发现提取出来的坡度值范围正常但小数部分全是整数的样子先怀疑是不是像素类型本身就是短整型存储。最后说采样方式。当点位置恰好在像素边缘时ST_Value默认采用最近邻采样也就是取坐标所在像素的值不做插值。这个行为对大多数“提取对应位置”的需求是合理的因为栅格本身就是离散的。但如果你希望更平滑的结果可以用ST_Value的重采样变体比如双线性、三次卷积。不过需要注意变体的用法因PostGIS版本而异至少在2.5之后可以传入采样算法参数。实际工作中除非做连续变量插值分析否则默认最近邻就足够了。4. 批量取值实战从ArcMap点表到PostGIS输出前面的基础都准备好了这一节进入真正能落地的场景。我拿“样地坡度提取”作为完整案例一步步演示从ArcMap导出点表、载入PostGIS、关联提取栅格值、结果校验的全过程。这个案例非常典型几乎涵盖ST_Value批量使用的主要环节包括坐标转换、空间连接、性能优化和结果对比。4.1 场景设定一次地形样地坡度提取假设我在某个山区做过实地样地调查外业收了200多个样地点存在Excel表格里字段包括样地编号、经度、纬度。手里有一张坡度栅格SRID是32650UTM投影30米分辨率NODATA值是-9999。任务很明确算出每个样地点的坡度值用来做后续的物种分布模型。如果你从ArcMap那边接手通常的做法是先在ArcMap里把点坐标用添加XY数据变成点要素然后导出成shp或gdb要素类。这一步我不建议直接用数据库的几何生成工具去建。在实际工作流里ArcMap在坐标转换和符号化上依然方便所以数据可以在ArcMap里准备好但最终成果要放数据库。4.2 建点表、载入数据、关联取值首先在PostGIS里建一张样地表字段设计好CREATE TABLE sample_points ( id SERIAL PRIMARY KEY, sample_no TEXT, geom_4326 GEOMETRY(Point, 4326), geom_utm GEOMETRY(Point, 32650) );这里我刻意保留了两个坐标系的字段。geom_4326用来跟业务系统对接geom_utm用来跟栅格做空间关联。虽然你也可以在查询时动态转换但批量执行时预转换显然更高效而且建表的时候一次性算好SQL写起来也更清爽。然后通过shp导入工具shp2pgsql或者直接SQL编写INSERT把点数据导进去。坐标转换一步到位UPDATE sample_points SET geom_utm ST_Transform(geom_4326, 32650);接下来就是核心的关联查询。用ST_Value配合ST_Intersects限定范围确保只在栅格覆盖范围内取数SELECT p.id AS sample_id, ST_Value(r.rast, 1, p.geom_utm) AS slope_val FROM sample_points p JOIN dem_slope r ON ST_Intersects(r.rast::geometry, p.geom_utm) ORDER BY p.id;很多朋友问为什么查询慢或者说返回的行数跟点数对不上。先说行数问题如果某个点落在栅格覆盖范围之外JOIN就匹配不到任何瓦片这条记录就消失了。你要是强制保留所有点就得用LEFT JOINSELECT p.id AS sample_id, ST_Value(r.rast, 1, p.geom_utm) AS slope_val FROM sample_points p LEFT JOIN dem_slope r ON ST_Intersects(r.rast::geometry, p.geom_utm) ORDER BY p.id;这样写之后范围外的点就会返回NULL或NODATA值。我一般倾向用LEFT JOIN宁可后续标记为无效也不愿意静默丢失数据。实际处理时还会遇到一个点在多个瓦片上都命中的情况。如果瓦片之间有重叠或者ST_Intersects的范围跨越多个切片JOIN会产生重复行。这时候要么用DISTINCT ON要么用LATERAL子查询精确控制每个点只取一行。4.3 性能优化执行计划与索引检查好如果数据量到了十几万点上面这种直接JOIN的写法大概率会慢得让你怀疑人生。ST_Intersects写的条件虽然能用上栅格的GiST索引但ST_Value本身是在每条命中的记录上再做一次像素定位。糟糕的执行顺序会导致栅格扫描次数爆炸。我在生产环境里更常用的写法是LATERAL子查询配合LIMIT 1让每个点只取一个栅格瓦片SELECT p.id AS sample_id, s.slope_val FROM sample_points p LEFT JOIN LATERAL ( SELECT ST_Value(r.rast, 1, p.geom_utm) AS slope_val FROM dem_slope r WHERE ST_Intersects(r.rast::geometry, p.geom_utm) ORDER BY ST_Distance(r.rast::geometry, p.geom_utm) LIMIT 1 ) s ON true ORDER BY p.id;这个写法的主要思路是对于每个点数据库只去找最近的一个栅格瓦片然后取这个瓦片上的像素值避免重复匹配和重叠瓦片干扰。ORDER BY ST_Distance在这种场景下看起来奇怪但它保证了取到的是最接近点的那个切片在有重叠瓦片或者边缘跨片时结果更可靠。实测下来十几万点的量级也能在几秒到几十秒内完成瓶颈主要在栅格瓦片的I/O上。还有两个实用技巧。第一个是确认栅格表有没有GiST索引SELECT indexname, indexdef FROM pg_indexes WHERE tablename dem_slope;如果没有手动创建CREATE INDEX ON dem_slope USING gist (st_convexhull(rast));第二个技巧是缩小范围。如果你的栅格数据全国铺开而样地点集中在几个省可以用ST_Intersects先做一次粗筛再对点做缓冲区判断避免全表扫描。4.4 结果校验和ArcMap多值提取结果对比SQL跑完了怎么确认结果是对的最可靠的方式是拿ArcMap的“提取多值至点”工具结果做交叉比对。这是很多人的盲区——写完了代码执行结果是个数但对不对心里没底。我的做法是在ArcMap里加载坡度栅格和样地点要素用多值提取至点工具输出一个新的点要素类导出属性表为dbf或Excel。然后把ArcMap提取的坡度值和PostGIS提取结果放到同一张表里对比。通常你会发现大部分点是一致的但总有少数点有差异。差异的来源一般有两个一是坐标精度二是数据版本。坐标精度问题很好理解ArcMap对点的定位是显示坐标精度的四舍五入PostGIS直接用几何对象定位两者可能相差亚像素级。数值差异只要在合理范围内就没问题。数据版本问题就严重了如果ArcMap里的栅格文件是旧版本而你载入PostGIS的是更新后的版本结果自然不同。所以做对比前先核对数据时间戳和来源。有一个经验可以省很多事不要拿所有点完全一致作为目标而应该看差异分布。差异在个位数百分比以内说明SQL没问题如果大面积都是某个固定偏移量比如所有点都差一个像元八成是坐标系或者栅格对齐出了问题从数据准备阶段排查。5. 常见问题速查与避坑清单写这一部分的时候我把这几年在技术群里看到的高频问题汇总了一下特别是那些有明确典型特征的。这些问题单个看都不起眼但积累起来非常耗时间。5.1 为什么ST_Value返回NULLST_Value返回NULL的情况主要有四种。第一点是空几何或坐标缺失这种情况多发生在数据导入时没有清洗数据。第二点位于栅格范围之外。第三栅格瓦片本身存在缝隙投影边缘或者相邻瓦片没覆盖到。第四你查询的时候没有LEFT JOIN点没有匹配到任何瓦片。遇到NULL第一步别急着看函数先去检查点几何本身是不是有效的ST_IsValid然后再用ST_Intersects逐步验证。5.2 怎样把提取结果导出成Excel之前热词里有ArcMap栅格数据转化导出为Excel这里就顺带说清楚。PostGIS的结果在数据库里但最终分析报告通常要求Excel。最简单的办法是psql的COPY指令COPY ( SELECT p.id AS sample_id, p.sample_no, ST_Value(r.rast, 1, p.geom_utm) AS slope_val FROM sample_points p LEFT JOIN dem_slope r ON ST_Intersects(r.rast::geometry, p.geom_utm) ORDER BY p.id ) TO D:/slope_result.csv WITH CSV HEADER;生成CSV后用Excel打开另存为xlsx即可。如果不想用绝对路径也可以用pgAdmin的查询结果导出功能在查询结果面板右键导出选CSV或Excel格式。5.3 几个容易忽略的小细节第一点是关于ST_Value和栅格字段名。如果你的表有多个栅格列写函数的时候要写清楚用哪一列否则函数会拿默认列。第二点是关于栅格切块的边缘。当点正好落在瓦片边界时由于像素对齐的问题取值可能会受到相邻瓦片的影响这也是为什么我建议用LATERAL加ORDER BY ST_Distance的方式。第三点波段索引从1开始而不是0这个很多人第一次用的时候会踩坑。还有一个细节关于内存和临时文件。如果你的点表有几十万行JOIN产生的中间结果集可能非常大数据库需要临时空间。这种情况下需要关注temp_tablespaces的配置以及work_mem的设置。实测中一个五十万点的查询如果work_mem过小磁盘排序会拖垮整体性能至少把work_mem调到64MB以上并给数据库预留足够临时空间。6. 写在最后的一点体会这几次做ST_Value的实战下来我最深的体会是真正花时间的往往不是SQL怎么写而是数据本身是否规范。栅格数据来源多样坐标系五花八门NODATA值千奇百怪切块方式也各不相同每一个环节出问题最后体现在ST_Value上的结果就是在一堆NULL和错位值中排查。所以现在我拿到一个新任务第一件事不是写SQL而是先花时间把栅格数据的元数据、坐标系、NODATA、波段类型全部摸清楚再决定查询策略。我一直建议常用这套工作流的朋友把raster2pgsql的载入命令和校验SQL写成一个固定的脚本模板数据更新时一键重跑能省掉大量重复劳动。另外如果你经常从ArcMap那边收到点表可以考虑直接在PostGIS里用ST_Transform统一坐标别再两边来回导。这套流程稳定跑通之后我相信你也会和我一样彻底把取像素值这件事从桌面软件搬到数据库里。