Reusing arcgisimage AxesImage instance

[Gormador]

Hi!

First of all, thank you for your work on basemap!

I'm looking for a way of reusing the AxesImage instance returned by the arcgisimage function of basemap, in order to reapply it to following and identical maps on which only the plotted data differs.
I tried using the matplotlib.axes.Axes.add_image() function on a basemap instance, thinking it might do the trick but alas, it does not.

I may simply be missing something rather obvious that would fulfill my need, but I have no clue of what it might be.

Here's the function that generates basemap instances and contains my (failed) attempt described above, as well as the straightforward -- but very unoptimized -- fix:

def getBasemapsInstances(nbGraphs, axes, swCornerCoord, neCornerCoord, arcgisService):

    mAxes = []

    swCorLat = swCornerCoord[0]
    swCorLon = swCornerCoord[1]
    neCorLat = neCornerCoord[0]
    neCorLon = neCornerCoord[1]

    la0 = swCorLat + (neCorLat - swCorLat)/2
    lo0 = swCorLon + (neCorLon - swCorLon)/2

    for i in range(1, nbGraphs):

        bmap = Basemap(ax=axes[i], 
                        llcrnrlon=swCorLon,llcrnrlat=swCorLat,
                        urcrnrlon=neCorLon,urcrnrlat=neCorLat,
                        epsg=3035,lon_0=lo0,lat_0=la0)

        ''' Doesn't work
        if i==1:
            arcServInst = bmap.arcgisimage(service=arcgisService,
                                            xpixels = 2000,
                                            dpi=1500,
                                            verbose= True)
        else:
            bmap.add_image(arcServInst)
        '''

        arcServInst = bmap.arcgisimage(service=arcgisService,
                                        xpixels = 2000,
                                        dpi=1500,
                                        verbose= True)

        bmap.drawmapboundary()

        mAxes.append(bmap)

    return mAxes

Thanks in advance for any input you may have! :-)

[efiring]

What you need to do is look at the code for the arcgisimage method and extract the parts up to but not including the actual imshow. You would return the result of the imread. Then call your new function outside your loop, save its imread output, and supply it to a call to bmap.imshow inside your loop.

(Side note: it looks like you are accustomed to Matlab-style indexing; better to just switch consistently to Python-style. As it is, you are ignoring your first axes element and processing nbGraphs-1 elements.)

[Gormador]

OK, will do. I was hoping to avoid doing this but it's not the end of the world. Also, will see how the function works firsthand!

As for the indexing -- and you can't guess it from this bit of code -- I actually am using the full range of my axes array. The element at index 0 is just not a basemap instance and is handled somewhere else. So range(1, nbGraphs) is absolutely intentional :-) Also, never used Matlab in my life! The first line of code I wrote was in C, which doesn't forgive such silliness ;-)

If you are curious, axes[0] contains a colormap that serves as a reference for what I'm plotting.

Anyway, thanks a lot for your answer and thank you for your contributions to Matplotlib :-)

[Gormador]

A bit of a follow-up: I took the time to do as suggested by @efiring in his answer and, short of a few tweaks here and there, things work perfectly as expected.

This is potentially a huge time saver when trying to display several maps with the same background image... which is exactly what I'm doing. So thanks again!

Below, my modifications to the code to make it work outside the library.
Beware: this is intended for python 3.5. Conflicting bits concern the urllib import.

def getArcgisimage(bmapInstance, server='http://server.arcgisonline.com/ArcGIS',
                        service='ESRI_Imagery_World_2D', 
                        xpixels=400, ypixels=None, dpi=96, 
                        verbose=False, **kwargs):
    import pyproj
    import matplotlib.image as mlabIm
    from urllib.request import urlopen

    _cylproj = ['cyl','merc','mill','gall','cea']

    if not hasattr(bmapInstance,'epsg'):
        msg = dedent("""
        Basemap instance must be creating using an EPSG code
        (http://spatialreference.org) in order to use the wmsmap method""")
        raise ValueError(msg)
    # find the x,y values at the corner points.
    p = pyproj.Proj(init="epsg:%s" % bmapInstance.epsg, preserve_units=True)
    xmin,ymin = p(bmapInstance.llcrnrlon,bmapInstance.llcrnrlat)
    xmax,ymax = p(bmapInstance.urcrnrlon,bmapInstance.urcrnrlat)
    if bmapInstance.projection in _cylproj:
        Dateline =\
        _geoslib.Point(
                    bmapInstance(180.,0.5*(bmapInstance.llcrnrlat+bmapInstance.urcrnrlat))
                )
        hasDateline = Dateline.within(bmapInstance._boundarypolyxy)
        if hasDateline:
            msg=dedent("""
            arcgisimage cannot handle images that cross
            the dateline for cylindrical projections.""")
            raise ValueError(msg)
    if bmapInstance.projection == 'cyl':
        xmin = (180./np.pi)*xmin; xmax = (180./np.pi)*xmax
        ymin = (180./np.pi)*ymin; ymax = (180./np.pi)*ymax
    # ypixels not given, find by scaling xpixels by the map aspect ratio.
    if ypixels is None:
        ypixels = int(bmapInstance.aspect*xpixels)
    # construct a URL using the ArcGIS Server REST API.
    basemap_url = \
"%s/rest/services/%s/MapServer/export?\
bbox=%s,%s,%s,%s&\
bboxSR=%s&\
imageSR=%s&\
size=%s,%s&\
dpi=%s&\
format=png32&\
f=image" %\
(server,service,xmin,ymin,xmax,ymax,bmapInstance.epsg,bmapInstance.epsg,xpixels,ypixels,dpi)
    # print URL?
    if verbose: print(basemap_url)
    return mlabIm.imread(urlopen(basemap_url))

And the adapted function from OP:

def getBasemapsInstances(nbGraphs, axes, swCornerCoord, neCornerCoord, arcgisService):

    mAxes = []

    swCorLat = swCornerCoord[0]
    swCorLon = swCornerCoord[1]
    neCorLat = neCornerCoord[0]
    neCorLon = neCornerCoord[1]

    la0 = swCorLat + (neCorLat - swCorLat)/2
    lo0 = swCorLon + (neCorLon - swCorLon)/2

    argisimage = None

    for i in range(1, nbGraphs+1):

        bmap = Basemap(ax=axes[i], 
                        llcrnrlon=swCorLon,llcrnrlat=swCorLat,
                        urcrnrlon=neCorLon,urcrnrlat=neCorLat,
                        epsg=3035,lon_0=lo0,lat_0=la0)

        if i==1:
            arcgisimage = getArcgisimage(bmap, service=arcgisService,
                                            xpixels = 2000,
                                            dpi=1500,
                                            verbose= True)

        bmap.imshow(arcgisimage, origin='upper')

        bmap.drawmapboundary()

        mAxes.append(bmap)

    return mAxes