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

Question

我正在进行的项目要求我在特定地理（经度/纬度）位置检索 Landsat 栅格数据。在筛选了一些教程并尝试了 GDAL、PostGIS 和 QGIS 之后，我成功地将 GeoTIFF Landsat 图像导入到 PostGIS 栅格表中，并从该表中按地理位置访问了值。但是，结果中存在一些问题：

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

以下是有关我的流程的一些信息。总的来说，我对 GIS 相当陌生，所以我几乎可以肯定这里有一个明显的错误：

• 从 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 卫星混合图像以供参考。我使用的网址是 @987654321@{x}&y={y}&z={z}，最小和最大缩放分别为 0 和 19。

• 这是我注意到 lssample 图层在 Google Hybrid 地图上落在西班牙海岸以外的地方。

• 我确保两层都在 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_VALUE。尝试 QGIS 显示范围之外的坐标导致没有返回条目，这并不意外。

所以，首先，我不知道 QGIS 使用的是什么坐标系。绝对不是原始形式的经度和纬度，但据我了解，EPSG 4326 应该是地理投影。

其次，我不知道为什么 QGIS 将 Landsat 场景放错了位置，或者在处理过程中场景没有正确转换的地方。

Answer 1

加入我们GIS SE，这是 GIS 相关问答的地方！

在这里为您提供帮助：

• 确实，您的罪行是 CRS。顶级PROJCRS 标签是 这里的关键，它从数据中读出"WGS 84 / UTM zone 19N"，用 底部的 EPSG 参考 (AUTHORITY["EPSG","32619"]])。
  EPSG:32619 是基于 WGS84 大地水准面（基准）的UTM projected CRS，单位为米，定义为到相应参考子午线（东经）的投影距离和赤道（北向）。由于您在导入期间定义了错误的 CRS（即 EPSG:4326），因此栅格的固有坐标值被视为度数，并将整个事物放置在世界的另一端。运行 UpdateRasterSRID (SELECT UpdateRasterSRID(<shema_name>, <your_raster_table>, rast, 32619);) 将栅格元数据设置为正确的 CRS 并重新加载图层。
• 至于 QGIS：它使用您告诉它使用的 CRS。 QGIS 带有一个非常方便的 on-the-fly 重投影功能（OTF，请查看“使用投影”here 的一般手册页），它可以让您定义用于投影和显示数据的任意 CRS（即，它将数据的 CRS 重新投影到内存中定义的 CRS 中，数据的元数据保持不变）。
  您可以找到 OTF 的快速链接按钮 em> 设置在 GUI 的右下角；将其设置为所需的 SRID（例如 4326）（您会注意到数据的视觉表示如何根据所选投影发生变化。显示的坐标也将使用 CRS 单位，例如 WGS84 的十进制度数）。