Cartopy

Cartopy is a Python package designed for geospatial data processing in order to produce maps and other geospatial data analyses.

Cartopy makes use of the powerful PROJ.4, NumPy and Shapely libraries and includes a programmatic interface built on top of Matplotlib for the creation of publication quality maps.

Key features of cartopy are its object oriented projection definitions, and its ability to transform points, lines, vectors, polygons and images between those projections.

You will find cartopy especially useful for large area / small scale data, where Cartesian assumptions of spherical data traditionally break down. If you’ve ever experienced a singularity at the pole or a cut-off at the dateline, it is likely you will appreciate cartopy’s unique features!

Prerequisites

$ pip install cython matplotlib scipy pillow

Beautifully simple maps

Cartopy has exposed an interface to enable easy map creation using matplotlib. Creating a basic map is as simple as telling matplotlib to use a specific map projection, and then adding some coastlines to the axes:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
plt.show()

Offline cache example

Local cache for Open Street Map tiles #732

import os
import types

import cartopy.io.img_tiles as img_tiles
import requests
import PIL


class CachedTiler(object):
    def __init__(self, tiler):
        self.tiler = tiler

    # Python2:
    # def __getattr__(self, name):
    #     # Mimic the tiler interface, but for methods, ensure that the "self"
    #     # that is passed through continues to be CachedTiler, and not the
    #     # contained tiler instance.
    #     attr = getattr(self.tiler, name, None)
    #     if isinstance(attr, types.MethodType):
    #         attr = types.MethodType(attr.im_func, self)
    #     return attr

    def __getattr__(self, name):
        """
        Python3 version:
        """
        attr = getattr(self.tiler, name, None)
        if isinstance(attr, types.MethodType):
            attr = types.MethodType(attr.__func__, self)
        return attr

    def get_image(self, tile):
        tileset_name = '{}'.format(self.tiler.__class__.__name__.lower())
        cache_dir = os.path.expanduser(os.path.join('~/', 'image_tiles', tileset_name))
        if not os.path.exists(cache_dir):
            os.makedirs(cache_dir)
        tile_fname = os.path.join(cache_dir, '_'.join(str(v) for v in tile) + '.png')
        if not os.path.exists(tile_fname):
            response = requests.get(self._image_url(tile),
                                    stream=True)

            with open(tile_fname, "wb") as fh:
                for chunk in response:
                    fh.write(chunk)
        with open(tile_fname, 'rb') as fh:
            img = PIL.Image.open(fh)
            img = img.convert(self.desired_tile_form)
        return img, self.tileextent(tile), 'lower'



from cartopy.io.img_tiles import StamenTerrain
import matplotlib.pyplot as plt


actual_tiler = StamenTerrain()
tiler = CachedTiler(actual_tiler)

mercator = tiler.crs
ax = plt.axes(projection=mercator)
ax.set_extent([-90, -73, 22, 34])
ax.add_image(tiler, 6)
ax.coastlines('10m')
plt.show()

OSM Preview

Jupyter Notebook Viewer - Example of using cartopy.io.img_tiles.OSM for Cardiff, with 1:10m coastlines overlayed

import cartopy.crs as ccrs
from cartopy.io.img_tiles import OSM
import matplotlib.pyplot as plt

osm_tiles = OSM()

plt.figure(figsize=(16, 16))

# Use the tile's projection for the underlying map.
ax = plt.axes(projection=osm_tiles.crs)

# Specify a region of interest, in this case, Cardiff.
ax.set_extent([-3.3849, -3.0424, 51.4260, 51.5842],
              ccrs.PlateCarree())

# Add the tiles at zoom level 12.
ax.add_image(osm_tiles, 12)

ax.coastlines('10m')

plt.show()

Longitude and Latitude guidline

from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

ax = plt.axes(projection=request.crs)
gl = ax.gridlines(draw_labels=True, alpha=0.2)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER

Draw point and label

Map tile acquisition

Demonstrates cartopy’s ability to draw map tiles which are downloaded on demand from the MapQuest tile server. Internally these tiles are then combined into a single image and displayed in the cartopy GeoAxes.

import matplotlib.pyplot as plt
from matplotlib.transforms import offset_copy

import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt

# ...

# Add a marker for the Eyjafjallajökull volcano.
plt.plot(-19.613333, 63.62, marker='o', color='yellow', markersize=12, alpha=0.7, transform=ccrs.Geodetic())

# Use the cartopy interface to create a matplotlib transform object
# for the Geodetic coordinate system. We will use this along with
# matplotlib's offset_copy function to define a coordinate system which
# translates the text by 25 pixels to the left.
geodetic_transform = ccrs.Geodetic()._as_mpl_transform(ax)
text_transform = offset_copy(geodetic_transform, units='dots', x=-25)

