Postgis使用工具raster2pgsql批量导入栅格数据（二）

另外一种批量导入栅格数据的方式，是使用工具gdalbuildvrt生成vrt文件 。

1.栅格数据

        WorldClim是一个高空间分辨率的全球天气和气候数据的数据库。可以从该网站下载天气数据。网址如下：https://www.worldclim.org。

下面是12个月份温度统计数据

 2.生成rvt文件

gdalbuildvrt -separate tmax_multi.vrt tmax*.bil

验证vrt文件的正确性

gdalinfo tmax_multi.vrt

结果如下：

Driver: VRT/Virtual Raster
Files: tmax_multi.vrt
       tmax1.bil
       tmax10.bil
       tmax11.bil
       tmax12.bil
       tmax2.bil
       tmax3.bil
       tmax4.bil
       tmax5.bil
       tmax6.bil
       tmax7.bil
       tmax8.bil
       tmax9.bil
Size is 2160, 900
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000000000000057,90.000000000000000)
Pixel Size = (0.166666666666667,-0.166666666666667)
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000)
Lower Left  (-180.0000000, -60.0000000)
Upper Right ( 180.0000000,  90.0000000)
Lower Right ( 180.0000000, -60.0000000)
Center      (   0.0000000,  15.0000000)
Band 1 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-478.000 Max=418.000
  NoData Value=-9999
Band 2 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-265.000 Max=401.000
  NoData Value=-9999
Band 3 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-373.000 Max=401.000
  NoData Value=-9999
Band 4 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-452.000 Max=413.000
  NoData Value=-9999
Band 5 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-421.000 Max=414.000
  NoData Value=-9999
Band 6 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-422.000 Max=420.000
  NoData Value=-9999
Band 7 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-335.000 Max=429.000
  NoData Value=-9999
Band 8 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-190.000 Max=441.000
  NoData Value=-9999
Band 9 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-94.000 Max=467.000
  NoData Value=-9999
Band 10 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-59.000 Max=489.000
  NoData Value=-9999
Band 11 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-76.000 Max=474.000
  NoData Value=-9999
Band 12 Block=128x128 Type=Int16, ColorInterp=Undefined
  Min=-153.000 Max=441.000
  NoData Value=-9999

3.使用工具raster2pgsql导入栅格数据

首先，生成sql文件

raster2pgsql  -d  -I -C  -M -F -t 100x100 -s 4326 tmax_multi.vrt  public.tmax_multi > tmax_multi.sql

然后，执行sql文件

psql -d postgis_32_sample -U postgres -f tmax_multi.sql

4.查询表raster_columns 的信息

SELECT * FROM raster_columns where r_table_name = 'tmax_multi';

在pgAdmin中的显示结果如下：

 

5.统计经度12.49，纬度41.88处每个月的平均温度

SELECT 
(ST_VALUE(rast,1, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As jan,
(ST_VALUE(rast,2, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As Feb,
(ST_VALUE(rast,3, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As mar,
(ST_VALUE(rast,4, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As apr,
(ST_VALUE(rast,5, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As may,
(ST_VALUE(rast,6, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As jun,
(ST_VALUE(rast,7, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As jul,
(ST_VALUE(rast,8, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As aug,
(ST_VALUE(rast,9, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As sep,
(ST_VALUE(rast,10, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As oct,
(ST_VALUE(rast,11, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As nov,
(ST_VALUE(rast,12, ST_SetSRID(ST_Point(12.49, 41.88),4326))/10) As dec
FROM tmax_multi where rid IN 
 (SELECT rid FROM tmax_multi WHERE ST_Intersects(ST_Envelope(rast),  ST_SetSRID(ST_Point(12.49, 41.88),4326)))

在pgAdmin中的显示结果如下：