0

我正在进行的项目要求我在特定地理（经度/纬度）位置检索 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 卫星混合图像以供参考。我使用的 urlhttps://mt1.google.com/vt/lyrs=y&x={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 场景放错了位置，或者在这个过程中场景没有正确转换。

4

1 回答 1

0

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

在这里为您提供帮助：

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