如何使用cartopy和matplotlib创建比例尺?

问题描述 投票:0回答:2

相对于 stackoverflow 中之前的示例,我搜索了其他替代方案以创建比例尺。

在我的研究中,我验证了 mpl_toolkits.basemap 中的 Basemap 类请参阅此处。它有“drawmapscale”方法。此方法具有选项 barstyle = 'fancy' 以获得更有趣的比例尺绘图。

因此,我尝试将底图的“drawmapscale”转换为cartopy版本。

尽管如此,结果并不乐观,我从图中收到了错误消息。我认为错误在于数据的转换。

这是脚本:


import numpy as np

import matplotlib.pyplot as plt

import pyproj

import cartopy.crs as ccrs
from matplotlib import is_interactive
from cartopy.crs import (WGS84_SEMIMAJOR_AXIS, WGS84_SEMIMINOR_AXIS)
from pyproj import Transformer

class Scalebar():
    
    
    def __init__(self, 
                 ax,
                 suppress_ticks=True,
                 geographical_crs = 4326,
                  planar_crs = 5880,
                  
                 fix_aspect=True,
                 anchor='C',
                  celestial=False,
                  round=False,
                  noticks=False,
                  metric_ccrs=ccrs.TransverseMercator()):
        
        self.ax = ax
        self.fix_aspect = fix_aspect
        
        # setting metric ccrs for reference in the plotting
        
        self.metric_ccrs = metric_ccrs
        

        self.anchor = anchor
        # geographic or celestial coords?
        self.celestial = celestial
        # map projection.
        self.projection = ax.projection
        
        self.geographical_crs = geographical_crs
        self.planar_crs = planar_crs
        

        self._initialized_axes = set()
        

        self.round = round
        
        # map boundary not yet drawn.
        self._mapboundarydrawn = False
        
        self.rmajor = np.float(ax.projection.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS)
        self.rminor = np.float(ax.projection.globe.semiminor_axis or WGS84_SEMIMINOR_AXIS)
        
        # set instance variables defining map region.
        
        self.xmin = self.projection.boundary.bounds[0]
        self.xmax = self.projection.boundary.bounds[2]
        self.ymin = self.projection.boundary.bounds[1]
        self.ymax = self.projection.boundary.bounds[3]
        
        
        
        self._width = self.xmax - self.xmin
        self._height = self.ymax - self.ymin
        
        self.noticks = noticks
        
        
    def __call__(self,x,y,
                 
                 inverse=False
                 ):
        """
        Calling the class instance with the arguments lon, lat will
        convert lon/lat (in degrees) to x/y map projection
        coordinates (in meters).
        
        If optional keyword ``inverse`` is True (default is False), 
        the inverse transformation from x/y to lon/lat is performed.
        
        Input arguments:
            lon, lat can be either scalar floats, sequences, or numpy arrays.
        """
        
    
        if not inverse:
                
            transformer = Transformer.from_crs("epsg:{0}".format(self.geographical_crs), 
                                               "epsg:{0}".format(self.planar_crs))
            
        else:
            transformer = Transformer.from_crs("epsg:{0}".format(self.planar_crs), 
                                               "epsg:{0}".format(self.geographical_crs))
        
        
        return transformer.transform(x, y)

    
    def drawmapscale(self, 
                     lon,
                     lat,
                     length,
                     lon0=None,
                     lat0=None,
                     barstyle='simple',\
                     units='km',
                     fontsize=9,
                     yoffset=None,
                     labelstyle='simple',\
                     fontcolor='k',
                     fillcolor1='w',
                     fillcolor2='k',\
                     format='%d',
                     zorder=None,
                     linecolor=None,
                     linewidth=None):
        """
        Draw a map scale at ``lon,lat`` of length ``length``
        representing distance in the map
        projection coordinates at ``lon0,lat0``.
        .. tabularcolumns:: |l|L|
        ==============   ====================================================
        Keywords         Description
        ==============   ====================================================
        units            the units of the length argument (Default km).
        barstyle         ``simple`` or ``fancy`` (roughly corresponding
                         to the styles provided by Generic Mapping Tools).
                         Default ``simple``.
        fontsize         for map scale annotations, default 9.
        fontcolor            for map scale annotations, default black.
        labelstyle       ``simple`` (default) or ``fancy``.  For
                         ``fancy`` the map scale factor (ratio betwee
                         the actual distance and map projection distance
                         at lon0,lat0) and the value of lon0,lat0 are also
                         displayed on the top of the scale bar. For
                         ``simple``, just the units are display on top
                         and the distance below the scale bar.
                         If equal to False, plot an empty label.
        format           a string formatter to format numeric values
        yoffset          yoffset controls how tall the scale bar is,
                         and how far the annotations are offset from the
                         scale bar.  Default is 0.02 times the height of
                         the map (0.02*(self.ymax-self.ymin)).
        fillcolor1(2)    colors of the alternating filled regions
                         (default white and black).  Only relevant for
                         'fancy' barstyle.
        zorder           sets the zorder for the map scale.
        linecolor        sets the color of the scale, by default, fontcolor
                         is used
        linewidth        linewidth for scale and ticks
        ==============   ====================================================
        Extra keyword ``ax`` can be used to override the default axis instance.
        """
        
        
        # get current axes instance (if none specified).
        ax = self.ax

        # convert length to meters
        lenlab = length
        if units == 'km':
            length = length*1000
        elif units == 'mi':
            length = length*1609.344
        elif units == 'nmi':
            length = length*1852
        elif units == 'ft':
            length = length*0.3048
        elif units != 'm':
            msg = "units must be 'm' (meters), 'km' (kilometers), "\
            "'mi' (miles), 'nmi' (nautical miles), or 'ft' (feet)"
            raise KeyError(msg)
        
        # Setting the center coordinates of the axes:
        
        xmin, xmax, ymin, ymax = self.ax.get_extent()
        
        if lon0 == None:
        
            lon0 = np.mean([xmin, xmax])
            
        if lat0 == None:
            lat0 = np.mean([ymin, ymax])
    
        
        # reference point and center of scale.
            
        x0,y0 = self(lon0,lat0)
        
        
        print('\n\n Central coords prior to transform')
        print('lon0,lat0: ', [lon0,lat0])
        
        print('\n\n central coordinates after transform')
        print('x0,y0: ', [x0,y0])
        
        
        
        xc,yc = self(lon,lat)
        
        
        
        print('\n\n positional coordinates prior to transform')
        print('lon, lat: ', [lon,lat])
        
        print('\n\n central coordinates after transform')
        print('xc,yc: ', [xc,yc])
        
        print('-'*20, '\n')
        
        # make sure lon_0 between -180 and 180
        lon_0 = ((lon0+360) % 360) - 360
        if lat0>0:
            if lon>0:
                lonlatstr = u'%g\N{DEGREE SIGN}N, %g\N{DEGREE SIGN}E' % (lat0,lon_0)
            elif lon<0:
                lonlatstr = u'%g\N{DEGREE SIGN}N, %g\N{DEGREE SIGN}W' % (lat0,lon_0)
            else:
                lonlatstr = u'%g\N{DEGREE SIGN}, %g\N{DEGREE SIGN}W' % (lat0,lon_0)
        else:
            if lon>0:
                lonlatstr = u'%g\N{DEGREE SIGN}S, %g\N{DEGREE SIGN}E' % (lat0,lon_0)
            elif lon<0:
                lonlatstr = u'%g\N{DEGREE SIGN}S, %g\N{DEGREE SIGN}W' % (lat0,lon_0)
            else:
                lonlatstr = u'%g\N{DEGREE SIGN}S, %g\N{DEGREE SIGN}' % (lat0,lon_0)
        
        # left edge of scale
         
        lon1,lat1 = self(x0-length/2,y0, inverse=True)
        x1,y1 = self(lon1,lat1)
        # right edge of scale
        lon4,lat4 = self(x0+length/2,y0, inverse=True)
        x4,y4 = self(lon4,lat4)
        
        x1 = x1-x0+xc
        y1 = y1-y0+yc
        
        
        
        
        print('\n\n positional coordinates prior to transform')
        print('lon1,lat1: ', [lon1,lat1])
        
       
        
        print('\n\n positional coordinates prior to transform')
        print('x1, y1: ', [x1,y1])
        
        print() 
        
        print('\n\n central coordinates after transform')
        print('lon4,lat4: ', [lon4,lat4])
        
        print('-'*20, '\n')
        
        
        print('\n\n central coordinates after transform')
        print('x4,y4: ', [x4,y4])
        
        print('-'*20, '\n')
        
        
        x4 = x4-x0+xc
        y4 = y4-y0+yc
        if x1 > 1.e20 or x4 > 1.e20 or y1 > 1.e20 or y4 > 1.e20:
            raise ValueError("scale bar positioned outside projection limb")
        # scale factor for true distance
            
            
        gc = pyproj.Geod(a=self.rmajor,b=self.rminor)
        az12,az21,dist = gc.inv(lon1,lat1,lon4,lat4)
        scalefact = dist/length
        # label to put on top of scale bar.
        
        if labelstyle=='simple':
            labelstr = units
        
        elif labelstyle == 'fancy':
            labelstr = units+" (scale factor %4.2f at %s)"%(scalefact,lonlatstr)
        elif labelstyle == False:
            labelstr = ''
        else:
            raise KeyError("labelstyle must be 'simple' or 'fancy'")
        
        # default y offset is 2 percent of map height.
        if yoffset is None: 
            yoffset = 0.02*(self.ymax-self.ymin)
        
        rets = [] # will hold all plot objects generated.
        
        # set linecolor
        if linecolor is None:
            linecolor = fontcolor
        
        # 'fancy' style
        if barstyle == 'fancy':
            #we need 5 sets of x coordinates (in map units)
            #quarter scale
            lon2,lat2 = self(x0-length/4,y0,inverse=True)
            x2,y2 = self(lon2,lat2)
            x2 = x2-x0+xc; y2 = y2-y0+yc
            #three quarter scale
            lon3,lat3 = self(x0+length/4,y0,inverse=True)
            x3,y3 = self(lon3,lat3)
            x3 = x3-x0+xc; y3 = y3-y0+yc
            #plot top line
            ytop = yc+yoffset/2
            ybottom = yc-yoffset/2
            ytick = ybottom - yoffset/2
            ytext = ytick - yoffset/2
            
            lontext , lattext = self(lon0,ytext, inverse=True)
            
            #lon_top, lat_top = self(lon4,ytop,inverse=True)
            #lon_top, lat_bottom = self(lon4,ybottom,inverse=True)
            
            transform = self.metric_ccrs # this crs projection is meant to be for metric data
            
            rets.append(self.plot([x1,x4],
                                  [ytop,ytop],
                                  transform=transform,
                                  color=linecolor, 
                                  linewidth=linewidth)[0])
            
            #plot bottom line
            rets.append(self.plot([x1,x4],
                                  [ybottom,ybottom],
                                  transform=transform,
                                  color=linecolor, 
                                  linewidth=linewidth)[0])
            
            #plot left edge
            rets.append(self.plot([x1,x1],
                                  [ybottom,ytop],
                                  transform=transform,
                                  color=linecolor,
                                  linewidth=linewidth)[0])
            
            #plot right edge
            rets.append(self.plot([x4,x4],
                                  [ybottom,ytop],
                                  transform=transform,
                                  color=linecolor, 
                                  linewidth=linewidth)[0])
            
            #make a filled black box from left edge to 1/4 way across
            rets.append(ax.fill([x1,x2,x2,x1,x1],
                                [ytop,ytop,ybottom,ybottom,ytop],
                                transform=transform,
                                ec=fontcolor,
                                fc=fillcolor1)[0])
            
                #make a filled white box from 1/4 way across to 1/2 way across
            rets.append(ax.fill([x2,x0,x0,x2,x2],
                                [ytop,ytop,ybottom,ybottom,ytop],
                                transform=transform,
                                ec=fontcolor,
                                fc=fillcolor2)[0])
                
            #make a filled white box from 1/2 way across to 3/4 way across
            rets.append(ax.fill([x0,x3,x3,x0,x0],
                                [ytop,ytop,ybottom,ybottom,ytop],
                                transform=transform,
                                ec=fontcolor,
                                fc=fillcolor1)[0])
                
            #make a filled white box from 3/4 way across to end
            rets.append(ax.fill([x3,x4,x4,x3,x3],
                                [ytop,ytop,ybottom,ybottom,ytop],
                                transform=transform,
                                ec=fontcolor,
                                fc=fillcolor2)[0])
                
            #plot 3 tick marks at left edge, center, and right edge
            rets.append(self.plot([x1,x1],
                                  [ytick,ybottom],
                                  color=linecolor, 
                                  transform=transform,
                                  linewidth=linewidth)[0])
            
            rets.append(self.plot([x0,x0],
                                  [ytick,ybottom],
                                  transform=transform,
                                  color=linecolor, 
                                  linewidth=linewidth)[0])
            
            rets.append(self.plot([x4,x4],
                                  [ytick,ybottom],
                                  transform=transform,
                                  color=linecolor, 
                                  linewidth=linewidth)[0])
            
            #label 3 tick marks
            rets.append(ax.text(x1,lattext,format % (0),\
            horizontalalignment='center',\
            verticalalignment='top',\
            fontsize=fontsize,color=fontcolor))
            rets.append(ax.text(x0,lattext,format % (0.5*lenlab),\
            horizontalalignment='center',\
            verticalalignment='top',\
            fontsize=fontsize,color=fontcolor))
            rets.append(ax.text(x4,lattext,format % (lenlab),\
            horizontalalignment='center',\
            verticalalignment='top',\
            fontsize=fontsize,color=fontcolor))
            #put units, scale factor on top
            rets.append(ax.text(x0,ytop+yoffset/2,labelstr,\
            horizontalalignment='center',\
            verticalalignment='bottom',\
            fontsize=fontsize,color=fontcolor))
        # 'simple' style
        elif barstyle == 'simple':
            rets.append(self.plot([x1,x4],[yc,yc],color=linecolor, linewidth=linewidth)[0])
            rets.append(self.plot([x1,x1],[yc-yoffset,yc+yoffset],color=linecolor, linewidth=linewidth)[0])
            rets.append(self.plot([x4,x4],[yc-yoffset,yc+yoffset],color=linecolor, linewidth=linewidth)[0])
            rets.append(ax.text(xc,yc-yoffset,format % lenlab,\
            verticalalignment='top',horizontalalignment='center',\
            fontsize=fontsize,color=fontcolor))
            #put units, scale factor on top
            rets.append(ax.text(xc,yc+yoffset,labelstr,\
            horizontalalignment='center',\
            verticalalignment='bottom',\
            fontsize=fontsize,color=fontcolor))
        else:
            raise KeyError("barstyle must be 'simple' or 'fancy'")
        if zorder is not None:
            for ret in rets:
                try:
                    ret.set_zorder(zorder)
                except:
                    pass
        return rets

    def plot(self, *args, **kwargs):
        """
        Draw lines and/or markers on the map
        (see matplotlib.pyplot.plot documentation).
        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.
        Extra keyword ``ax`` can be used to override the default axis instance.
        Other \**kwargs passed on to matplotlib.pyplot.plot.
        """
        ax = self.ax
        self._save_use_hold(ax, kwargs)
        try:
            ret =  ax.plot(*args, 
                           **kwargs)
        finally:
            self._restore_hold(ax)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        return ret

    def _save_use_hold(self, ax, kwargs):
        h = kwargs.pop('hold', None)
        if hasattr(ax, '_hold'):
            self._tmp_hold = ax._hold
            if h is not None:
                ax._hold = h

    def _restore_hold(self, ax):
        if hasattr(ax, '_hold'):
            ax._hold = self._tmp_hold
    
    
    
    def set_axes_limits(self,ax=None):
        """
        Final step in Basemap method wrappers of Axes plotting methods:
        Set axis limits, fix aspect ratio for map domain using current
        or specified axes instance.  This is done only once per axes
        instance.
        In interactive mode, this method always calls draw_if_interactive
        before returning.
        """
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()

        # If we have already set the axes limits, and if the user
        # has not defeated this by turning autoscaling back on,
        # then all we need to do is plot if interactive.
        if (hash(ax) in self._initialized_axes
                                 and not ax.get_autoscalex_on()
                                 and not ax.get_autoscaley_on()):
            if is_interactive():
                import matplotlib.pyplot as plt
                plt.draw_if_interactive()
            return

        self._initialized_axes.add(hash(ax))
        # Take control of axis scaling:
        ax.set_autoscale_on(False)
        # update data limits for map domain.
        corners = ((self.xmin, self.ymin), (self.xmax, self.ymax))
        ax.update_datalim(corners)
        ax.set_xlim((self.xmin, self.xmax))
        ax.set_ylim((self.ymin, self.ymax))
        # if map boundary not yet drawn for elliptical maps, draw it with default values.

        # make sure aspect ratio of map preserved.
        # plot is re-centered in bounding rectangle.
        # (anchor instance var determines where plot is placed)
        if self.fix_aspect:
            ax.set_aspect('equal',anchor=self.anchor)
        else:
            ax.set_aspect('auto',anchor=self.anchor)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        # force draw if in interactive mode.
        if is_interactive():
            import matplotlib.pyplot as plt
            plt.draw_if_interactive()
    
    
    def _cliplimb(self,ax,coll):
        if not self._mapboundarydrawn:
            return coll, None
        c = self._mapboundarydrawn
        if c not in ax.patches:
            p = ax.add_patch(c)
            #p.set_clip_on(False)
        try:
            coll.set_clip_path(c)
        except:
            for item in coll:
                item.set_clip_path(c)
        return coll,c
    

