Re: [Matplotlib-users] shadedrelief basemap question

[Jeff W. <jef...@no...>]

On 8/1/12 4:09 PM, Scott Henderson wrote:
> Hi Jeff,
>
> I'm continuing to enjoy using basemap, but have a question about the shaded relief background. I frequently use greyscale shaded relief on the continents, but blue or white for the oceans. the shadedrelief() function is really convenient, but it includes shading for the oceans. Is there a way to afterwards just shade oceans? I've included my modified shaded_relief function in this email that uses some of the other natural earth products in case you need to see it.
>
> Example:
> bmap.drawmapboundary(fill_color='aqua')
> bmap.shadedrelief()
>
> Or something like:
> bmap.shadedrelief()
> **bmap.filloceans('aqua')
>
>
> Thanks,
> Scott

Scott:  You could overlay the shaded relief image on a land-sea mask, 
where the land part of the mask is transparent.  Like this:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
bmap = Basemap(projection='ortho',lat_0=45,lon_0=-100,resolution='l')
bmap.drawmapboundary()
bmap.drawmeridians(np.arange(0,360,30))
bmap.drawparallels(np.arange(-90,90,30))
bmap.shadedrelief()
bmap.drawlsmask(ocean_color='aqua',land_color=(255,255,255,1))
plt.show()

>
> PS. is there a basemap mailing list I should send questions like this to?

No - just send questions to matplotlib-users (I've cc'ed the list)

-Jeff
>
>
> from PIL import Image
> from matplotlib.image import pil_to_array
>
>      def shadedrelief_highres(self,style='hypsometric', region=(-75.0, -35.0, -65.0, -15.0),ax=None,scale=None,**kwargs):
>          ''' can avoid Memory Error with warpimage() by pre-cutting high-res arrays with GDAL?
>              region = (minLon, minLat, maxLon, maxLat)'''
>          print 'NOTE: only set up for cyl coordinates with extent defined by region'
>          basedir = '/Users/scotthenderson/data/natural_earth/'
>          if style == 'hypsometric':
>              path = os.path.join(basedir,'HYP_HR_SR_OB_DR.tif')
>          elif style == 'natural':
>              path = os.path.join(basedir,'NE1_HR_LC_SR_W_DR.tif')
>          elif style == 'hillshade':
>              path = os.path.join(basedir,'SR_HR.tif')
>
>          outfile = os.path.join(basedir,'tmp.tif')
>          if os.path.isfile(outfile): os.remove(outfile)
>          os.system('gdalwarp -te {0} {1} {2} {3} {infile} {out}'.format(*region, infile=path, out=outfile))
>
>          #copy only relevant commands from warpimage imshow() overrides colorbar for rgb or PIL arrays
>          pilImage = Image.open(outfile)
>          if scale is not None:
>              w, h = pilImage.size
>              width = int(np.round(w*scale))
>              height = int(np.round(h*scale))
>              pilImage = pilImage.resize((width,height),Image.ANTIALIAS)
>          self._bm_rgba = pil_to_array(pilImage)
>          # if pil_to_array returns a 2D array, it's a grayscale image.
>          # create an RGB image, with R==G==B.
>          if self._bm_rgba.ndim == 2:
>              tmp = np.empty(self._bm_rgba.shape+(3,),np.uint8)
>              for k in range(3):
>                  tmp[:,:,k] = self._bm_rgba
>              self._bm_rgba = tmp
>
>          im = self.imshow(self._bm_rgba, ax=ax, **kwargs)
>          return im


-- 
Jeffrey S. Whitaker         Phone  : (303)497-6313
Meteorologist               FAX    : (303)497-6449
NOAA/OAR/PSD  R/PSD1        Email  : Jef...@no...
325 Broadway                Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web    : http://tinyurl.com/5telg