Python colors.LightSource类代码示例

def test_light_source_shading_color_range():
    # see also
    #http://matplotlib.org/examples/pylab_examples/shading_example.html

    from matplotlib.colors import LightSource
    from matplotlib.colors import Normalize

    refinput = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
    norm = Normalize(vmin=0, vmax=50)
    ls = LightSource(azdeg=0, altdeg=65)
    testoutput = ls.shade(refinput, plt.cm.jet, norm=norm)
    refoutput = np.array([
        [[0., 0., 0.58912656, 1.],
        [0., 0., 0.67825312, 1.],
        [0., 0., 0.76737968, 1.],
        [0., 0., 0.85650624, 1.]],
        [[0., 0., 0.9456328, 1.],
        [0., 0., 1., 1.],
        [0., 0.04901961, 1., 1.],
        [0., 0.12745098, 1., 1.]],
        [[0., 0.22156863, 1., 1.],
        [0., 0.3, 1., 1.],
        [0., 0.37843137, 1., 1.],
        [0., 0.45686275, 1., 1.]]
        ])
    assert_array_almost_equal(refoutput, testoutput)

def main():

    fig, ax = plt.subplots(3, 2, figsize=(7, 10))
    fig.tight_layout()    

    if 1:
        # Same terrain and data.
        data = make_test_data('circles', noise_factor=0.05) * 2 - 7
        terrain = data
    else:
        # Different terrain and data. Matplotlib hill shading can only show the data.
        data = make_test_data('hills') * -1
        terrain = make_test_data('circles', noise_factor=0.05) * 2
        
    assert terrain.shape == data.shape, "{} != {}".format(terrain.shape, data.shape)
    print("data range: {} {}".format(np.min(data), np.max(data)))
    
    if len(sys.argv) > 1:
        cmap_name = sys.argv[1]
    else:
        #cmap_name = 'copper' # from http://matplotlib.org/examples/pylab_examples/shading_example.html?highlight=codex%20shade
        #cmap_name = 'gist_earth' # works reasonably fine with all of them
        cmap_name = 'bwr'        # shows that mpl & pegtop don't work when there is no increasing intensity
        #cmap_name = 'cubehelix'  # shows that HSV blending does not work were color is black
        #cmap_name = 'rainbow'    # shows that mpl & pegtop don't work when there is no increasing intensity
        #cmap_name = 'Paired_r'   # is nice to inspect when the data is different from the terrain
        
    print ("Using colormap: {!r}".format(cmap_name))
    cmap = plt.cm.get_cmap(cmap_name)
    cmap.set_bad('cyan')
    cmap.set_over('cyan')
    cmap.set_under('cyan')
    
    abs_max = np.max(np.abs(data)) # force color bar to be symmetrical. 
    norm = mpl.colors.Normalize(vmin=-abs_max, vmax=abs_max)
    #norm = mpl.colors.Normalize(vmin=-2, vmax=3)
    
    draw(ax[0, 0], cmap=cmap, norm=norm, title='No shading', 
         image_data = data)

    draw(ax[0, 1], cmap=plt.cm.gray, title='Matplotlib intensity', 
         image_data = mpl_surface_intensity(terrain))

    ls = LightSource(azdeg=DEF_AZIMUTH, altdeg=DEF_ELEVATION)
    draw(ax[1, 0], cmap=cmap, norm=norm, title='Matplotlib hill shading', 
         image_data = ls.shade(data, cmap=cmap, norm=norm))
    
    draw(ax[1, 1], cmap=cmap, norm=norm, title='Pegtop blending', 
         image_data = mpl_hill_shade(data, terrain=terrain,  
                                     cmap=cmap, norm=norm, blend_function=pegtop_blending))
    
    draw(ax[2, 0], cmap=cmap, norm=norm, title='RGB blending', 
         image_data = mpl_hill_shade(data, terrain=terrain, 
                                     cmap=cmap, norm=norm, blend_function=rgb_blending))    
    
    draw(ax[2, 1], cmap=cmap, norm=norm, title='HSV blending', 
         image_data = mpl_hill_shade(data, terrain=terrain, 
                                     cmap=cmap, norm=norm, blend_function=hsv_blending))    

    plt.show()

    def plot_shaded(self, cmap=plt.cm.terrain, lightsource_kwargs=None, *args, **kwargs):

        if lightsource_kwargs is None:
            lightsource_kwargs = {'azdeg':225, 'altdeg':5}

        extent = [self.x.min(), self.x.max(), self.y.min(), self.y.max()]

        arr = self.z.copy()
        nan_mask = np.isnan(arr)
        arr_min = arr[~nan_mask].min()
        if nan_mask.any():
            arr[nan_mask] = max(arr_min-10, 0)

        ls = LightSource(**lightsource_kwargs)
        shaded = ls.shade(arr, cmap=cmap)

        fig = plt.figure()
        ax = fig.add_subplot(111)

        im = ax.imshow(shaded, cmap=cmap, extent=extent, *args, **kwargs)
        plt.colorbar(im)

        ax.get_xaxis().get_major_formatter().set_useOffset(False)
        ax.get_yaxis().get_major_formatter().set_useOffset(False)

        plt.show()

 def background_map(self, ax):
     llcrnrlat, urcrnrlat, llcrnrlon, urcrnrlon, lat_ts = (31, 44, -126, -113, 37.5)
     m = Basemap(projection='merc', llcrnrlat=llcrnrlat,
                 urcrnrlat=urcrnrlat, llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon,
                 lat_ts=lat_ts, resolution='i', ax=ax)
     m.drawmapboundary(fill_color='lightblue', zorder=0)
     m.fillcontinents(zorder=0)
     etopofn = '/home/behry/uni/data/etopo1_central_europe_gmt.grd'
     etopodata = Dataset(etopofn, 'r')
     z = etopodata.variables['z'][:]
     x_range = etopodata.variables['x_range'][:]
     y_range = etopodata.variables['y_range'][:]
     spc = etopodata.variables['spacing'][:]
     lats = np.arange(y_range[0], y_range[1], spc[1])
     lons = np.arange(x_range[0], x_range[1], spc[0])
     topoin = z.reshape(lats.size, lons.size, order='C')
     # transform to nx x ny regularly spaced 5km native projection grid
     nx = int((m.xmax - m.xmin) / 5000.) + 1; ny = int((m.ymax - m.ymin) / 5000.) + 1
     topodat, x, y = m.transform_scalar(np.flipud(topoin), lons, lats, nx, ny, returnxy=True)
     ls = LightSource(azdeg=300, altdeg=15, hsv_min_sat=0.2, hsv_max_sat=0.3,
                      hsv_min_val=0.2, hsv_max_val=0.3)
     # shade data, creating an rgb array.
     rgb = ls.shade(np.ma.masked_less(topodat / 1000.0, 0.0), cm.gist_gray_r)
     m.imshow(rgb)
     m.drawmeridians(np.arange(6, 12, 2), labels=[0, 0, 0, 1], color='white',
                     linewidth=0.5, zorder=0)
     m.drawparallels(np.arange(44, 50, 2), labels=[1, 0, 0, 0], color='white',
                     linewidth=0.5, zorder=0)
     m.drawcoastlines(zorder=1)
     m.drawcountries(linewidth=1.5, zorder=1)
     m.drawstates()
     m.drawrivers(color='lightblue', zorder=1)
     return m