# now the test


from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

import geopandas as gpd

def get_standard_gdf():
    """ basic function for getting some geographical data in geopandas GeoDataFrame python's instance:
        An example data can be downloaded from Brazilian IBGE:
        ref: ftp://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2017/Brasil/BR/br_municipios.zip    
    """
    gdf_path = r'C:\my_file_path\Shapefile.shp'

    return gpd.read_file(gdf_path)

def format_ax(ax, projection, xlim, ylim):

    ax.set_xlim(xlim)
    ax.set_ylim(ylim)

    ax.set_global()
    ax.coastlines()

def main():
    
    
    fig = plt.figure(figsize=(8, 10))

    # Label axes of a Plate Carree projection with a central longitude of 180:
    
    #for enum, proj in enumerate(['Mercator, PlateCarree']):
    
    gdf = get_standard_gdf()
    
    xmin, ymin, xmax, ymax = gdf.total_bounds
    xlim = [xmin,  xmax]
    ylim = [ymin,  ymax]
    
    
    lon_c = np.mean(xlim)
    lat_c = np.mean(ylim)
    
    projection = ccrs.PlateCarree(central_longitude=0)
    
    ax1 = fig.add_subplot(3, 1, 1,
                          projection=projection,
                          xlim=[xmin, xmax],
                          ylim=[ymin, ymax])
    
    
    gdf.plot(ax=ax1, transform=projection)
    
    format_ax(ax1, projection, xlim, ylim)
    
    
    Grider = ax1.gridlines(draw_labels=True)
    Grider.xformatter = LONGITUDE_FORMATTER
    Grider.yformatter = LATITUDE_FORMATTER
    Grider.xlabels_top  = False
    Grider.ylabels_right  = False
    
    
    # Label axes of a Mercator projection without degree symbols in the labels
    # and formatting labels to include 1 decimal place:
    ax2 = fig.add_subplot(3, 1, 2, 
                          projection=ccrs.Mercator(),
                          xlim=[xmin, xmax],
                          ylim=[ymin, ymax])
    gdf.plot(ax=ax2, transform=projection)
    
    
    format_ax(ax2, projection, xlim, ylim)
    
    
    Grider = ax2.gridlines(draw_labels=True)
    Grider.xformatter = LONGITUDE_FORMATTER
    Grider.yformatter = LATITUDE_FORMATTER
    Grider.xlabels_top  = False
    Grider.ylabels_right  = False
    
    
    ax3 = fig.add_subplot(3, 1, 3, 
                          projection=ccrs.Robinson(central_longitude=lon_c,
                                                   #central_latitude=lat_c
                                                   ),
                          xlim=[xmin, xmax],
                          ylim=[ymin, ymax])
    
    gdf.plot(ax=ax3, transform=projection)
    
    format_ax(ax3, projection, xlim, ylim)
    
    
    ax3.set_xticks([-180, -120, -60, 0, 60, 120, 180])
    ax3.set_yticks([-78.5, -60, -25.5, 25.5, 60, 80])

    ax3.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
    ax3.yaxis.set_major_formatter(LATITUDE_FORMATTER)
    plt.draw()
    return fig, fig.get_axes()