# Add text 25 pixels to the left of the volcano.
plt.text(-19.613333, 63.62, u'Eyjafjallajökull',
         verticalalignment='center', horizontalalignment='right',
         transform=text_transform,
         bbox=dict(facecolor='wheat', alpha=0.5, boxstyle='round'))

GoogleMap Cache example

한국/오송 출력.

import os
import types
import six
import requests

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

import cartopy.io.img_tiles as img_tiles
import cartopy.crs as ccrs

from multiprocessing import Process, Lock
from PIL import Image
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER


mutex = Lock()


class CachedTiler(object):
    def __init__(self, tiler):
        self.tiler = tiler

    # Python2:
    # def __getattr__(self, name):
    #     # Mimic the tiler interface, but for methods, ensure that the "self"
    #     # that is passed through continues to be CachedTiler, and not the
    #     # contained tiler instance.
    #     attr = getattr(self.tiler, name, None)
    #     if isinstance(attr, types.MethodType):
    #         attr = types.MethodType(attr.im_func, self)
    #     return attr

    def __getattr__(self, name):
        """
        Python3 version:
        """
        attr = getattr(self.tiler, name, None)
        if isinstance(attr, types.MethodType):
            attr = types.MethodType(attr.__func__, self)
        return attr

    def get_image(self, tile):
        tileset_name = '{}'.format(self.tiler.__class__.__name__.lower())
        cache_dir = os.path.expanduser(os.path.join(os.getcwd(), 'image_tiles', tileset_name))
        if not os.path.exists(cache_dir):
            os.makedirs(cache_dir)
        tile_fname = os.path.join(cache_dir, '_'.join(str(v) for v in tile) + '.png')
        if not os.path.exists(tile_fname):
            #with mutex:
            url = self._image_url(tile)
            print('Download Image URL: {}'.format(url))
            response = requests.get(url, stream=True)
            print('Result: {}'.format(response))
            with open(tile_fname, "wb") as fh:
                for chunk in response:
                    fh.write(chunk)
        with open(tile_fname, 'rb') as fh:
            img = Image.open(fh)
            img = img.convert(self.desired_tile_form)
        return img, self.tileextent(tile), 'lower'


def new_get_image(self, tile):
    if six.PY3:
        from urllib.request import urlopen, Request
    else:
        from urllib2 import urlopen
        url = self._image_url(tile)  # added by H.C. Winsemius
        req = Request(url) # added by H.C. Winsemius
        # req.add_header('User-agent', 'your bot 0.1')
        # fh = urlopen(url)  # removed by H.C. Winsemius
        fh = urlopen(req)
        im_data = six.BytesIO(fh.read())
        fh.close()
        img = Image.open(im_data)
        img = img.convert(self.desired_tile_form)
        return img, self.tileextent(tile), 'lower'


#img_tiles.GoogleWTS.get_image = new_get_image
#request = img_tiles.GoogleTiles()
#extent = [-89, -88, 41, 42]
#ax = plt.axes(projection=request.crs)
#ax.set_extent(extent)
#ax.add_image(request, 8)
#plt.show()

import cartopy.crs as ccrs
from cartopy.io.img_tiles import OSM
import matplotlib.pyplot as plt

OSM.get_image = new_get_image

#actual_tiler = OSM()
actual_tiler = img_tiles.GoogleTiles()
osm_tiles = CachedTiler(actual_tiler)

plt.figure(figsize=(16, 16))
# fig = plt.figure(figsize=(9, 13))
# Use the tile's projection for the underlying map.
ax = plt.axes(projection=osm_tiles.crs)

gl = ax.gridlines(draw_labels=True, alpha=0.2)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER

#extent = [-39, -38.25, -13.25, -12.5]
#extent = [36.6412, 127.2951, 36.6258, 127.3273]
extent = [127.2951, 127.3273, 36.6412, 36.6258]
#extent = [36.38, 127.17, 36.37, 127.19]

# Specify a region of interest, in this case, Cardiff.
# ax.set_extent([-3.3849, -3.0424, 51.4260, 51.5842], ccrs.PlateCarree())
# ax.set_extent([36.6412, 127.2951, 36.6258, 127.3273], ccrs.PlateCarree())
#ax.set_extent(extent, osm_tiles.crs)
ax.set_extent(extent, ccrs.PlateCarree())
#print(ax.get_extent())
#ax.set_xlim([extent[0], extent[2]])
#ax.set_ylim([extent[1], extent[3]])

# Add the tiles at zoom level 12.
ax.add_image(osm_tiles, 16)
ax.coastlines('10m')
plt.show()

Troubleshooting

Unable to determine GEOS version. Ensure you have 3.3.3 or later installed, or installation may fail

Stackoverflow - Using pip install to install Cartopy but missing Proj version at least 4.9.0

pip로 설치시 아래와 같은 에러가 발생할 수 있다.