def SurfPlot(data):
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib import cbook
    from matplotlib import cm
    from matplotlib.colors import LightSource
    import matplotlib.pyplot as plt
    import numpy as np

    ncols, nrows = data.shape
    if ncols * nrows > 20000: print 'Warning, surface plotting things of this size is slow...'
    z = data
    x = np.linspace(0, ncols, ncols)
    y = np.linspace(0, nrows, nrows)
    x, y = np.meshgrid(x, y)
    
    region = np.s_[0:min(ncols,nrows), 0:min(ncols,nrows)]
    x, y, z = x[region], y[region], z[region]

    fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(10,10))

    ls = LightSource(270, 45)
    # To use a custom hillshading mode, override the built-in shading and pass
    # in the rgb colors of the shaded surface calculated from "shade".
    rgb = ls.shade(z, cmap=cm.CMRmap, vert_exag=0.01, blend_mode='soft')
    surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=rgb, linewidth=0, antialiased=True, shade=False)
    plt.show()

    def apply(self):
        from mpl_toolkits.basemap import shiftgrid, cm
        common.ShowQuestion
        display = self.VpyartDisplay.value
        if (isinstance(display, pyart.graph.RadarMapDisplay) or
            isinstance(display, pyart.graph.GridMapDisplay)):
            pass
        elif (isinstance(display, pyart.graph.RadarDisplay) or
              isinstance(display, pyart.graph.AirborneRadarDisplay)):
            common.ShowWarning(
                "Topography require a MapDisplay, be sure to "
                "check the 'use MapDisplay' box in the 'Display Options' Menu")
            return
        else:
            common.ShowWarning(
                "Need a pyart display instance, be sure to "
                "link this components (%s), to a radar or grid display" %
                self.name)
            return

        filename = str(self.lineEdit.text())
        if filename != self.current_open:
            if filename.startswith("http"):
                resp = common.ShowQuestion(
                    "Loading a file from the internet may take long." +
                    " Are you sure you want to continue?")
                if resp != QtWidgets.QMessageBox.Ok:
                    return
            self.etopodata = Dataset(filename)
            self.current_open = filename

        topoin = np.maximum(0, self.etopodata.variables['ROSE'][:])
        lons = self.etopodata.variables['ETOPO05_X'][:]
        lats = self.etopodata.variables['ETOPO05_Y'][:]
        # shift data so lons go from -180 to 180 instead of 20 to 380.
        topoin, lons = shiftgrid(180., topoin, lons, start=False)

        # plot topography/bathymetry as an image.

        # create the figure and axes instances.
        # setup of basemap ('lcc' = lambert conformal conic).
        # use major and minor sphere radii from WGS84 ellipsoid.
        m = self.VpyartDisplay.value.basemap
        # transform to nx x ny regularly spaced 5km native projection grid
        nx = int((m.xmax - m.xmin)/500.) + 1
        ny = int((m.ymax - m.ymin)/500.) + 1
        topodat = m.transform_scalar(topoin, lons, lats, nx, ny)
        # plot image over map with imshow.

        # draw coastlines and political boundaries.

        ls = LightSource(azdeg=90, altdeg=20)
        # convert data to rgb array including shading from light source.
        # (must specify color map)
        rgb = ls.shade(topodat, cm.GMT_relief)
        im = m.imshow(rgb)

        self.VpyartDisplay.update(strong=False)

    def draw_figure(self):
        ls = LightSource(azdeg=315, altdeg=45)
        self.surf = self.data

        # try to plot the whole array
        try:
            blended_surface = ls.shade(self.surf, cmap=self.colormap, vert_exag=5, blend_mode=b"overlay",
                                       vmin=np.nanmin(self.surf), vmax=np.nanmax(self.surf))
        # too big, two attempts for sub-sampling
        except MemoryError:
            rows = self.data.shape[0]
            cols = self.data.shape[1]

            # 1st attempt: <= 1024x1024
            try:
                max_elements = 1024
                row_stride = rows // max_elements + 1
                col_stride = cols // max_elements + 1
                self.surf = self.data[::row_stride, ::col_stride]
                blended_surface = ls.shade(self.surf, cmap=self.colormap, vert_exag=5, blend_mode=b"overlay",
                                           vmin=np.nanmin(self.surf), vmax=np.nanmax(self.surf))

            except MemoryError:
                # 2st attempt: <= 512x512
                max_elements = 512
                row_stride = rows // max_elements + 1
                col_stride = cols // max_elements + 1
                self.surf = self.data[::row_stride, ::col_stride]
                blended_surface = ls.shade(self.surf, cmap=self.colormap, vert_exag=5, blend_mode=b"overlay",
                                           vmin=np.nanmin(self.surf), vmax=np.nanmax(self.surf))

            log.debug("too big: %s x %s > subsampled to %s x %s"
                      % (self.data.shape[0], self.data.shape[1], self.surf.shape[0], self.surf.shape[1]))

        self.axes.coastlines(resolution='50m', color='gray', linewidth=1)
        img = self.axes.imshow(blended_surface, origin='lower', cmap=self.colormap,
                               extent=self.geo_extent, transform=ccrs.PlateCarree())
        img.set_clim(vmin=np.nanmin(self.surf), vmax=np.nanmax(self.surf))
        # add gridlines with labels only on the left and on the bottom
        grl = self.axes.gridlines(crs=ccrs.PlateCarree(), color='gray', draw_labels=True)
        grl.xformatter = LONGITUDE_FORMATTER
        grl.yformatter = LATITUDE_FORMATTER
        grl.xlabel_style = {'size': 8}
        grl.ylabel_style = {'size': 8}
        grl.ylabels_right = False
        grl.xlabels_top = False

        if self.cb:
            self.cb.on_mappable_changed(img)
        else:
            self.cb = plt.colorbar(img, ax=self.axes)
        self.cb.ax.tick_params(labelsize=8)

 def plotRainbow(self):
     
     self.setUp()
     
     ls = LightSource(270, 45)
     # To use a custom hillshading mode, override the built-in shading and pass
     # in the rgb colors of the shaded surface calculated from "shade".
     rgb = ls.shade(self.Z, cmap=cm.gist_earth, vert_exag=0.1, blend_mode='soft')
     
     self.ax.plot_surface(self.X, self.Y, self.Z, facecolors=rgb)
     
     self.setLabels()
     plt.show()

