Search code examples
pythonmatplotlibcartopy

Saving cartopy map images for use in future plots


I am using cartopy to plot an animated trajectory on top of a GoogleTiles satellite image. I am generating thousands of frames independently and then ffmpeg-ing them together, and it seems like it would be more efficient to save a map image and to load in and use for future frames rather than having to reload the map from Google each time.

At the moment the code looks like this:

import cartopy.io.img_tiles as cimgt

google = cimgt.GoogleTiles(style='satellite')

ax = plt.axes(projection=google.crs)
ax.set_extent([west, east, south, north])
ax.add_image(google, 13)
ax.plot(data[t])

plt.savefig(fname = f'frame{t}.png')

and this code is run in a batch with t representing each individual frame. Is there a way to save and then load in a map image to avoid calling cimgt.GoogleTiles() each time?

EDIT:

I have applied the solutionbelow and made some progress. Now I have:

import cartopy.io.img_tiles as cimgt

mapbox = cimgt.MapboxTiles(access_token=token, map_id = 'satellite-v9', cache=True)

ax = plt.axes(projection=mapbox.crs)
ax.set_extent(extent)
ax.add_image(mapbox, 14)

im_arr = ax.get_children()[0].get_array()

plt.imshow(img_arr, extent=extent, origin='lower')
plt.tripcolor(triangulation, data)
plt.xlim(west,east)
plt.ylim(south,north)
plt.show()

enter image description here

You can see the issue with the projection here, the saved image looks distorted. Is there a step I missed that can display the image in the correct projection?


Solution

  • First, my comment on the relevant lines of code:-

      google = cimgt.GoogleTiles(style='satellite')
      ax = plt.axes(projection=google.crs)
    
    1. This will request several Google tiled images every time it runs for a plot.
    2. The tiled images (always having square shape) altogether will cover an extent larger than the specified plotting extent.
    3. The plot process will clip the tiled-images to the specified plotting extent.
    4. The process that always use the above code is not optimized for plotting an array of many plots that use a common base map.

    My answer will attempt to optimize the above process.

    1. The Google tile images will be requested once and cached (for next use)
    2. The tiled images will be clipped to the plot extent in the first attempt, then ...
    3. The clipped images data is extracted along with other relevant info (colormap)
    4. The obtained image data array and its colormap are the essence that we need as basemap for plotting an array of maps

    Following is the demo code that works.

    import matplotlib.pyplot as plt
    import cartopy.crs as ccrs
    import cartopy.io.img_tiles as cimgt
    import matplotlib.cm as cm
    import numpy as np
    import matplotlib
    import cartopy
    import matplotlib.path as mpath
    
    # For demo purposes
    extentd = [99.65, 101.25, 12.5, 13.75] #west,east,south,north
    zoomlev = 9
    noproj = ccrs.PlateCarree()
    #geodproj = ccrs.Geodetic()
    mercator = ccrs.Mercator()
    
    # Setup style of map
    google = cimgt.GoogleTiles(style='satellite', cache=True)
    # Setup proper projection that can do clipping
    #  to meet the requirement
    # There may be some warning at this step
    
    # Create a plot axes of choice
    ax = plt.axes(projection=mercator)     
    ax.set_extent(extentd, crs=noproj)
    # Request and process the images
    imgfactory = ax.add_image(google, zoomlev)
    
    # Check plot of the extent 
    ax.plot([extentd[0], extentd[1]], [extentd[2], extentd[3]], 'or-', 
            transform=ccrs.PlateCarree(), markersize=12, zorder=21)
    
    plt.title('Clipped Images from Google Tiles')
    plt.show()
    

    step1

    Once the first plot is done, we can grab the hidden information that are useful in the real plotting.

    # Get 'child0' object that contains the image data array
    for ix,child in enumerate(ax.get_children()):
        #print(child, type(child))
        if isinstance(child, matplotlib.image.AxesImage):
            print("* Found image object.")
            child0 = child
            ccmap = child0.get_cmap()
            img_array = child0.get_array()
            break
        pass
    

    At this step, we have

    • ccmap, the colormap of the data array
    • img_array, the data array representing the clipped images.

    Both will be used repetitively in several plots of maps.

    They can be saved to files for future uses.

    Next, we use the data array for check plot.

    # Check plot of the image data array
    # Note the use of 'img_array' and 'ccmap'
    # This uses basic matplotlib axes
    #  and use degrees as basic units
    ax2 = plt.subplot(111)
    ax2.imshow(img_array, extent=extentd , origin='lower', cmap=ccmap)
    
    plt.title('Image_data array aligned to lat-long extent')
    plt.show()
    

    step2

    And here the last step, the final array of plots that requires a common basemap is obtained.

    # Plot 3x2 array of maps
    rows = 3
    cols = 2
    
    # Prep (x,y) data for extra features on each plot
    xss = np.linspace(99.65, 101.25, 36)
    # yss = 13+np.random.random(36)/2
    
    # This uses PlateCarree projection, other projection can be problematic
    fig, ax = plt.subplots(nrows=rows, ncols=cols, figsize=(13, 15), 
                           subplot_kw={'projection': ccrs.PlateCarree()})
    
    # This uses Mercator 
    # Warning, image array takes v. long time to do image warping
    #fig, ax = plt.subplots(nrows=rows, ncols=cols, figsize=(5.2, 6), 
    #                       subplot_kw={'projection': ccrs.Mercator()})
    
    for i in range(rows):
        for j in range(cols):
            # Plot base raster map, using data array
            ax[i][j].imshow(child0.get_array(), origin='lower', 
                            cmap=child0.get_cmap(), extent=extentd, zorder=10)
    
            # Plot individual data sets on subplots
            ax[i][j].plot(xss, 13+np.random.random(len(xss))/2, 
                          color='cyan', zorder=20, 
                          transform=ccrs.PlateCarree(), linewidth=2)
            # Write map title
            ax[i][j].set_title("Row:"+str(i)+", Col:"+str(j))
            ax[i][j].set_extent(extentd)
            gl=ax[i][j].gridlines(crs=ccrs.PlateCarree(), draw_labels=[1,1,0,0], 
                               lw=1, color='white', linestyle='--', zorder=20)
            gl.top_labels = False
            gl.right_labels = False
        pass
    
    plt.tight_layout()
    plt.show()
    

    step3