使用Landsat和PostGIS按地理位置检索光栅数据

我正在从事的项目要求我检索特定地理位置(lon/lat)的陆地卫星光栅数据。在筛选了一些教程并试用了GDAL、PostGIS和QGIS之后,我成功地将GeoTIFF Landsat图像导入到PostGIS光栅表中,并从该表中按地理位置访问值。然而,结果中存在一些问题:

• 我不了解QGIS在其界面中使用的坐标系,因为它们的范围是十万个
• 光栅加载到西班牙海岸以外的QGIS中,而不是在美国缅因州的顶部。

以下是有关我的流程的一些信息。总的来说,我对地理信息系统相当陌生,因此我几乎可以肯定,在这里会发现一个明显的错误:

• 从USGS GloVis下载Landsat 8 GeoTIFF文件
• 将5级图像重命名为更友好的图像,以便于指挥忍者。
• 为光栅表创建postgres数据库并运行 CREATE EXTENSION postgis;

• 跑 gdalinfo LSSampleB5.TIF ,打印以下输出:

  Driver: GTiff/GeoTIFF Files: LSSampleB5Test2.TIF Size is 7871, 7971 Coordinate System is: PROJCS["WGS 84 / UTM zone 19N", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-69], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["Easting",EAST], AXIS["Northing",NORTH], AUTHORITY["EPSG","32619"]] Origin = (318285.000000000000000,5216715.000000000000000) Pixel Size = (30.000000000000000,-30.000000000000000) Metadata: AREA_OR_POINT=Point Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 318285.000, 5216715.000) ( 71d23'37.53"W, 47d 4'44.12"N) Lower Left ( 318285.000, 4977585.000) ( 71d18' 9.77"W, 44d55'42.53"N) Upper Right ( 554415.000, 5216715.000) ( 68d16'58.41"W, 47d 6' 6.11"N) Lower Right ( 554415.000, 4977585.000) ( 68d18'36.69"W, 44d56'58.62"N) Center ( 436350.000, 5097150.000) ( 69d49'20.56"W, 46d 1'29.87"N) Band 1 Block=7871x1 Type=UInt16, ColorInterp=Gray

• 我将此输出解释为EPSG 4326格式(这可能是我的错),因此我运行以下命令将GeoTIFF作为PostGIS光栅导入:

  raster2pgsql -s 4326 -I LSSampleB5.TIF -F -t 50x50 -d | psql -U postgres rastertest

• 这成功导入了一个新表。然后,我使用QGIS获得了对正在发生的事情的视觉直觉。

• 在下面 Database -> DB Manager -> PostGIS -> rastertest -> public 我将lssampleb5添加到画布中。

• 我在QGIS中创建了一个新的XYZ连接,以添加Google satillite混合图像供参考。我使用的url是 https://mt1.google.com/vt/lyrs=y&x= {x}&y={y}&z={z} 最小和最大缩放分别为0和19。

• 在这里,我注意到了一个事实,即在谷歌混合地图上,lssample层降落在西班牙海岸附近。

• 我确保两层都在EPSG 4326投影上,没有变化。

• 我并没有气馁地继续,我尝试了一个数据库查询来获取单个像素值。由于我的样本数据降落在西班牙附近,因此我使用QGIS对附近的有效坐标对进行了采样以进行查询。查询是:

  SELECT rid, ST_Value(rast, 1, ST_SetSRID(ST_Point(448956,5041439), 4326)) as b5 FROM lssampleb5 WHERE ST_Intersects(rast, ST_SetSRID(ST_Point(448956,5041439), 4326)::geometry, 1);

• 这返回了有效的行ID和5776的ST\U值。尝试超出QGIS显示的范围的坐标时,没有返回任何条目,这并非意外。

所以,首先,我不知道QGIS在坐标系中使用的是什么。这绝对不是原始形式的经纬度,但根据我的理解,EPSG 4326应该是一个地理投影。

第二,我不知道为什么QGIS把陆地卫星场景放错了地方,或者在这个过程中场景没有正确转换。