def getTopoRGB(topogrid):
    topotmp = topogrid.getData().copy()
    # make a masked array
    topotmp = np.ma.array(topotmp)
    topodat = np.ma.masked_where(np.isnan(topotmp), topotmp)
    cy = topogrid.geodict["ymin"] + (topogrid.geodict["ymax"] - topogrid.geodict["ymin"]) / 2.0
    # flag the regions where topography is less than 0 (we'll color this ocean later)
    i = np.where(topodat == SEA_LEVEL)

    # do shaded relief stuff
    # keys are latitude values
    # values are multiplication factor
    zdict = {
        0: 0.00000898,
        10: 0.00000912,
        20: 0.00000956,
        30: 0.00001036,
        40: 0.00001171,
        50: 0.00001395,
        60: 0.00001792,
        70: 0.00002619,
        80: 0.00005156,
    }
    # find the mean latitude of the map we're making, and use that to come up with a zfactor
    mlat = abs(int(round(cy / 10) * 10))
    zfactor = zdict[mlat]
    ls = LightSource(azdeg=AZDEFAULT, altdeg=ALTDEFAULT)

    # draw the ocean in light blue
    water_color = [0.47, 0.60, 0.81]
    palette1 = copy.deepcopy(cm.binary)

    # draw the light shaded-topography
    rgbdata = topodat * zfactor  # apply the latitude specific zfactor correction
    rgb = ls.shade(rgbdata, cmap=palette1)  # apply the light shading to our corrected topography

    # this is an rgb data set now, so masking the pixels won't work, but explicitly setting all of the
    # "bad" pixels to our water color will
    red = rgb[:, :, 0]
    green = rgb[:, :, 1]
    blue = rgb[:, :, 2]
    red[i] = water_color[0]
    green[i] = water_color[1]
    blue[i] = water_color[2]
    rgb[:, :, 0] = red
    rgb[:, :, 1] = green
    rgb[:, :, 2] = blue

    rgb = np.flipud(rgb)

    return (rgb, palette1)

