本文整理汇总了Python中mpl_toolkits.basemap.Basemap.is_land方法的典型用法代码示例。如果您正苦于以下问题:Python Basemap.is_land方法的具体用法?Python Basemap.is_land怎么用?Python Basemap.is_land使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类mpl_toolkits.basemap.Basemap
的用法示例。
在下文中一共展示了Basemap.is_land方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: do_calc
# 需要导入模块: 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())
示例2: cleanContinents
# 需要导入模块: 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
示例3: propagate_towards_land
# 需要导入模块: 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
示例4: latBucketing
# 需要导入模块: 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
示例5: cleanContinents
# 需要导入模块: 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
示例6: bucketOceans
# 需要导入模块: 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
示例7: get_ocean_mask
# 需要导入模块: 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)
示例8: print
# 需要导入模块: 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()
示例9: get_fvcom
# 需要导入模块: 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:
#.........这里部分代码省略.........
示例10: zip
# 需要导入模块: 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 '########################################'
示例11: __init__
# 需要导入模块: 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)
示例12: Basemap
# 需要导入模块: 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)
示例13: main
# 需要导入模块: 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)
示例14: zip
# 需要导入模块: 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 '########################################'
示例15: isOcean
# 需要导入模块: 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])