$ pip install cartopy
Collecting cartopy
  Downloading https://files.pythonhosted.org/packages/e5/92/fe8838fa8158931906dfc4f16c5c1436b3dd2daf83592645b179581403ad/Cartopy-0.17.0.tar.gz (8.9MB)
     |████████████████████████████████| 8.9MB 848kB/s
  Installing build dependencies ... done
  Getting requirements to build wheel ... error
  ERROR: Complete output from command /Users/your/.pyenv/versions/3.6.8/envs/opy-your-3.6.8/bin/python /Users/your/.pyenv/versions/3.6.8/envs/opy-your-3.6.8/lib/python3.6/site-packages/pip/_vendor/pep517/_in_process.py get_requires_for_build_wheel /var/folders/6v/598t63jj15lbcn_nq9nm12ww0000gn/T/tmprp39fedd:
  ERROR: setup.py:171: UserWarning: Unable to determine GEOS version. Ensure you have 3.3.3 or later installed, or installation may fail.
    '.'.join(str(v) for v in GEOS_MIN_VERSION), ))
  Proj 4.9.0 must be installed.
  ----------------------------------------
ERROR: Command "/Users/your/.pyenv/versions/3.6.8/envs/opy-your-3.6.8/bin/python /Users/your/.pyenv/versions/3.6.8/envs/opy-your-3.6.8/lib/python3.6/site-packages/pip/_vendor/pep517/_in_process.py get_requires_for_build_wheel /var/folders/6v/598t63jj15lbcn_nq9nm12ww0000gn/T/tmprp39fedd" failed with error code 1 in /private/var/folders/6v/598t63jj15lbcn_nq9nm12ww0000gn/T/pip-install-y2vdys54/cartopy

geos와 proj를 설치해야 한다.

Ubuntu는 아래와 같이 설치한다.

$ apt-get install libproj-dev proj-data proj-bin
$ apt-get install libgeos-dev
$ pip install cython
$ pip install cartopy

MacOS는 brew를 사용하면 된다.

$ brew install geos
$ brew install proj

ValueError: A non-empty list of tiles should be provided to merge

.io.img_tiles returns "ValueError: A non-empty list of tiles should be provided to merge." #1341

import numpy as np 
import matplotlib as mpl        
import matplotlib.pyplot as plt 
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt
import six
from PIL import Image

def new_get_image(self, tile):
    if six.PY3:
        from urllib.request import urlopen, Request
    else:
        from urllib2 import urlopen
    url = self._image_url(tile)  # added by H.C. Winsemius
    req = Request(url) # added by H.C. Winsemius
    req.add_header('User-agent', 'your bot 0.1')
    # fh = urlopen(url)  # removed by H.C. Winsemius
    fh = urlopen(req)
    im_data = six.BytesIO(fh.read())
    fh.close()
    img = Image.open(im_data)
    img = img.convert(self.desired_tile_form)

    return img, self.tileextent(tile), 'lower'

cimgt.GoogleWTS.get_image = new_get_image
request = cimgt.GoogleTiles()
extent = [-89, -88, 41, 42]

ax = plt.axes(projection=request.crs)
ax.set_extent(extent)

ax.add_image(request, 8)
plt.show()

OSM Example:

from cartopy.io.img_tiles import OSM
OSM.get_image = new_get_image
osm_tiles = OSM()

GeoAxes.clear() 이후 그림이 그려지지 않는 현상

GeoAxes.clear()이후, 이미지를 다시 그릴 경우 cartopy.mpl.GeoAxes._done_img_factory변수가 True로 설정되면서 cartopy.mpl.GeoAxes.draw() 함수에서 self.imshow()호출부가 SKIP 된다. 따라서 해당 플래그를 False로 강제 변경하면 된다.

class PlotCanvasPanel(wx.Panel):
    def __init__(self, *args, **kwds):
        wx.Panel.__init__(self, *args, **kwds)

        self.figure = Figure()
        self.figure.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0, wspace=0.0, hspace=0.0)
        self.axes = self.figure.add_subplot(111, projection=GoogleTiles().crs)
        # ...
        self.canvas = FigureCanvas(self, -1, self.figure)
        self.sizer = wx.BoxSizer(wx.VERTICAL)
        self.sizer.Add(self.canvas, proportion=1, border=0, flag=wx.ALL | wx.GROW)

        self.SetSizer(self.sizer)
        self.Fit()

    def draw(self):
        self.axes: GeoAxes
        self.axes.clear()  # CLEAR !!
        self.axes._done_img_factory = False  # Force flag setting.

        self.axes.add_image(CachedTiler(GoogleTiles()), self.zoom_level)
        self.axes.set_extent(self.extent, ccrs.PlateCarree())
        self.axes.coastlines('10m')
        # ...
        self.canvas.draw()