一、 数据准备
从地理空间数据云下载 DEM 数据,选择 SRTM 90 米分辨率,选取覆盖湖北省范围,查询到共有 6 幅影像,也可以从别的网络渠道获取免费开源 DEM 数据。
二、 查看数据
下载完数据后使用 QGIS 打开 6 幅图像,设置渲染波段,下载的数据统一 WGS84 坐标系,通常就可以用GIS软件来做镶嵌、合并、裁剪等操作来抽取研究范围内的栅格数据,本系列文章主要用PostGIS来处理分析栅格影像数据,用于满足系统底层开发需求。
三、 数据导入
PostGIS 作为 PostgreSQL 数据库的一个扩展,让其具备了针对空间数据的存储、查询、分析能力,因此首先要安装好PostgreSQL数据库,本文已经创建好了【gisdata】数据库,可根据项目实际业务需求单独建立空间数据库,然后安装 PostGIS 扩展,为了支持栅格数据,还需要创建 postgis_raster 扩展。
- 创建用户、数据库
-- 创建用户、密码createusergisdatawithpassword'gisdata';-- 创建数据库createdatabasegisdatawithowner gisdata;-- 为用户授予数据库的权限grantallprivilegesondatabasegisdatatogisdata;-- 授予所有表操作权限grantallprivilegesonalltablesinschemapublictogisdata;-- 切换到新创建的数据库\c gisdata-- 查看扩展列表\dx- 查询是否已经安装扩展
-- 查看可用扩展,是否已安装postgis,postgis_rasterselect*frompg_available_extensions;-- 如果未安装,在开启扩展创建postgis_raster扩展createextension postgis_raster;- 上传文件至数据库服务器
注意:当前数据库用户需要有文件的读写权限
- 批量导入数据
- 当前执行用户需要有raster2pgsql执行权限,要安装配置好postgis
- 执行命令的路径在文件所在路径
- 导入前原始影像文件坐标系保持一致,如果不一致会报错,因为有空间坐标约束,但是也可以正常导入,但是强烈建议不要这么做
/** * 坐标系未指定,默认源文件坐标系 * -I 创建索引 * -C 添加约束 * -M 导入后执行vacuum analyze * -P 保证瓦片的大小一致 * -F 保留文件名 * -t 瓦片大小(宽和高)为256x256 * *.img 导入所有img文件 * public.srtm_90 生成一个在public模式下的srtm_90表 * gisdata导入的数据库 */raster2pgsql-I-C-M-P-F-t256x256*.imgpublic.srtm_90|psql-d gisdata注意:按照下载下来的分块的DEM影像数据,导入完成后,将每一个栅格文件按照指定瓦片大小(宽和高)进行切分,每一行数据是一个rast数据,如果不指定切片大小,则按照默认一个文件为一个rast数据.
四、 基本查询
数据入库完毕后,会将生成的表注册到系统视图【raster_columns】中,该表维护了栅格数据的基本元数据信息,可以看到包括坐标系、分辨率、瓦片大小、栅格数据类型、空间范围等等详细信息
一般在实践中,在进行栅格数据的空间查询、裁剪、转换等操作之前,先通过 【raster_columns】视图获取相关栅格的信息,有助于确定合适的操作方法和参数。
数据导入后生成一个 【public.srtm_90】 表,查询表的基本情况,为后续的栅格数据查询分析提供基础
-- 1. 数据注册查询select*fromraster_columns;-- 返回结果:Name|Value|----------------+-----------------------------------------------------------------+r_table_catalog|gisdata|r_table_schema|public|r_table_name|srtm_90|r_raster_column|rast|srid|4326|scale_x|0.0008333333|scale_y|-0.0008333333|blocksize_x|256|blocksize_y|256|same_alignment|true|regular_blocking|false|num_bands|1|pixel_types|{16BSI}|nodata_values|{-32768.0}|out_db|{false}|extent|POLYGON((10524.88,10535,120.1235,120.1224.88,10524.88))|spatial_index|true|-- 2. 查询记录,理解每个字段意义有助于更好的理解对栅格表的处理分析,尤其是rast字段,是Postgis处理栅格的基本核心select*fromsrtm_90limit3;-- 返回结果-- rid:瓦片行id,-- rast:代表一个被切成的256x256栅格瓦片,-- filename:瓦片所属的文件rid|rast|filename|-----+-------------------+--------------+1|01000001004F...77A0|srtm_58_06.img|2|01000001004F...7FA0|srtm_58_06.img|3|01000001004F...6D80|srtm_58_06.img|-- 3. 查看元数据信息select(ST_MetaData(rast)).*fromsrtm_90whererid=1;-- 返回结果,返回瓦片的基本信息:Name|Value|----------+-------------------+upperleftx|114.12458333333363|upperlefty|31.362916666666656|width|256|height|256|scalex|0.0002777777777778|scaley|-0.0002777777777778|skewx|0.0|skewy|0.0|srid|4326|numbands|1|-- 4. 瓦片基本统计信息,均值、方差、最值等select(ST_SummaryStats(rast,1,true)).*fromsrtm_90whererid=1;-- 返回结果,单个瓦片信息的统计信息,因为是256x256大小,总像素数为65536count|sum|mean|stddev|min|max|-----+---------+-----------------+------------------+------+------+65536|114254751|1743.389144897461|134.71134814142528|1392.0|2193.0|-- 5. 按文件名分组查询每个文件被分成了多少瓦片selectfilename,count(1)fromsrtm_90groupbyfilename;-- 返回结果,基本上每个文件被分成了576个256x256大小的瓦片:filename|count|--------------+-----+srtm_59_06.img|576|srtm_58_07.img|576|srtm_60_07.img|569|srtm_60_06.img|572|srtm_58_06.img|576|srtm_59_07.img|576|--6. 按像素值统计个数select(ST_ValueCount(rast)).*fromsrtm_90whererid=1limit10-- 返回结果:按照dem值统计个数value|count|------+-----+1954.0|66|1939.0|90|1940.0|79|1948.0|84|1969.0|48|1975.0|50|1935.0|101|1913.0|110|1892.0|126|1881.0|134|