Plot RGB satellite image in Cassini-Soldner projection on a map in Python
<p>I have a geo-referenced RGB satellite image from the MODIS instrument in geotiff format. <strong>What is the correct way to plot it on a map using cartopy and preserve the RGB colours?</strong></p>

<p>The main obstacle I guess is the projection of the image which is Cassini-Soldner:</p>

<pre><code>import numpy as np
from osgeo import gdal, osr

ds = gdal.Open('modis_201303261252_rgb.tif')

print(ds.GetGeoTransform())
(-1669791.8857914428, 250.0, 0.0, 1669792.327327792, 0.0, -250.0)

proj = ds.GetProjection()

inproj = osr.SpatialReference()
inproj.ImportFromWkt(proj)
print(inproj)

PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Cassini_Soldner"],
    PARAMETER["latitude_of_origin",72],
    PARAMETER["central_meridian",-4],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
</code></pre>

<p>I tried to follow this example <a href="https://ocefpaf.github.io/python4oceanographers/blog/2015/03/02/geotiff/" rel="nofollow">https://ocefpaf.github.io/python4oceanographers/blog/2015/03/02/geotiff/</a>
and use cartopy to define a projection from an EPSG code. So I googled for Cassini-Soldner EPSG code (9806), but cartopy's <code>ccrs.epsg()</code> doesn't recognise it.</p>

<p>I want to use <code>plt.imshow()</code> method, but I'm a bit confused what to use as a projection keyword when the axis is created and what to pass as a <code>transform=</code> argument in <code>imshow</code>.</p>

<p>(I would add this as a comment but I need 50 rep.)</p>

<p>Out of interest, how did you come by a Cassini projection MODIS RGB?</p>

<p>An alternative is to use a Cartopy pull request I've made, which allows for the time dimension to be used with GeoAxes.add_wmts(), so that you can pull MODIS RGB imagery tiles for your time from NASA GIBS straight into a Cartopy GeoAxes: <a href="https://github.com/SciTools/cartopy/pull/788" rel="nofollow">https://github.com/SciTools/cartopy/pull/788</a></p>

python matplotlib gdal imshow cartopy