if __name__ == '__main__':
    
    length = 1000
    
    fig, axes = main()
    
    gdf = get_standard_gdf()
    
    xmin, ymin, xmax, ymax = gdf.total_bounds
    xoff = 0.3 * (xmax - xmin)
    yoff = 0.2 * (ymax - ymin)
    
    for ax in axes:
        if hasattr(ax, 'projection'):
            x0, x1, y0, y1 = np.ravel(ax.get_extent())
            
            Scaler = Scalebar(ax=ax,
                              metric_ccrs=ccrs.Geodetic())
            Scaler.drawmapscale(lon = xmin+xoff, 
                                lat = ymin + yoff,
                                
                                length=length,
                                units = 'km',
                                barstyle='fancy',
                                yoffset=0.2 * (ymax - ymin)
                                )
            
    fig.suptitle('Using Cartopy')
    fig.show()
    
    

运行上述代码时,比例尺在地理轴中错位。比例尺 xticks 错位,y 轴高度比例也不正确。

这是一个示例:geopandas 绘制为蓝色。请注意,比例尺仅在第二个和第三个地理轴中可见。

Global view with three different projections

matplotlib matplotlib-basemap geopandas cartopy pyproj
2个回答
0
投票

我找到了当前问题的解决方案。

为了简洁起见,代码显示在here

请随意检查一下。该算法仍然需要一些调整才能支持其他 cartopy 投影。

同时可以应用于PlateCarree投影。


0
投票

对于那些通过某种搜索引擎遇到这个问题的人:我最近开发了一个综合包来帮助解决这个问题:matplotlib-map-utils

它具有专用函数 (

scale_bar()
) 和类 (
ScaleBar
),用于创建具有大量自定义功能的对象。

# Importing
from matplotlib_map_utils.core.scale_bar import ScaleBar, scale_bar
# Setting up a plot
fig, ax = matplotlib.pyplot.subplots(1,1, figsize=(5,5), dpi=150)
# Adding a scale bar to the upper-right corner of the axis, 
# in the same projection as whatever geodata you plotted
# Here, this scale bar will have the "boxes" style
scale_bar(ax=ax, location="upper right", style="boxes", bar={"projection":3857})

rendered scale bar on a plot

这并没有真正触及多少东西可以自定义的表面:查看

docs\howto_scale_bar
中的 GitHub 存储库 以了解更多信息。

© www.soinside.com 2019 - 2024. All rights reserved.