GOES-16 geolocation seems off when saving as TIFF

[simonrp84]

Describe the bug
Hello,
One of my students has noticed that the geolocation seems to be incorrect for GOES-16. A GeoTIFF saved by SatPy has an offset compared to a coastline overlay in QGIS.

To Reproduce
The following python script will produce a GeoTIFF with the incorrect location info

#GOES-16
from      datetime          import datetime,timedelta
from      satpy             import Scene
from      glob             import glob
import   os

BASE_DIR      =  {input_dir}
OUT_DIR      =  {output_dir}

count         =   0
fnames   =   glob(os.path.join(BASE_DIR, 'OR*-Rad*-*C01*_s2018350*.nc'))
fnames.sort()

print("Number of Scenes: ", len(fnames))

incount         =   1

for f in fnames:

	infs   =   glob(f.replace('C01', 'C*')[:len(BASE_DIR) + 50] + '*.nc')
	scn   =   Scene(reader='abi_l1b', filenames=infs)

	print('   -   Processing',incount,'of',len(fnames))
	scn.load(['true_color'])
	new	=	scn.resample(scn.min_area(), resampler='native')
	outf   =   OUT_DIR+'GOES_'+str(incount).zfill(4)+'_RGB.tiff'
	new.save_dataset('true_color', outf)
	incount      +=   1

Then load the image in QGIS (or similar) and add a coastline overlay. Some areas, particularly California, South-Eastern South America and the Atlantic Islands, have a 10+pixel shift in geolocation.

Expected behavior
That a coastline overlaid in GIS software will match the actual coastlines visible in the L1 sat data.

Actual results
This is the gdalinfo output for a GeoTIFF saved by SatPy:
PROJCS["unnamed",
GEOGCS["unnamed ellipse",
DATUM["unknown",
SPHEROID["unnamed",6378137,298.2572220960422]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Geostationary_Satellite"],
PARAMETER["central_meridian",-75],
PARAMETER["satellite_height",35786023],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["Meter",1],
EXTENSION["PROJ4","+a=6378137.0 +b=6356752.31414 +lon_0=-75.0 +h=35786023.0 +proj=geos +units=m +sweep=x"]]

Screenshots

Please see these two images for an example:
https://www2.physics.ox.ac.uk/sites/default/files/profiles/proud/bob1-44130.png
https://www2.physics.ox.ac.uk/sites/default/files/profiles/proud/bob2-44131.png

Environment Info:
Ubuntu 14.04 + Mint 19 + Ubuntu 18.04
Latest SatPy pulled from git.

[djhoese]

As a first guess, I have a feeling this is due to QGIS not supporting the geos projection's +sweep parameter. It isn't QGIS's fault, the GeoTIFF format itself doesn't officially support the parameter. Most people have been getting around this issue in two ways:

1. Resample the data again using gdalwarp to a projection with the default sweep parameter.
2. Using newer versions of gdal to both write (create the geotiffs) and read the data (QGIS). GDAL has added a "hack" to include the sweep information so that it can understand it.

The gdalinfo shows that the version of gdal used to create the images is new enough that GDAL is "hacking" in the sweep parameter in the "EXTENSION" part. I'm guessing QGIS's libproj or gdal version is older though. To prove that this is the problem, I'm guessing that the distortion/shift is not is bad in the upper-left and lower-right portions of the image compared to the Baja California and eastern Canada regions.

I can give the gdalwarp command to use if you want to try method 1 above.

[simonrp84]

Hah, that's simple! Updating QGIS to the latest build (3.4) fixed the problem, coastlines now appear in the right places. Apologies for blaming SatPy for this issue.

[djhoese]

Great to hear it was that simple and good to know for future my reference. Thanks for filing the bug.