使用 Landsat 和 PostGIS 按地理位置检索栅格数据
Retrieving raster data by geographic location using Landsat and PostGIS
我正在从事的项目要求我检索特定地理 (lon/lat) 位置的 Landsat 栅格数据。在筛选了一些教程并尝试使用 GDAL、PostGIS 和 QGIS 之后，我成功地将 GeoTIFF Landsat 图像导入到 PostGIS 栅格 table 中，并从 table 按地理位置访问值。但是，结果中存在一些问题：
- 我不明白 QGIS 在其界面中使用的坐标系，因为它们 运行 有数十万个
- 栅格加载到西班牙海岸外的 QGIS 中，而不是原先应该加载到美国缅因州的顶部。
这里有一些关于我的过程的信息。一般来说，我对 GIS 还很陌生，所以我几乎可以肯定这里有一个明显的错误：
- 从 USGS GloVis 下载 Landsat 8 GeoTIFF 文件
- 将 band 5 图像重命名为更便于命令忍者使用的名称。
- 为光栅 tables 和 运行
CREATE EXTENSION postgis; 创建 postgres 数据库
运行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
这成功导入了一个新的table。然后我使用 QGIS 来直观地了解正在发生的事情。
在 Database -> DB Manager -> PostGIS -> rastertest -> public 下，我将我的 lssampleb5 添加到 canvas。
我在 QGIS 中新建了一个 XYZ 连接来添加 Google 卫星混合图像以供参考。我使用的 url 是 <a href="https://mt1.google.com/vt/lyrs=y&x=" rel="nofollow noreferrer">https://mt1.google.com/vt/lyrs=y&x=</a>{x}&y={y}&z={z} 最小和最大缩放分别为 0 和 19。
这里是我注意到 lssample 图层在 Google 混合地图上降落在西班牙海岸的地方。
我确保两层都在 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 显示的 运行ge 之外的坐标导致没有返回条目，这并不意外.
所以，首先，我不知道QGIS用的是什么坐标系。这绝对不是原始形式的经度和纬度，但根据我的理解，EPSG 4326 应该是地理投影。
其次，我不知道为什么 QGIS 会把 Landsat 场景放错地方，或者在这个过程中哪里没有正确形成场景。运行
加入我们 GIS SE，那是与 GIS 相关的地方 Q/A！
在这里为您提供帮助：
- 的确，你的罪过就是CRS。顶级
PROJCRS 标签是
键在这里，它从数据中读出 "WGS 84 / UTM zone 19N"，用
底部的 EPSG 参考 (AUTHORITY["EPSG","32619"]]).
EPSG:32619 is a UTM projected CRS based on the WGS84 geoid (datum) and with units in meter, defined as the projected distance to the corresponding reference meridians (Easting) and the equator (Northing). Since you defined the wrong CRS during import (i.e. EPSG:4326), the inherent coordinate values of the raster were treated as degree´s and the whole thing placed to the other end of the world. Run UpdateRasterSRID (SELECT UpdateRasterSRID(<shema_name>, <your_raster_table>, rast, 32619);) 将栅格的元数据设置为正确的 CRS 并重新加载图层。
- 至于 QGIS：它使用您告诉它使用的 CRS。 QGIS 带有一个非常方便的 on-the-fly 重投影功能（OTF，查看 'working with projections' here) 允许您为要投影和显示的数据定义任意 CRS（即将数据的 CRS 重新投影到内存中定义的 CRS，数据的元数据保持不变）。
您可以快速找到link 按钮到 GUI 右下角的 OTF 设置；将其设置为您想要的 SRID（例如 4326）（您会注意到数据的视觉表示如何根据所选投影发生变化。显示的坐标也将使用 CRS 单位，例如 WGS84 的十进制度数）。
我正在从事的项目要求我检索特定地理 (lon/lat) 位置的 Landsat 栅格数据。在筛选了一些教程并尝试使用 GDAL、PostGIS 和 QGIS 之后，我成功地将 GeoTIFF Landsat 图像导入到 PostGIS 栅格 table 中，并从 table 按地理位置访问值。但是，结果中存在一些问题：
- 我不明白 QGIS 在其界面中使用的坐标系，因为它们 运行 有数十万个
- 栅格加载到西班牙海岸外的 QGIS 中，而不是原先应该加载到美国缅因州的顶部。
这里有一些关于我的过程的信息。一般来说，我对 GIS 还很陌生，所以我几乎可以肯定这里有一个明显的错误：
- 从 USGS GloVis 下载 Landsat 8 GeoTIFF 文件
- 将 band 5 图像重命名为更便于命令忍者使用的名称。
- 为光栅 tables 和 运行
CREATE EXTENSION postgis; 创建 postgres 数据库
运行
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这成功导入了一个新的table。然后我使用 QGIS 来直观地了解正在发生的事情。
在
Database -> DB Manager -> PostGIS -> rastertest -> public下，我将我的 lssampleb5 添加到 canvas。我在 QGIS 中新建了一个 XYZ 连接来添加 Google 卫星混合图像以供参考。我使用的 url 是
<a href="https://mt1.google.com/vt/lyrs=y&x=" rel="nofollow noreferrer">https://mt1.google.com/vt/lyrs=y&x=</a>{x}&y={y}&z={z}最小和最大缩放分别为 0 和 19。这里是我注意到 lssample 图层在 Google 混合地图上降落在西班牙海岸的地方。
我确保两层都在 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 显示的 运行ge 之外的坐标导致没有返回条目，这并不意外.
所以，首先，我不知道QGIS用的是什么坐标系。这绝对不是原始形式的经度和纬度，但根据我的理解，EPSG 4326 应该是地理投影。
其次，我不知道为什么 QGIS 会把 Landsat 场景放错地方，或者在这个过程中哪里没有正确形成场景。运行
加入我们 GIS SE，那是与 GIS 相关的地方 Q/A！
在这里为您提供帮助：
- 的确，你的罪过就是CRS。顶级
PROJCRS标签是 键在这里，它从数据中读出"WGS 84 / UTM zone 19N"，用 底部的 EPSG 参考 (AUTHORITY["EPSG","32619"]]).
EPSG:32619 is a UTM projected CRS based on the WGS84 geoid (datum) and with units in meter, defined as the projected distance to the corresponding reference meridians (Easting) and the equator (Northing). Since you defined the wrong CRS during import (i.e. EPSG:4326), the inherent coordinate values of the raster were treated as degree´s and the whole thing placed to the other end of the world. Run UpdateRasterSRID (SELECT UpdateRasterSRID(<shema_name>, <your_raster_table>, rast, 32619);) 将栅格的元数据设置为正确的 CRS 并重新加载图层。 - 至于 QGIS：它使用您告诉它使用的 CRS。 QGIS 带有一个非常方便的 on-the-fly 重投影功能（OTF，查看 'working with projections' here) 允许您为要投影和显示的数据定义任意 CRS（即将数据的 CRS 重新投影到内存中定义的 CRS，数据的元数据保持不变）。
您可以快速找到link 按钮到 GUI 右下角的 OTF 设置；将其设置为您想要的 SRID（例如 4326）（您会注意到数据的视觉表示如何根据所选投影发生变化。显示的坐标也将使用 CRS 单位，例如 WGS84 的十进制度数）。