def shade_other_data():
    """Demonstrates displaying different variables through shade and color."""
    y, x = np.mgrid[-4:2:200j, -4:2:200j]
    z1 = np.sin(x**2)  # Data to hillshade
    z2 = np.cos(x**2 + y**2)  # Data to color

    norm = Normalize(z2.min(), z2.max())
    cmap = plt.cm.RdBu

    ls = LightSource(315, 45)
    rgb = ls.shade_rgb(cmap(norm(z2)), z1)

    fig, ax = plt.subplots()
    ax.imshow(rgb)
    ax.set_title('Shade by one variable, color by another', size='x-large')

    def _getShaded(self,ptopo):
        maxvalue = self.contour_colormap.vmax
        ls1 = LightSource(azdeg = 120, altdeg = 45)
        ls2 = LightSource(azdeg = 225, altdeg = 45)
        intensity1 = ls1.hillshade(ptopo, fraction = 0.25, vert_exag = VERT_EXAG)
        intensity2 = ls2.hillshade(ptopo, fraction = 0.25, vert_exag = VERT_EXAG)
        intensity = intensity1*0.5 + intensity2*0.5

        ptoposc = ptopo/maxvalue
        rgba = self.contour_colormap.cmap(ptoposc)
        rgb = np.squeeze(rgba)

        draped_hsv = ls1.blend_hsv(rgb,np.expand_dims(intensity,2))
        
        return draped_hsv

