Python Basemap.is_land方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def do_calc(LATLIMS_AM, LONLIMS_AM, indir, outdir):
    land_checker = Basemap()
    if land_checker.is_land(LATLIMS_AM, LONLIMS_AM):
        print "SOS! Sorry you have selected a land pixel!"
        pygame.mixer.music.load("SOS.midi")
        pygame.mixer.music.play()
        while pygame.mixer.music.get_busy():
            # plot animado?
            time.sleep(1)
    else:
        dataAM = extract_series(LATLIMS_AM, LONLIMS_AM, indir)
        data_am = np.double(dataAM["Series"])
        if all(np.isnan(a) for a in data_am):
            print "THE SOUND OF SILENCE. Also, BATMAN. Everything is Rest and NaN"
            pygame.mixer.music.load("Batman_song.midi")
            pygame.mixer.music.play()
            while pygame.mixer.music.get_busy():
                # Anim plot? See Matplotlib.Animation
                time.sleep(1)
        else:
            am = get_music(data_am)

            music = pygame.mixer.Sound("Oc.midi")
            pygame.mixer.music.load("Oc.midi")
            pygame.mixer.music.play()
            anim = plot_animation(
                data_am,
                (u"Music from Lat = %.2f Lon = %.2f" % (dataAM["Lat"], dataAM["Lon"])),
                "serie.png",
                t_max=36000,
            )  # music.get_length())

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def cleanContinents(m, xpoints, ypoints):

    """
	Takes in the projection m and the grid of x and y
	clears out the points that are inside continents

	returns a matrix where each column corresponds to a value in xpoints
	"""

    seaMap = [[j for j in ypoints if not Basemap.is_land(m, m(i, j)[0], m(i, j)[1])] for i in xpoints]

    return seaMap

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def propagate_towards_land(proj, res, xmin, xmax, ymin, ymax, lon, lat):
		rcParams['figure.subplot.hspace'] = 0.1 # less height between subplots
		fig = plt.figure()
		# Initialise map
		startTime = datetime.now()
		map = Basemap(llcrnrlon=xmin,llcrnrlat=ymin,
					  urcrnrlon=xmax, urcrnrlat=ymax, 
					  resolution=res, projection=proj)
		ax = fig.add_subplot(211)
		map.drawcoastlines()
		map.drawmapboundary(fill_color='aqua')
		map.fillcontinents(color='coral',lake_color='aqua')
		initTime = datetime.now()-startTime

		# Propagate single particle until hitting coast
		lonHistory = np.array([lon])
		latHistory = np.array([lat])
		startTime = datetime.now()
		while not map.is_land(lon, lat):
			lon = lon + 0.01
			lat = lat + 0.001
			lonHistory = np.append(lonHistory, lon)
			latHistory = np.append(latHistory, lat)
		calcTime = datetime.now()-startTime

		# Plot and save figure
		plt.text(xmin, ymax, 'initialisation time: ' + 
						str(initTime.total_seconds()) + 's\n'
					'calculation time: ' + 
						str(calcTime.total_seconds()) + 's\n'
					'iterations: ' + str(len(latHistory)) +
					',  resolution: ' + res)
		map.plot(lonHistory, latHistory, '.b')
		map.plot(lon, lat, '*y')
		filename = 'figures_test_landmask_map_size/%s-%s-%s-%s-%s-%s.png' \
				% (proj, res, xmin, xmax, ymin, ymax)

		# Fullres zoom around landing point
		ax = fig.add_subplot(212)
		subsetmap = Basemap(llcrnrlon=lon-0.5,llcrnrlat=lat-.5,
					  urcrnrlon=lon+0.5, urcrnrlat=lat+.5, 
					  resolution='f', projection=proj)
		subsetmap.drawcoastlines()
		subsetmap.drawmapboundary(fill_color='aqua')
		subsetmap.fillcontinents(color='coral',lake_color='aqua')
		subsetmap.plot(lonHistory, latHistory, '*b')
		subsetmap.plot(lon, lat, '.y')
		
		#plt.show()
		plt.savefig(filename)
		print 'Saved figure: ' + filename

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def latBucketing(rawData,dx):

	buckets = {-90 + i*dx:[] for i in range(180/dx+1)}
	for bucket in buckets:
		for i in range(len(rawData['lon'])):
			x = rawData['lon'][i]
			y = rawData['lat'][i]
			if bucket + dx > y >= bucket:
				if not Basemap.is_land(m,m(x,y)[0],m(x,y)[1]):
					buckets[bucket].append(rawData['classif'][i])

	rates = ratesInBuckets(buckets)

	return rates

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def cleanContinents(m,xpoints,ypoints,predictions):

	cleanPredictions = []
	seaMap = []
	for i in range(len(xpoints)):
		xi = xpoints[i]
		pred = []
		y = []
		for j in range(len(ypoints)):
			yj = ypoints[j]
			if not Basemap.is_land(m,m(xi,yj)[0],m(xi,yj)[1]):
				pred.append(predictions[i][j])
				y.append(yj)
		cleanPredictions.append(pred)
		seaMap.append(y)

	return cleanPredictions,seaMap

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def bucketOceans(m,rawData):

	oceanDict = {ocean:{} for ocean in oceans}

	for i in range(len(rawData['lat'])):
		x = rawData['lon'][i]
		y = rawData['lat'][i]
		classif = rawData['classif'][i]
		ocean = oceanBoundaries(x,y)
		if not Basemap.is_land(m,m(x,y)[0],m(x,y)[1]):
			if ocean in oceans:
				if classif not in oceanDict[ocean]:
					oceanDict[ocean][classif] = 1
				else:
					oceanDict[ocean][classif] += 1
		print float(i)/len(rawData['lat'])
	print 'Bucketing done with a bucket of len '+str(len(oceanDict))
	return oceanDict

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def get_ocean_mask():
    print('Getting ocean mask...')
    from mpl_toolkits.basemap import Basemap
    latmin = grd[1].min()
    latmax = grd[1].max()
    lonmin = grd[0].min()
    lonmax = grd[0].max()
    print("Latitude {}--{},\n\
Longitude {}--{}".format(
    round(latmin,2),
    round(latmax,2),
    round(lonmin,2),
    round(lonmax,2)))
    m = Basemap(rsphere=6378137,resolution=coastres,projection='cea',
        llcrnrlat=latmin,urcrnrlat=latmax,
        llcrnrlon=lonmin,urcrnrlon=lonmax)
    (east,north) = m(grd[0],grd[1])

    ocean_mask = [not m.is_land(x,y) for (x,y) in zip(east,north)]#list(map(lambda x,y: not m.is_land(x,y),zip(x,y)))
    return np.array(ocean_mask)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
    print("computing adjacency...", file=sys.stderr, flush=True)
    adjacency = pd.DataFrame(False, index=regions, columns=regions)
    regions_set = set(regions)
    land_regions = set()
    covered = set()
    num_neighbors = args.num_neighbors
    for i in range(num_lats):
        for j in range(num_lngs):
            neighbors = [(i + lat_iterator, j + lng_iterator) for lat_iterator in
                         list(range(-num_neighbors, num_neighbors + 1)) for lng_iterator in
                         list(range(-num_neighbors, num_neighbors + 1))]
            for (n1, n2) in neighbors:
                if (n1, n2) in regions_set and \
                                (n1, n2) != (i, j) and \
                                n1 >= 0 and n2 >= 0 and \
                        m.is_land(country_lngs_m[n1, n2], country_lats_m[n1, n2]):
                    adjacency.ix[(i, j), (n1, n2)] = True
                    adjacency.ix[(n1, n2), (i, j)] = True
                    land_regions.add((n1, n2))

    land_regions = list(land_regions)
    print("done", file=sys.stderr, flush=True)
    print('%s land regions (out of %s)' % (len(land_regions), len(regions)), file=sys.stderr, flush=True)

