【应用实践】湖北省 DEM 数据下载导入PostGIS并实现基本查询

一、 数据准备从地理空间数据云下载 DEM 数据选择 SRTM 90 米分辨率选取覆盖湖北省范围查询到共有 6 幅影像也可以从别的网络渠道获取免费开源 DEM 数据。二、 查看数据下载完数据后使用 QGIS 打开 6 幅图像设置渲染波段下载的数据统一 WGS84 坐标系通常就可以用GIS软件来做镶嵌、合并、裁剪等操作来抽取研究范围内的栅格数据本系列文章主要用PostGIS来处理分析栅格影像数据用于满足系统底层开发需求。三、 数据导入PostGIS 作为 PostgreSQL 数据库的一个扩展让其具备了针对空间数据的存储、查询、分析能力因此首先要安装好PostgreSQL数据库本文已经创建好了【gisdata】数据库可根据项目实际业务需求单独建立空间数据库然后安装 PostGIS 扩展为了支持栅格数据还需要创建 postgis_raster 扩展。创建用户、数据库-- 创建用户、密码createusergisdatawithpasswordgisdata;-- 创建数据库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_90whererid1;-- 返回结果返回瓦片的基本信息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_90whererid1;-- 返回结果单个瓦片信息的统计信息因为是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_90whererid1limit10-- 返回结果按照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|