The project I am working on requires that I retrieve Landsat raster data at specific geographic (lon/lat) locations. After sifting through some tutorials and experimenting with GDAL, PostGIS, and QGIS, I successfully imported a GeoTIFF Landsat image into a PostGIS raster table and accessed values by geographic location from that table. However, there were a few issues in the result:
Here's some information about my process. I am fairly new to GIS in general, so I am almost certain theres a blatant error to be found here:
CREATE EXTENSION postgis;Run gdalinfo LSSampleB5.TIF, printing the following output:
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
I interpretted this output as EPSG 4326 format (which may be my crime), so I ran the following command to import the GeoTIFF as a PostGIS raster:
raster2pgsql -s 4326 -I LSSampleB5.TIF -F -t 50x50 -d | psql -U postgres rastertest
This successfully imported a new table. I then used QGIS to get a visual intuition of what was going on.
Under Database -> DB Manager -> PostGIS -> rastertest -> public I added my lssampleb5 to the canvas.
I created a new XYZ Connection in QGIS to add Google satillite hybrid images for reference. The url I used was https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z} with min and max zoom of 0 and 19 respectively.
Here is where I took note of the fact that the lssample layer landed off the coast of Spain on the Google Hybrid map.
I made sure both layers were on EPSG 4326 projection, no change.
Not too discouraged to move on, I tried a database query to get a single pixel value. Since my sample data landed near Spain, I used QGIS to sample a valid coordinate pair near there for the query. The query was:
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);
This returned a valid row ID and an ST_VALUE of 5776. Trying coordinates outside the range displayed by QGIS resulted in no returned entries, which isn't unexpected.
So, first of all, I do not know what QGIS is using for its coordinate system. It's definitely not longitude and latitude in a raw form, but from my understanding, EPSG 4326 is supposed to be a geographic projection.
Second, I don't know why QGIS is misplacing the Landsat scene in the wrong place, or where in the process the scene was not transformed properly.
Join us at GIS SE, that´s the place for GIS related Q/A!
To help you out here:
PROJCRS tag is the
key here, it reads out "WGS 84 / UTM zone 19N" from the data, with
the EPSG reference at the bottom (AUTHORITY["EPSG","32619"]]).SELECT UpdateRasterSRID(<shema_name>, <your_raster_table>, rast, 32619);) to set the raster's metadata to the correct CRS and reload the layer.