else:
    regions = ['F', 'M']

region_name2int = dict([(name, i) for (i, name) in enumerate(regions)])
land_region_name2int = dict([(name, i) for (i, name) in enumerate(land_regions)])
review_frequency = Counter()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
class get_fvcom():
    def __init__(self, mod):
        self.modelname = mod
    
    def get_url(self, starttime, trackdays):
        '''
        get different url according to starttime and endtime.
        urls are monthly.
        '''
        #self.hours = int(round((endtime-starttime).total_seconds()/60/60))
        endtime = starttime + timedelta(trackdays)
        self.starttime = starttime; self.endtime = endtime
        self.trackdays = trackdays
        #int(round((endtime-starttime).total_seconds()/60/60/24))
        #print self.hours
                
        if self.modelname == "global":
            turl = '''http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_GLOBAL_FORECAST.nc'''
            
            try:
                self.tdata = netCDF4.Dataset(turl).variables
                #MTime = tdata['time'] 
                MTime = self.tdata['Times']
            except:
                print '"massbay" database is unavailable!'
                raise Exception()
            #Times = netCDF4.num2date(MTime[:],MTime.units)
            Times = []
            for i in MTime:
                strt = '201'+i[3]+'-'+i[5]+i[6]+'-'+i[8]+i[9]+' '+i[11]+i[12]+':'+i[14]+i[15]
                Times.append(datetime.strptime(strt,'%Y-%m-%d %H:%M'))#'''
            fmodtime = Times[0]; emodtime = Times[-1]+timedelta(days=1)    
            if starttime<fmodtime or starttime>emodtime or endtime<fmodtime or endtime>emodtime:
                #print 'Time: Error! Model(global) only works between %s with %s(UTC).'%(fmodtime,emodtime)
                #raise Exception()
                url = 'error' 
                return url,fmodtime,emodtime
            '''
            npTimes = np.array(Times)
            tm1 = npTimes-starttime; #tm2 = mtime-t2
            index1 = np.argmin(abs(tm1))#
            #index1 = netCDF4.date2index(starttime,MTime,select='nearest')
            index2 = index1 + self.days#
            #print 'index1,index2',index1,index2
            #url = url.format(index1, index2)
            self.mTime = Times[index1:index2] #'''
            self.mTime = np.array(Times) 
            #self.url = turl
            #loncs = self.tdata['lonc'][:]; self.latc = self.tdata['latc'][:]  #quantity:165095
            # Basic model data.
            self.basicdata = np.load('/var/www/cgi-bin/ioos/track/FVCOM_global_basic_data.npz')            
            self.lonc, self.latc = self.basicdata['lonc'], self.basicdata['latc']  #quantity:165095
            self.lons, self.lats = self.basicdata['lon'], self.basicdata['lat']
            self.h = self.basicdata['h']; self.siglay = self.basicdata['siglay']#; print '3'
            
            # Real-time model data.
            self.iy = [i for i in range(len(self.mTime)) if self.mTime[i].day==self.starttime.day]
            
            self.mtime = np.array(Times[self.iy[0]:self.iy[0]+int(trackdays)+1])
            
        
        return turl,fmodtime,emodtime
        
    def get_globaltrack(self,lon,lat,depth,track_way): #,b_index,nvdepth,,bcon 
        '''
        Get forecast points start at lon,lat
        '''
        
        modpts = dict(lon=[lon], lat=[lat])
        # For boundary.
        self.gmap=Basemap(projection='cyl',llcrnrlat=lat-self.trackdays, urcrnrlat=lat+self.trackdays, llcrnrlon=lon-self.trackdays, urcrnrlon=lon+self.trackdays,resolution='l')
        
        nodeindex = np.argwhere((self.lons >= lon-0.3) & (self.lons <= lon+0.3) & (self.lats >= lat-0.3) & (self.lats <= lat+0.3))        
        if len(nodeindex)==0:
            print 'No model data around. Out of Model area.'
            return modpts 
        waterdepth = self.h[[nodeindex]]; #print waterdepth
        if np.mean(waterdepth)<(abs(depth)): 
            print 'This point is too shallow.Less than %d meter.'%abs(depth)
            return modpts #raise Exception()
        depth_total = self.siglay[:,nodeindex[0]]*np.mean(waterdepth); #print depth_total
        layer = np.argmin(abs(depth_total+depth)); #print 'layer',layer
        
        u = self.tdata['u'][self.iy[0]:self.iy[0]+int(self.trackdays)+1,layer,:]; #print '5'
        v = self.tdata['v'][self.iy[0]:self.iy[0]+int(self.trackdays)+1,layer,:]; 
        #modpts = [(lon,lat)]#;st = []
        
        ld = self.trackdays*12
        for j in xrange(int(ld)): 
            if self.gmap.is_land(lat,lon):
                return modpts
            inds = np.argwhere((self.lonc >= lon-0.3) & (self.lonc <= lon+0.3) & (self.latc >= lat-0.3) & (self.latc <= lat+0.3))
            if len(inds) == 0:
                return modpts#print 'inds',inds

            if track_way=='backward' : # backwards case
                tratime = self.endtime-timedelta(hours=j*2)
                iday = [i for i in range(len(self.mtime)) if self.mtime[i].day==tratime.day]; #print iday
                u_t1 = np.mean(u[iday[0]][inds])*(-1); v_t1 = np.mean(v[iday[0]][inds])*(-1)
            else:
#.........这里部分代码省略.........

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
		print '\t%6.2f seconds to initialise map' % \
				(datetime.now()-startTime).total_seconds()

		# Extract polygons for faster checking of stranding
		polys = [p.boundary for p in map.landpolygons]
		
		# Check beaching of random points within map bounds
		npoints = 10000
		isLand = [0]*npoints
		np.random.seed(1)
		x = np.random.uniform(xmin, xmax, npoints)
		y = np.random.uniform(ymin, ymax, npoints)

		# using Basemap.is_land()
		startTime = datetime.now()
		isLand = [map.is_land(X, Y) for X,Y in zip(x,y)]
		print '\t%6.2f seconds to check that %s points out of %s are ' \
				'stranded, using Basemap.is_land()' % \
				((datetime.now()-startTime).total_seconds(), \
				sum(isLand), npoints)

		# using polygons
		startTime = datetime.now()
		isLand = points_in_polys(np.c_[x, y], polys)
		print '\t%6.2f seconds to check that %s points out of %s are ' \
				'stranded, using matplotlib polygons' % \
				((datetime.now()-startTime).total_seconds(), \
				sum(isLand), npoints)

	print 
print '########################################'

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
class Grib:
  def __init__(self, file = 'gribs/20171128_151424_.grb'):
    self.file = file
    self.gribfile = pygrib.open(file)
    
  def plot(self, idx = 0, show = False):
    grb = self.gribfile.select()[idx]
    grbU = self.gribfile.select(name='10 metre U wind component')[idx]
    grbV = self.gribfile.select(name='10 metre V wind component')[idx]
    lons = np.linspace(float(grb['longitudeOfFirstGridPointInDegrees']), \
              float(grb['longitudeOfLastGridPointInDegrees']), int(grb['Ni']) )
    lats = np.linspace(float(grb['latitudeOfFirstGridPointInDegrees']), \
                float(grb['latitudeOfLastGridPointInDegrees']), int(grb['Nj']) )
    data = np.sqrt(np.square(grbU.values)+np.square(grbV.values))
    grid_lon, grid_lat = np.meshgrid(lons, lats) #regularly spaced 2D grid
    self.map = Basemap(projection='cyl', llcrnrlon=lons.min(), \
        urcrnrlon=lons.max(),llcrnrlat=lats.min(),urcrnrlat=lats.max(), \
        resolution='c')
     
    x, y = self.map(grid_lon, grid_lat)
     
    cs =self.map.pcolormesh(x,y,data,shading='flat',cmap=plt.cm.rainbow)
    
    self.map.drawcoastlines()
    self.map.drawmapboundary()
    self.map.drawparallels(np.arange(-90.,120.,30.),labels=[1,0,0,0])
    self.map.drawmeridians(np.arange(-180.,180.,60.),labels=[0,0,0,1])
     
    plt.colorbar(cs,orientation='vertical', shrink=0.5)
    plt.title('Predicted wind strength') 
    plt.savefig('grib.png')
    if show:
        plt.show()
    plt.close()


  def interpolator(self):
    grb = self.gribfile.select()[0]
    grb1 = self.gribfile.select()[1]
    grb_l = self.gribfile.select()[-1]
    lons = np.linspace(float(grb['longitudeOfFirstGridPointInDegrees']), \
              float(grb['longitudeOfLastGridPointInDegrees']), int(grb['Ni']) )
    lats = np.linspace(float(grb['latitudeOfFirstGridPointInDegrees']), \
                float(grb['latitudeOfLastGridPointInDegrees']), int(grb['Nj']) )
    grbU = self.gribfile.select(name='10 metre U wind component')
    grbV = self.gribfile.select(name='10 metre V wind component')
    times = np.array([i for i,g in enumerate(grbU) if g.values.shape == grb.values.shape])
    U = np.dstack([g.values for g in grbU if g.values.shape == grb.values.shape])
    V = np.dstack([g.values for g in grbV if g.values.shape == grb.values.shape])
    print(lats.shape, lons.shape, times.shape, U.shape)
    if lats[0] > lats[-1]:
        lats = np.flip(lats,0)
        U = np.flip(U, 0)
        V = np.flip(V,0)
    self.interpolatorU = RegularGridInterpolator((lats, lons, times), U)
    self.interpolatorV = RegularGridInterpolator((lats, lons, times), V)

  def getwind(self, p, time):
    U = self.interpolatorU((p.lat, p.lon, time))
    V = self.interpolatorV((p.lat, p.lon, time))
    wind = sqrt(U**2 + V**2)
    winddir = atan2(U,V)*180/pi
    return wind, winddir
  
  def is_land(self,p):
    return self.map.is_land(p.lat, p.lon)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap(projection='aeqd', lon_0 = 10, lat_0 = 50, resolution='h')

