另外一种批量导入栅格数据的方式,是使用工具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中的显示结果如下:
|