class LightFilter(BaseFilter):
    "simple gauss filter"

    def __init__(self, sigma, fraction=0.5):
        self.gauss_filter = GaussianFilter(sigma, alpha=1)
        self.light_source = LightSource()
        self.fraction = fraction
        #hsv_min_val=0.5,hsv_max_val=0.9,
        #                                hsv_min_sat=0.1,hsv_max_sat=0.1)

    def get_pad(self, dpi):
        return self.gauss_filter.get_pad(dpi)

    def process_image(self, padded_src, dpi):
        t1 = self.gauss_filter.process_image(padded_src, dpi)
        elevation = t1[:,:,3]
        rgb = padded_src[:,:,:3]

        rgb2 = self.light_source.shade_rgb(rgb, elevation,
                                           fraction=self.fraction)

        tgt = np.empty_like(padded_src)
        tgt[:,:,:3] = rgb2
        tgt[:,:,3] = padded_src[:,:,3]

        return tgt

def display_depth_matplotlib(z):
    """
    Same as above but using matplotlib instead.
    """
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib.colors import LightSource
    
    m, n = z.shape
    x, y = np.mgrid[0:m, 0:n]
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ls = LightSource(azdeg=0, altdeg=65)
    greyvals = ls.shade(z, plt.cm.Greys)
    ax.plot_surface(x, y, z, rstride=1, cstride=1, linewidth=0, antialiased=False, facecolors=greyvals)
    plt.axis('off')
    plt.axis('equal')
    plt.show()

def display_colorbar():
    """Display a correct numeric colorbar for a shaded plot."""
    y, x = np.mgrid[-4:2:200j, -4:2:200j]
    z = 10 * np.cos(x**2 + y**2)

    cmap = plt.cm.copper
    ls = LightSource(315, 45)
    rgb = ls.shade(z, cmap)

    fig, ax = plt.subplots()
    ax.imshow(rgb)

    # Use a proxy artist for the colorbar...
    im = ax.imshow(z, cmap=cmap)
    im.remove()
    fig.colorbar(im)

    ax.set_title('Using a colorbar with a shaded plot', size='x-large')

def main3():
    ax = plt.axes(projection=ccrs.PlateCarree())

    elev, crs, extent = io_srtm.srtm_composite(-5, 52, 2, 2)

    elev = np.ma.masked_less_equal(elev, 0, copy=False)

    use_mpl_light_source = False
    
    if use_mpl_light_source:
        from matplotlib.colors import LightSource
    
        ls = LightSource(azdeg=90, altdeg=80,
                         hsv_min_val=0.6, hsv_min_sat=0.9,
                         hsv_max_val=0.8, hsv_max_sat=1,
                         )
    
        rgb = ls.shade(elev, plt.get_cmap('Greens', 3))
    else:
        import matplotlib.colors as mcolors
        rgb = set_shade(elev,
                        cmap=mcolors.ListedColormap([plt.get_cmap('Greens', 3)(0.5)])
                    )


    ax.imshow(rgb,
                extent=extent,
                transform=crs
                )

    x = np.linspace(extent[0], extent[1], elev.shape[0])
    y = np.linspace(extent[2], extent[3], elev.shape[1])
#
#    ax.contour(x, y, elev, 100,
#               linestyles='-',
#               colors='blue',
#               linewidths=0.3,
#               alpha=0.4,
#               transform=crs,
#               )

    plt.show()