x, y = map(0, 0)

print map.is_land(x, y)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def main():
    usage = """
    python create_grid.py
    """
    num_args= 0
    parser = OptionParser(usage=usage)

    parser.add_option('', '--center-x', dest='center_x', default=None,
            help="Center the grid on x (longitude) coordinate", type=float)
    parser.add_option('', '--center-y', dest='center_y', default=None,
            help="Center the grid on y (latitude) coordinate", type=float)
    parser.add_option('', '--grid_width', dest='grid_width', default=60,
            help='How wide the grid is in degrees', type=float)

    parser.add_option('-r', '--resolution', dest='resolution', default=4, 
            help="The resolution we want the grid to be at", type='int')
    parser.add_option('-o', '--output-file', dest='output_file', default=None,
            help="A file to dump the output to", type='str')
    parser.add_option('', '--min-x', dest='min_x', default=None,
            help='The minimum longitude', type='float')
    parser.add_option('', '--max-x', dest='max_x', default=None,
            help='The maximum longitude', type='float')
    parser.add_option('', '--min-y', dest='min_y', default=None,
            help='The minimum latitude', type='float')
    parser.add_option('', '--max-y', dest='max_y', default=None,
            help='The maximum latitude', type='float')

    #parser.add_option('-o', '--options', dest='some_option', default='yo', help="Place holder for a real option", type='str')
    #parser.add_option('-u', '--useless', dest='uselesss', default=False, action='store_true', help='Another useless option')

    (options, args) = parser.parse_args()


    if options.center_x is not None:
        options.min_x = options.center_x - options.grid_width / 2;
        options.max_x = options.center_x + options.grid_width / 2;

    if options.center_y is not None:
        options.min_y = options.center_y - options.grid_width / 2;
        options.max_y = options.center_y + options.grid_width / 2;

    if options.max_y > 90:
        options.max_y = 90 
    if options.min_y < -90:
        options.min_y = 90

    xs = np.linspace(options.min_x, options.max_x, options.resolution)
    ys = np.linspace(options.min_y, options.max_y, options.resolution)

    from mpl_toolkits.basemap import Basemap

    my_map = Basemap(llcrnrlon=options.min_x,llcrnrlat=options.min_y,
                  urcrnrlon=options.max_x,urcrnrlat=options.max_y,\
                rsphere=(6378137.00,6356752.3142),\
                resolution='l',area_thresh=1000.,projection='merc',\
                lat_1=50.,lon_0=-107.)# draw coastlines, country boundaries, fill continents.

    is_land = dict()

    for x, y in it.product(xs, ys):
        mx, my = my_map(x, y)

        if my_map.is_land(mx, my):
            print x, y

    (options, args) = parser.parse_args()

    if len(args) < num_args:
        parser.print_help()
        sys.exit(1)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
		print '\t%6.2f seconds to initialise map' % \
				(datetime.now()-startTime).total_seconds()

		# Extract polygons for faster checking of stranding
		polys = [p.boundary for p in map.landpolygons]
		
		# Check beaching of random points within map bounds
		npoints = 1000
		isLand = [0]*npoints
		np.random.seed(1)
		x = np.random.uniform(xmin, xmax, npoints)
		y = np.random.uniform(ymin, ymax, npoints)

		# using Basemap.is_land()
		startTime = datetime.now()
		isLandBasemap = [map.is_land(X, Y) for X,Y in zip(x,y)]
		print '\t%6.2f seconds to check that %s points out of %s are ' \
				'stranded, using Basemap.is_land()' % \
				((datetime.now()-startTime).total_seconds(), \
				sum(isLandBasemap), npoints)

		# using polygons
		startTime = datetime.now()
		isLandPoly = points_in_polys(np.c_[x, y], polys)
		print '\t%6.2f seconds to check that %s points out of %s are ' \
				'stranded, using matplotlib polygons' % \
				((datetime.now()-startTime).total_seconds(), \
				sum(isLandPoly), npoints)

	print 
print '########################################'

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import is_land [as 别名]
def isOcean(x,y):

	return not Basemap.is_land(m,m(x,y)[0],m(x,y)[1])