def numpy_band_stats(ds, tif, band_no, partitions=100):
band = ds.GetRasterBand(band_no)
data = band.ReadAsArray()
no_data_val = band.GetNoDataValue()
data_type = get_datatype(band)
mask_data = ma.masked_where(data == no_data_val, data)
if data_type is 'Categorical':
no_categories = np.max(mask_data) - np.min(mask_data) + 1
else:
no_categories = np.nan
image_source = geoio.RasterioImageSource(tif)
l = [basename(tif), band_no, no_data_val,
ds.RasterYSize, ds.RasterXSize,
np.min(mask_data), np.max(mask_data),
np.mean(mask_data), np.std(mask_data),
data_type, float(no_categories),
image_nans(image_source, partitions)]
ds = None
return [str(a) for a in l]
python类masked_where()的实例源码
def test_mem_masked_where(self,level=rlevel):
# Ticket #62
from numpy.ma import masked_where, MaskType
a = np.zeros((1, 1))
b = np.zeros(a.shape, MaskType)
c = masked_where(b, a)
a-c
def test_mem_masked_where(self,level=rlevel):
# Ticket #62
from numpy.ma import masked_where, MaskType
a = np.zeros((1, 1))
b = np.zeros(a.shape, MaskType)
c = masked_where(b, a)
a-c
def test_mem_masked_where(self,level=rlevel):
# Ticket #62
from numpy.ma import masked_where, MaskType
a = np.zeros((1, 1))
b = np.zeros(a.shape, MaskType)
c = masked_where(b, a)
a-c
test_regression.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 25
收藏 0
点赞 0
评论 0
def test_mem_masked_where(self,level=rlevel):
# Ticket #62
from numpy.ma import masked_where, MaskType
a = np.zeros((1, 1))
b = np.zeros(a.shape, MaskType)
c = masked_where(b, a)
a-c
def test_mem_masked_where(self,level=rlevel):
# Ticket #62
from numpy.ma import masked_where, MaskType
a = np.zeros((1, 1))
b = np.zeros(a.shape, MaskType)
c = masked_where(b, a)
a-c
def plot_distribution(ax, lon_samples, lat_samples, to_plot='d', resolution=30, **kwargs):
if 'cmap' in kwargs:
cmap = kwargs.pop('cmap')
else:
cmap = next(cmaps)
artists = []
if 'd' in to_plot:
lon_grid, lat_grid, density = density_distribution(
lon_samples, lat_samples, resolution)
density = ma.masked_where(density <= 0.05*density.max(), density)
a = ax.pcolormesh(lon_grid, lat_grid, density, cmap=cmap,
transform=ccrs.PlateCarree(), **kwargs)
artists.append(a)
if 'e' in to_plot:
lon_grid, lat_grid, cumulative_density = cumulative_density_distribution(
lon_samples, lat_samples, resolution)
a = ax.contour(lon_grid, lat_grid, cumulative_density, levels=[
0.683, 0.955], cmap=cmap, transform=ccrs.PlateCarree())
artists.append(a)
if 's' in to_plot:
a = ax.scatter(lon_samples, lat_samples, color=cmap(
[0., 0.5, 1.])[-1], alpha=0.1, transform=ccrs.PlateCarree(), edgecolors=None, **kwargs)
artists.append(a)
return artists
def test_mem_masked_where(self,level=rlevel):
# Ticket #62
from numpy.ma import masked_where, MaskType
a = np.zeros((1, 1))
b = np.zeros(a.shape, MaskType)
c = masked_where(b, a)
a-c
def test_mem_masked_where(self,level=rlevel):
# Ticket #62
from numpy.ma import masked_where, MaskType
a = np.zeros((1, 1))
b = np.zeros(a.shape, MaskType)
c = masked_where(b, a)
a-c
def import_all_year_data(tile):
temp = np.zeros([46, 2400*2400], np.dtype(float))
if int(tile[5:6])<2:
temp[:]=np.nan
for doy in range(1, 369, 8):
evifile = buildVrtFile(root, doy, tile, 'evi')
cloudfile = buildVrtFile(root, doy, tile, 'cloudmask')
aerosolfile = buildVrtFile(root, doy, tile, 'aerosolmask')
#if no file found for this DOY
if evifile == 0: continue
#doyList.append(doy)
#build vrt for EVI
vrtEVI = os.path.join(os.path.dirname(evifile), str(1000+doy)[1:]+tile+'EVI_vrt.vrt')
print "Building the vrt file: ", evifile
os.system('gdalbuildvrt -separate -input_file_list '+evifile+' '+vrtEVI)
inEVI = gdal.Open(vrtEVI)
EVI = inEVI.ReadAsArray()
#build vrt for cloudmask
vrtcloud = os.path.join(os.path.dirname(cloudfile), str(1000+doy)[1:]+tile+'cloud_vrt.vrt')
print "Building the vrt file: ", cloudfile
os.system('gdalbuildvrt -separate -input_file_list '+cloudfile+' '+vrtcloud)
incloud = gdal.Open(vrtcloud)
cloud = incloud.ReadAsArray()
#build vrt for aerosol
vrtaerosol = os.path.join(os.path.dirname(aerosolfile), str(1000+doy)[1:]+tile+'aerosol_vrt.vrt')
print "Building the vrt file: ", aerosolfile
os.system('gdalbuildvrt -separate -input_file_list '+aerosolfile+' '+vrtaerosol)
inaerosol = gdal.Open(vrtaerosol)
aerosol = inaerosol.ReadAsArray()
global rows, cols, geoProj, geoTran
rows = 2400
cols = 2400
geoTran = inEVI.GetGeoTransform()
geoProj = inEVI.GetProjection()
#mask for bad quality
EVIgood = ma.masked_where((cloud != 1)|(aerosol == 0)|(EVI < 0)|(EVI > 10000), EVI)
EVIgood = EVIgood.reshape(EVIgood.size/2400/2400, 2400*2400)
medianEVI = np.nanmedian(EVIgood, axis=0)
EVI = None
aerosol = None
cloud = None
EVIgood = None
#assign to the 46 layer of matrix
temp[(doy-1)/8, :] = medianEVI
meanEVI = None
return temp
def barbs(self, x, y, u, v, *args, **kwargs):
"""
Make a wind barb plot (u, v) with on the map.
(see matplotlib.pyplot.barbs documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \*args and \**kwargs passed on to matplotlib.pyplot.barbs
Returns two matplotlib.axes.Barbs instances, one for the Northern
Hemisphere and one for the Southern Hemisphere.
"""
if _matplotlib_version < '0.98.3':
msg = dedent("""
barb method requires matplotlib 0.98.3 or higher,
you have %s""" % _matplotlib_version)
raise NotImplementedError(msg)
ax, plt = self._ax_plt_from_kw(kwargs)
# allow callers to override the hold state by passing hold=True|False
b = ax.ishold()
h = kwargs.pop('hold',None)
if h is not None:
ax.hold(h)
lons, lats = self(x, y, inverse=True)
unh = ma.masked_where(lats <= 0, u)
vnh = ma.masked_where(lats <= 0, v)
ush = ma.masked_where(lats > 0, u)
vsh = ma.masked_where(lats > 0, v)
try:
retnh = ax.barbs(x,y,unh,vnh,*args,**kwargs)
kwargs['flip_barb']=True
retsh = ax.barbs(x,y,ush,vsh,*args,**kwargs)
except:
ax.hold(b)
raise
ax.hold(b)
# Because there are two collections returned in general,
# we can't set the current image...
#if plt is not None and ret.get_array() is not None:
# plt.sci(retnh)
# clip for round polar plots.
if self.round:
retnh,c = self._clipcircle(ax,retnh)
retsh,c = self._clipcircle(ax,retsh)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
return retnh,retsh
def barbs(self, x, y, u, v, *args, **kwargs):
"""
Make a wind barb plot (u, v) with on the map.
(see matplotlib.pyplot.barbs documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \*args and \**kwargs passed on to matplotlib.pyplot.barbs
Returns two matplotlib.axes.Barbs instances, one for the Northern
Hemisphere and one for the Southern Hemisphere.
"""
if _matplotlib_version < '0.98.3':
msg = dedent("""
barb method requires matplotlib 0.98.3 or higher,
you have %s""" % _matplotlib_version)
raise NotImplementedError(msg)
ax, plt = self._ax_plt_from_kw(kwargs)
# allow callers to override the hold state by passing hold=True|False
b = ax.ishold()
h = kwargs.pop('hold',None)
if h is not None:
ax.hold(h)
lons, lats = self(x, y, inverse=True)
unh = ma.masked_where(lats <= 0, u)
vnh = ma.masked_where(lats <= 0, v)
ush = ma.masked_where(lats > 0, u)
vsh = ma.masked_where(lats > 0, v)
try:
retnh = ax.barbs(x,y,unh,vnh,*args,**kwargs)
kwargs['flip_barb']=True
retsh = ax.barbs(x,y,ush,vsh,*args,**kwargs)
except:
ax.hold(b)
raise
ax.hold(b)
# Because there are two collections returned in general,
# we can't set the current image...
#if plt is not None and ret.get_array() is not None:
# plt.sci(retnh)
# clip for round polar plots.
if self.round:
retnh,c = self._clipcircle(ax,retnh)
retsh,c = self._clipcircle(ax,retsh)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
return retnh,retsh
def get_vowel_segments(media_path, n_fft=2048):
downsample = 1
samplerate = 44100 // downsample
win_s = n_fft // downsample # fft size
hop_s = n_fft // downsample # hop size
s = source(media_path, samplerate, hop_s)
samplerate = s.samplerate
tolerance = 0.6
pitch_o = pitch("yin", win_s, hop_s, samplerate)
pitch_o.set_unit("Hz")
pitch_o.set_tolerance(tolerance)
pitches = []
confidences = []
# total number of frames read
total_frames = 0
samples=[]
pitches=[]
while True:
samples, read = s()
pitch_ = pitch_o(samples)[0]
#pitch = int(round(pitch))
confidence = pitch_o.get_confidence()
#print("%f %f %f" % (total_frames / float(samplerate), pitch, confidence))
pitches += [pitch_]
confidences += [confidence]
total_frames += read
if read < hop_s: break
pitches = np.array(pitches)
confidences = np.array(confidences)
cleaned_pitches = ma.masked_where(confidences < tolerance, pitches)
cleaned_pitches = ma.masked_where(cleaned_pitches > 1000, cleaned_pitches)
try: output = list(np.logical_not(cleaned_pitches.mask))
except: output = []
return output