Retrieving raster data by geographic location using Landsat and PostGIS

Question

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:

I do not understand the coordinate system being used by QGIS in its interface, as they range in the hundred thousands
The raster loaded into QGIS off the coast of Spain, rather than on top of Maine, USA as it was supposed to.

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:

Download Landsat 8 GeoTIFF file from USGS GloVis
Rename the band 5 image to something more friendly to command ninja with.
Create postgres database for raster tables and run 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.

geozelot geozelot · Accepted Answer · 2018-03-19T10:14:16

Join us at GIS SE, that´s the place for GIS related Q/A!

To help you out here:

Indeed, your crime was the CRS. The top level 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"]]).
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);) to set the raster's metadata to the correct CRS and reload the layer.
As for QGIS: it uses the CRS that you tell it to use. QGIS comes with a very handy on-the-fly reprojection feature (OTF, check out the general man page on 'working with projections' here) that lets you define an arbitrary CRS for your data to be projected and displayed (i.e. it reprojects the data's CRS into the defined one in memory, the data's metadata stays untouched).
You can find the quick link button to the OTF settings in the bottom right corner of the GUI; set it to your desired SRID (e.g. 4326) (you´ll notice how the visual representation of your data changes according to the chosen projection. also the displayed coordinates will use the CRS units, e.g. decimal degrees for WGS84).

Retrieving raster data by geographic location using Landsat and PostGIS

1 Answers