Python Basemap.readshapefile方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def build_basemap():
    shapefile = DUBLIN_DED
    shp = fiona.open(shapefile+".shp")
    bds = shp.bounds
    shp.close()
    extra = 0.01
    ll = (bds[0], bds[1])
    ur = (bds[2], bds[3])
    coords = list(chain(ll, ur))
    w, h = coords[2] - coords[0], coords[3] - coords[1]

    m = Basemap(
    projection='tmerc',
    lon_0=-2.,
    lat_0=49.,
    ellps = 'WGS84',
    llcrnrlon=coords[0] - extra * w,
    llcrnrlat=coords[1] - extra + 0.01 * h,
    urcrnrlon=coords[2] + extra * w,
    urcrnrlat=coords[3] + extra + 0.01 * h,
    lat_ts=0,
    resolution='i',
    suppress_ticks=True)

    m.readshapefile(
        shapefile,
        'dublin',
        color='none',
        zorder=2)

    return m, coords

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def ResearchRegion_surface():
    """
    在地图上画出柱表面混合比图
    :return:
    """

    fig = plt.figure(figsize=(11, 8), facecolor="white")
    # data = np.loadtxt('seasonAvr_data/SurfaceMixingRatio/1_seasonAvr.txt')
    data = np.loadtxt("allYearAvr_data/SurfaceMixingRatio/allYearAvr.txt")
    arr = np.zeros((180, 360))
    for i in range(180):
        arr[i, :] = data[179 - i, :]

    longitude = np.loadtxt("lonlat_data/longitude.txt")
    latitude = np.loadtxt("lonlat_data/latitude.txt")

    m = Basemap(llcrnrlon=70, llcrnrlat=15, urcrnrlon=138, urcrnrlat=55, projection="mill", resolution="h")
    m.drawparallels(np.arange(5.5, 90.5, 1.0), color="w", linewidth=0.5, dashes=[1, 1], labels=[0, 0, 0, 0])
    m.drawmeridians(np.arange(60.5, 181.5, 1.0), color="w", linewidth=0.5, dashes=[1, 1], labels=[0, 0, 0, 0])
    m.drawmapboundary(fill_color="0.3")
    m.readshapefile("shp/CHINA", "CHINA", drawbounds=1, color="black")

    topo = maskoceans(longitude, latitude, arr)
    im = m.pcolormesh(longitude, latitude, topo, shading="flat", cmap=plt.cm.jet, latlon=True, vmin=0, vmax=500)
    m.drawlsmask(ocean_color="w", lsmask=0)

    cbar = m.colorbar()
    cbar.ax.set_ylabel("SurfaceMixingRatio", color="black", fontsize="14", rotation=90)
    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def graph(data_dict):
    # normalize data
    norm_data_dict = normalize_dict(data_dict)

    # Set up the basic map
    map = Basemap()
    map.drawmapboundary(fill_color=WATER_COLOR)
    # map.fillcontinents(color='#ddaa66',lake_color=WATER_COLOR)

    # Read in the counties
    map.readshapefile(SHAPEFILE, 'us_counties', drawbounds = False)

    for info, shape in zip(map.us_counties_info, map.us_counties):
        ax.add_patch(
            Polygon(
                np.array(shape),
                True,
                facecolor=get_color(
                    norm_data_dict.get(
                        (info['STATEFP'], info['COUNTYFP']),
                        -1
                    )
                ),
                edgecolor='k',
                linewidth=1.
            )
        )

    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
    def draw_samples(self, n_samples, burnin=100, thin=10):
        if not self.kde_model:
            self.fit_kde()

        # pick starting point
        start_x = np.array([[0., 31.]])

        # get region for restricting samples
        margin = 0.5
        m = Basemap(llcrnrlon=self.uganda_data.LONGITUDE.min() - margin,
            llcrnrlat=self.uganda_data.LATITUDE.min() - margin,
            urcrnrlon=self.uganda_data.LONGITUDE.max() + margin,
            urcrnrlat=self.uganda_data.LATITUDE.max() + margin,
            resolution='l',
            area_thresh=10000)
        m.readshapefile("data/regions/UGA_adm0", "regions", drawbounds=True)
        for xy, info in zip(m.regions, m.regions):
            p = path.Path(xy)

        slice_samples, _ = SliceSampler.SliceSampler.mvslice(self.kde_model.pdf, 
            start_x,
            sample_size=n_samples, 
            burnin=burnin, 
            thin=thin, 
            confirm_region=p)
        
        # we get lat, long; we want x, y
        return np.fliplr(slice_samples)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def main():
    m = Basemap(
        projection='merc',
        ellps='WGS84',
        llcrnrlon=lower_left[0],
        llcrnrlat=lower_left[1],
        urcrnrlon=upper_right[0],
        urcrnrlat=upper_right[1],
        lat_ts=0,
        resolution='i')
    m.readshapefile(
        'taxi_zone/taxi_zones_lat_lon',
        'nyc',
        color='none',
        zorder=2)

    predicted_values = read_file("./predict/ground_truth.txt")

    #pickups = pickle.load(open('data/30_unnormalized.p', 'rb'))

    #ground_truth = read_file('./predict/ground_truth.txt')


    gen_heat_map(m, predicted_values)
    
    print 'finished'

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def build_map_obj(shapefile):
    shp = fiona.open(shapefile + '.shp')
    bds = shp.bounds

    shp.close()
    extra = 0.01
    ll = (bds[0], bds[1])
    ur = (bds[2], bds[3])
    coords = list(chain(ll, ur))

    w, h = coords[2] - coords[0], coords[3] - coords[1]

    m = Basemap(
                projection='tmerc',
                lon_0=(coords[0] + coords[2])/2,
                lat_0=(coords[1] + coords[3])/2,
                ellps='helmert',
                llcrnrlon=coords[0] - w * 0.01,
                llcrnrlat=coords[1] - extra + 0.01 * h,
                urcrnrlon=coords[2] + w * 0.01,
                urcrnrlat=coords[3] + extra + 0.01 * h,
                resolution='i',
                suppress_ticks=True
                )

    m.readshapefile(shapefile,
                    'SF',
                    color='black',
                    zorder=2
                    )
    return m, coords

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
    def crop(self):
        cc = [49.4, 2.420480, 51.726419, 6.589545]  # belgium
        x0 = cc[1]
        y0 = cc[0]
        x1 = cc[3]
        y1 = cc[2]

        bmap = Basemap(resolution='i', llcrnrlat=y0, llcrnrlon=x0, urcrnrlat=y1, urcrnrlon=x1)
        bmap.land = None

        for path in self.shapefile_paths:
            bmap.readshapefile(path, 'land')
            self.polys += [shape2poly(shape) for shape in bmap.land]

        with open(self.in_file, 'r') as csvfile:
            csvreader = csv.reader(csvfile, delimiter=',')
            header = True
            for row in csvreader:
                if header:
                    self.header = ','.join(row)
                    header = False
                else:
                    p1 = Point(float(row[0]), float(row[1]))
                    p2 = Point(float(row[2]), float(row[3]))
                    self.segments.append((p1, p2, row))

        self.segments = [self.fix(segment) for segment in self.segments]
        # flatten
        self.segments = [item for sublist in self.segments for item in sublist]    

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def plot_ppi(sweep_dict, parm, **kwargs):
	radar_loc=kwargs.get('radar_loc', [-12.2492,  131.0444])
	fig_name=kwargs.get('fig_name', 'ppi_'+parm+'_.png')
	fig_path=kwargs.get('fig_path', getenv('HOME')+'/bom_mds/output/')
	f=figure()
	#gp_loc=[-12.2492,  131.0444]
	Re=6371.0*1000.0
	rad_at_radar=Re*sin(pi/2.0 -abs(radar_loc[0]*pi/180.0))#ax_radius(float(lat_cpol), units='degrees')
	lons=radar_loc[1]+360.0*sweep_dict['xar']/(rad_at_radar*2.0*pi)
	lats=radar_loc[0] + 360.0*sweep_dict['yar']/(Re*2.0*pi)
	def_loc_dict={'lat_0':lats.mean(), 'lon_0':lons.mean(),'llcrnrlat':lats.min(), 'llcrnrlon':lons.min(), 'urcrnrlat':lats.max() , 'urcrnrlon':lons.max(), 'lat_ts':lats.mean()}
	loc_dict=kwargs.get('loc_dict', def_loc_dict)
	map= Basemap(projection='merc', resolution='l',area_thresh=1., **loc_dict)
	xx, yy = map(lons, lats)
	#map.drawcoastlines()
	#map.drawcountries()
	map.drawmapboundary()
	map.readshapefile(getenv('HOME')+'/bom_mds/shapes/cstntcd_r','coast',drawbounds=True, linewidth=0.5,color='k',antialiased=1,ax=None)
	map.drawmeridians(array([129,130,131,132,133]), labels=[1,0,0,1])
	map.drawparallels(array([-14,-13,-12,-11,-10]), labels=[1,0,0,1])
	levs_dict={'VR':linspace(-15,15,31), 'CZ': linspace(-8,64,10), 'PH': linspace(0,185,255), "RH": linspace(0,1.5,16), "SW":linspace(0, 5, 11), "ZD":linspace(-10,10,21), 'VE':linspace(-30,30,31), 'TI':linspace(-30,30,31), 'KD':linspace(-1.0,6.0,30)}
	titles_dict={'VR': 'Velocity m/s', 'CZ':'Corrected Reflectivity dBz', 'PH': 'Differential Prop Phase (degrees)', 'RH':'Correlation Co-ef', 'SW':'Spectral Width (m/s)', 'ZD':'Differentail Reflectivity dBz', 'VE':'Edited Velocity (m/s)', 'TI':'Simualated winds,(m/s)', 'KD':'Specific differential Phase (Degrees/m)'}
	map.contourf(xx,yy,sweep_dict[parm], levels=levs_dict[parm])
	p=sweep_dict['date']
	dtstr='%(#1)02d-%(#2)02d-%(#3)04d %(#4)02d%(#5)02dZ ' %{"#1":p.day, "#2":p.month, "#3":p.year, "#4":p.hour, '#5':p.minute}
	title(sweep_dict['radar_name']+' '+dtstr+titles_dict[parm])
	colorbar()
	savefig(fig_path+fig_name)
	close(f)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
    def plotmap(self,fig,ax):
        """ This function will plot the map of Alaska. The data will be plotted
            over it and will use the basemap class to position everything.
            Input
                fig - The figure handle for the plots.
                ax - The axes handle that the map will be plotted over.
            Output
                m - This is the handle for the basemap object.
        """
        latlim2 = self.params['latbounds']
        lonlim2 = self.params['lonbounds']
        m = Basemap(projection='merc',lon_0=sp.mean(lonlim2),lat_0=sp.mean(latlim2),\
        lat_ts=sp.mean(latlim2),llcrnrlat=latlim2[0],urcrnrlat=latlim2[1],\
        llcrnrlon=lonlim2[0],urcrnrlon=lonlim2[1],\
        rsphere=6371200.,resolution='i',ax=ax)
        # draw coastlines, state and country boundaries, edge of map.
        #m.drawcoastlines()
    #    m.drawstates()
    #    m.drawcountries()
        m.readshapefile('st99_d00','states',drawbounds=True)

        merstep = sp.round_((lonlim2[1]-lonlim2[0])/5.)
        parstep = sp.round_((latlim2[1]-latlim2[0])/5.)
        meridians=sp.arange(lonlim2[0],lonlim2[1],merstep)
        parallels = sp.arange(latlim2[0],latlim2[1],parstep)
        m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
        m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
        plt.hold(True)
        return m

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def create_map():
    poverty_data, error = world_bank_api.most_recent_poverty_data()
    if error is not None:
        print(error)
        return

    poverty_data = utils.dictify_list_of_dicts(poverty_data, 'country_code')

    world_map = Basemap()
    world_map.readshapefile('borders', 'countries')
    country_to_color, country_names = decide_colors_and_countries(world_map, poverty_data)

    axes = plt.gca() # get current axes instance
    for nshape, seg in enumerate(world_map.countries):
        country_name = country_names[nshape]
        if country_name not in country_to_color:
            print('could not find: ' + country_name)
            continue
        color = country_to_color[country_name]
        poly = Polygon(seg, facecolor=color, edgecolor=color)
        axes.add_patch(poly)

    add_legend()
    fig = plt.gcf()
    fig.set_size_inches(30, 15)

    plt.axis('off')
    fig.savefig('countries_by_poverty_rate_world_bank_data.png', dpi=100,
                bbox_inches='tight',
                pad_inches=0)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def load_basemap_shapefile(shapefile, shapefile_reference_name, **kwargs):
    '''
    Creates a basemap, then loads a shapefile into it. If an axis
    is not specified, it creates a new figure with a single axis.

    Arguments:
    shapefile: The filename of the shapefile (as a string).
    shapefile_reference_name: The name to give the shapefile
        in Basemap.
    **kwargs: arguments to the Basemap constructor.

    Returns:
    bmap: The Basemap instance.
    shapefile_shapes: The shape reference of the shapefile in the
        Basemap instance.
    shapefile_info: The reference for the shape information in the
        shapefile Basemap instance.
    ax: The axis the Basemap instance was loaded into.
    '''
    #Basemap needs shapefile paths with no extensions, because
    #apparently doing this one-liner inside the code would be
    #too much effort.
    shapefile_no_ext = os.path.splitext(shapefile)[0]

    #Create an axis object if one wasn't passed as a kwarg:
    if 'ax' not in kwargs:
        ax = plt.figure().add_subplot(111)
        kwargs['ax'] = ax
    
    bmap = Basemap(**kwargs)
    bmap.readshapefile(shapefile_no_ext,shapefile_reference_name)
    shapefile_shapes = getattr(bmap,shapefile_reference_name)
    shapefile_info = getattr(bmap,shapefile_reference_name+'_info')
    return bmap, shapefile_shapes, shapefile_info, kwargs['ax']

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def shapefile_boundary(boundary_type, region_names):
    '''
    :param boundary_type: The type of spatial subset boundary
    :type boundary_type: :mod:'string'

    :param region_names: An array of regions for spatial subset
    :type region_names: :mod:'list'
    '''
    # Read the shapefile
    map_read = Basemap()
    regions = []
    shapefile_dir = os.sep.join([os.path.dirname(__file__), 'shape'])
    map_read.readshapefile(os.path.join(shapefile_dir, boundary_type),
                           boundary_type, drawbounds=False)
    # Note: The shapefile may contain countries or states with latin letters.
    # Hence, the output of readshapefile can be a mix of ascii and unicode
    # strings. Because python 2 and 3 treat them differently, we must
    # explicitly check.
    if boundary_type == 'us_states':
        for region_name in region_names:
            region_name = _force_unicode(region_name)
            for iregion, region_info in enumerate(map_read.us_states_info):
                state1 = _force_unicode(region_info['st'], 'latin-1')
                state2 = _force_unicode(region_info['state'], 'latin-1')
                if state1 == region_name or state2 == region_name:
                    regions.append(np.array(map_read.us_states[iregion]))
    elif boundary_type == 'countries':
        for region_name in region_names:
            region_name = _force_unicode(region_name)
            for iregion, region_info in enumerate(map_read.countries_info):
                country = _force_unicode(region_info['COUNTRY'], 'latin-1')
                if (country.replace(" ", "").lower() ==
                    region_name.replace(" ", "").lower()):
                    regions.append(np.array(map_read.countries[iregion]))
    return regions

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
def create_map():
    """
    create the base map for the choropleth.
    """
    logging.debug('create_map() -- Please wait while I create the world.')
    degrees_width = 118.0
    degrees_height = 55.0
    center_lat = 44.5
    center_lon = -110.0
    my_map = Basemap(  # ax=ax,
        projection='merc',  # default is cyl
        ellps='WGS84',
        lat_0=center_lat, lon_0=center_lon,
        llcrnrlat=center_lat - degrees_height / 2.0,
        llcrnrlon=center_lon - degrees_width / 2.0,
        urcrnrlat=center_lat + degrees_height / 2.0,
        urcrnrlon=center_lon + degrees_width / 2.0,
        resolution='i',  # 'c', 'l', 'i', 'h', 'f'
    )
    logging.debug('created map')
    logging.debug('loading shapes...')
    for section_name in CONTEST_SECTIONS.keys():
        # logging.debug('trying to load shape for %s', section_name)
        try:
            my_map.readshapefile('shapes/%s' % section_name, section_name, drawbounds=False)
        except IOError, err:
            logging.error('Could not load shape for %s' % section_name)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
