Changeset 2742 for TI05-delivery


Ignore:
Timestamp:
25/07/07 12:42:08 (12 years ago)
Author:
spascoe
Message:

Using Basemap to transform the grid into Mercator projection.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • TI05-delivery/ows_framework/branches/ows_framework-ddp/ows_server/ows_server/lib/ddp_render.py

    r2710 r2742  
    55 
    66import pylab, cdms 
     7from matplotlib.toolkits.basemap import Basemap 
    78from cdms import selectors 
    89import Numeric as N 
     
    1920 
    2021def grid_corners_to_mesh(grid_corner_lat, grid_corner_lon, 
    21                          bl_index=0): 
     22                         bl_index=0, basemap=None): 
    2223    """ 
    2324    Converts a scrip-style grid description into a matplotlib 
     
    3536    @param bl_index: The index in grid_corner_* representing the 
    3637        minimum lat/lon corner of the rotated grid. 
     38    @param basemap: If not None this is a Basemap instance used to 
     39        transform the mesh to projection coordinates. 
    3740 
    3841    @param return: (mesh_x, mesh_y) 
     
    6972    mesh_y[-1,-1] = grid_corner_lat[-1,-1, tr_index] 
    7073 
    71     return (mesh_x, mesh_y) 
     74    if basemap: 
     75        return basemap(mesh_x, mesh_y) 
     76    else: 
     77        return (mesh_x, mesh_y) 
    7278 
    7379class RotatedGridRenderer(object): 
     
    8288 
    8389    """ 
    84     def __init__(self, shape, grid_corner_lat, grid_corner_lon): 
     90    def __init__(self, shape, grid_corner_lat, grid_corner_lon, basemap=None): 
    8591        """ 
    8692        @param shape: The shape of required grid 
     
    96102        g_lon = N.reshape(grid_corner_lon, g_shape) 
    97103 
    98         self.mesh_x, self.mesh_y = grid_corners_to_mesh(g_lat, g_lon, 3) 
     104        self.mesh_x, self.mesh_y = grid_corners_to_mesh(g_lat, g_lon, 3, 
     105                                                        basemap=basemap) 
    99106 
    100107    def plotSlice(self, var, minValue, maxValue): 
     
    141148 
    142149def main(filename, grid_filename, output_prefix): 
     150 
     151    bbox = (-16.5, 47.0, 4.5, 61.5) 
     152 
     153    # Create the map projection 
     154    bm = Basemap(projection='merc', 
     155                 resolution='l', 
     156                 llcrnrlon=bbox[0], 
     157                 llcrnrlat=bbox[1], 
     158                 urcrnrlon=bbox[2], 
     159                 urcrnrlat=bbox[3], 
     160                 lat_ts=0, 
     161                 ) 
     162 
    143163    # Some configuration options 
    144164    dpi = 100 
    145     xlim = [-16.5,4.5] 
    146     ylim = [47,61.5] 
     165    xlim, ylim = bm(bbox[::2], bbox[1::2]) 
    147166     
    148167    # Initialise logging 
     
    154173    r = None 
    155174    gridshape = None 
     175 
    156176 
    157177    for variable in f.listvariables(): 
     
    176196            logger.info('Initialising gridbox coordinates for grid %s' % (gridshape,)) 
    177197            r = RotatedGridRenderer(gridshape, 
    178                                     g['grid_corner_lat'], g['grid_corner_lon']) 
     198                                    g['grid_corner_lat'], g['grid_corner_lon'], 
     199                                    basemap=bm) 
    179200 
    180201 
     
    192213        # Loop over each time point 
    193214        for i in r.iterPlotTimes(var, minValue, maxValue): 
     215            bm.drawcoastlines() 
    194216            logger.info('Plotted index %d of %s' % (i, variable)) 
    195217            figname = '%s_%s_%s.png' % (output_prefix, variable, i) 
Note: See TracChangeset for help on using the changeset viewer.