class Mapper:
    def __init__(self):
        # Great Bay 43.0669N, 70.8686W
        self.m = Basemap(llcrnrlon=-70.96,llcrnrlat=43.,urcrnrlon=-70.55,urcrnrlat=43.2, projection='merc', resolution ='f')
        self.fig = plt.figure()
        self.ax = self.fig.add_subplot(111)
        self.m.drawmapboundary(fill_color='#ffffff')
        self.m.readshapefile('./n130', 'n130')
        self.m.fillcontinents(color='#cccccc',lake_color='#ffffff')

        self.patches = []
        for info, shape in zip(self.m.n130_info, self.m.n130):
            self.patches.append(Polygon(np.array(shape), True))
        self.ax.add_collection(PatchCollection(self.patches, facecolor= '#ffffff', edgecolor='#ffffff', linewidths=0., zorder=2))

    def plot(self, pointA, pointB):
        self.m.plot([pointA[1], pointB[1]], [pointA[0], pointB[0]], latlon=True, linewidth=1, color='#666666')

    def trace(self, points):
        for i, value in enumerate(points):
            try:
                self.plot(points[i], points[i+1])
            except (IndexError, ValueError) as e:
                print(e)
        self.m.plot(points[0][1], points[0][0], latlon=True, label='', marker='*', linewidth=.5, color='k')
        self.m.plot(points[i][1], points[i][0], latlon=True, label='', marker='8', linewidth=.5, color='k')

    def display(self, name):
        # draw coastlines, meridians and parallels.
        plt.savefig('.'.join([name, 'pdf']), format='pdf')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import readshapefile [as 别名]
class Mapper:
    def __init__(self, state):
        with open(state_boundries_file_loc) as f:
            for line in f:
                line = line.rstrip().split(',')
                if line[0].lower() == state.lower():
                    ur_lon, ll_lat = float(line[1]), float(line[3])
                    ll_lon, ur_lat = float(line[2]), float(line[4])
                    center_lon = ll_lon + ur_lon
                    if math.copysign(1, ll_lon) != math.copysign(1, ur_lon):
                        center_lon = -180 + center_lon
                    else:
                        center_lon = center_lon / 2
                    center_lat = (ll_lat + ur_lat) / 2

                    self.m = Basemap(llcrnrlon = ll_lon, llcrnrlat = ll_lat,
                                urcrnrlon = ur_lon, urcrnrlat = ur_lat,
                                projection = 'lcc', lat_1 = center_lat,
                                lon_0 = center_lon)
    def setBorders(self):
        self.m.readshapefile(states_shapefile_loc,
                            'states', drawbounds = True)

    def plotPoint(self, lat, lon):
        self.m.plot(lat, lon, latlon=True, marker='o', color='m')

    def displayMap(self):
        plt.show()