Here are the examples of the python api numpy.ma.filled taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.
56 Examples
3
Source : recfunctions.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def _fix_output(output, usemask=True, asrecarray=False):
"""
Private function: return a recarray, a ndarray, a MaskedArray
or a MaskedRecords depending on the input parameters
"""
if not isinstance(output, MaskedArray):
usemask = False
if usemask:
if asrecarray:
output = output.view(MaskedRecords)
else:
output = ma.filled(output)
if asrecarray:
output = output.view(recarray)
return output
def _fix_defaults(output, defaults=None):
3
Source : test_old_ma.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def test_testBasic1d(self):
# Test of basic array creation and properties in 1 dimension.
(x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
assert_(not isMaskedArray(x))
assert_(isMaskedArray(xm))
assert_equal(shape(xm), s)
assert_equal(xm.shape, s)
assert_equal(xm.dtype, x.dtype)
assert_equal(xm.size, reduce(lambda x, y:x * y, s))
assert_equal(count(xm), len(m1) - reduce(lambda x, y:x + y, m1))
assert_(eq(xm, xf))
assert_(eq(filled(xm, 1.e20), xf))
assert_(eq(x, xm))
def test_testBasic2d(self):
3
Source : test_old_ma.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def test_testMasked(self):
# Test of masked element
xx = arange(6)
xx[1] = masked
assert_(str(masked) == '--')
assert_(xx[1] is masked)
assert_equal(filled(xx[1], 0), 0)
def test_testAverage1(self):
3
Source : mstats_basic.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def msign(x):
"""Returns the sign of x, or 0 if x is masked."""
return ma.filled(np.sign(x), 0)
def pearsonr(x,y):
3
Source : NESOSIM.py
with GNU General Public License v3.0
from akpetty
with GNU General Public License v3.0
from akpetty
def fillMaskAndNaNWithZero(arr):
""" Helper function: Fill masked and nan values in an array
with 0, in place
Args:
arr (var): A numpy ndarray
returns:
None (performs operation in place)
"""
arr[np.isnan(arr)] = 0.
arr = ma.filled(arr, 0.)
arr[~np.isfinite(arr)] = 0.
def fill_nan_no_negative(arr, region_maskG, negative_to_zero=True):
3
Source : transforms.py
with GNU General Public License v3.0
from Artikash
with GNU General Public License v3.0
from Artikash
def get_points(self):
if self._invalid:
points = self._transform.transform(self._bbox.get_points())
points = np.ma.filled(points, 0.0)
self._points = points
self._invalid = 0
return self._points
get_points.__doc__ = Bbox.get_points.__doc__
3
Source : fmm_planner.py
with MIT License
from devendrachaplot
with MIT License
from devendrachaplot
def set_goal(self, goal, auto_improve=False):
traversible_ma = ma.masked_values(self.traversible * 1, 0)
goal_x, goal_y = int(goal[0] / (self.scale * 1.)), \
int(goal[1] / (self.scale * 1.))
if self.traversible[goal_x, goal_y] == 0. and auto_improve:
goal_x, goal_y = self._find_nearest_goal([goal_x, goal_y])
traversible_ma[goal_x, goal_y] = 0
dd = skfmm.distance(traversible_ma, dx=1)
dd = ma.filled(dd, np.max(dd) + 1)
self.fmm_dist = dd
return
def set_multi_goal(self, goal_map):
3
Source : fmm_planner.py
with MIT License
from devendrachaplot
with MIT License
from devendrachaplot
def set_multi_goal(self, goal_map):
traversible_ma = ma.masked_values(self.traversible * 1, 0)
traversible_ma[goal_map == 1] = 0
dd = skfmm.distance(traversible_ma, dx=1)
dd = ma.filled(dd, np.max(dd) + 1)
self.fmm_dist = dd
return
def get_short_term_goal(self, state):
3
Source : grid.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def actnum_array(self):
"""Returns the 3D ndarray which for active cells.
Values are 1 for active, 0 for inactive, in C order (read only).
"""
actnumv = self.get_actnum().values
actnumv = ma.filled(actnumv, fill_value=0)
return actnumv
@property
3
Source : _well_roxapi.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def _get_roxlog(wlogtypes, wlogrecords, roxlrun, lname): # pragma: no cover
roxcurve = roxlrun.log_curves[lname]
tmplog = roxcurve.get_values().astype(np.float64)
tmplog = npma.filled(tmplog, fill_value=np.nan)
tmplog[tmplog == -999] = np.nan
if roxcurve.is_discrete:
wlogtypes[lname] = "DISC"
wlogrecords[lname] = roxcurve.get_code_names()
else:
wlogtypes[lname] = "CONT"
wlogrecords[lname] = None
return tmplog
def export_well_roxapi(
3
Source : mstats_basic.py
with MIT License
from ktraunmueller
with MIT License
from ktraunmueller
def ttest_ind(a, b, axis=0):
a, b, axis = _chk2_asarray(a, b, axis)
(x1, x2) = (a.mean(axis), b.mean(axis))
(v1, v2) = (a.var(axis=axis, ddof=1), b.var(axis=axis, ddof=1))
(n1, n2) = (a.count(axis), b.count(axis))
df = n1+n2-2
svar = ((n1-1)*v1+(n2-1)*v2) / float(df)
svar == 0
t = (x1-x2)/ma.sqrt(svar*(1.0/n1 + 1.0/n2)) # N-D COMPUTATION HERE!!!!!!
t = ma.filled(t, 1) # replace NaN t-values with 1.0
probs = betai(0.5*df,0.5,float(df)/(df+t*t)).reshape(t.shape)
return t, probs.squeeze()
ttest_ind.__doc__ = stats.ttest_ind.__doc__
3
Source : mstats_basic.py
with BSD 3-Clause "New" or "Revised" License
from mspass-team
with BSD 3-Clause "New" or "Revised" License
from mspass-team
def msign(x):
"""Returns the sign of x, or 0 if x is masked."""
return ma.filled(np.sign(x), 0)
def pearsonr(x, y):
3
Source : pfp_io.py
with BSD 3-Clause "New" or "Revised" License
from OzFlux
with BSD 3-Clause "New" or "Revised" License
from OzFlux
def smap_docarbonfluxes(cf,ds,smap_label,si,ei):
ncname = cf["Variables"][smap_label]["name"]
data,flag,attr = pfp_utils.GetSeriesasMA(ds,ncname,si=si,ei=ei)
data = data*12.01*1800/1E6
data = numpy.ma.filled(data,float(-9999))
return data,flag
def smap_donetshortwave(ds,smap_label,si,ei):
3
Source : pfp_io.py
with BSD 3-Clause "New" or "Revised" License
from OzFlux
with BSD 3-Clause "New" or "Revised" License
from OzFlux
def smap_donetshortwave(ds,smap_label,si,ei):
ts = int(ds.globalattributes["time_step"])
# do the net shortwave radiation
Fsd,Fsd_flag,a = pfp_utils.GetSeriesasMA(ds,"Fsd",si=si,ei=ei)
Fsu,Fsu_flag,a = pfp_utils.GetSeriesasMA(ds,"Fsu",si=si,ei=ei)
# get the net shortwave radiation and convert to MJ/m2/day at the same time
Fnsw = ((Fsd - Fsu)*ts*60)/1E6
# now get the QC flag
Fnsw_flag = Fsd_flag+Fsu_flag
Fnsw = numpy.ma.filled(Fnsw,float(-9999))
return Fnsw,Fnsw_flag
def smap_dopressure(ds,smap_label,si,ei):
3
Source : pfp_io.py
with BSD 3-Clause "New" or "Revised" License
from OzFlux
with BSD 3-Clause "New" or "Revised" License
from OzFlux
def smap_dopressure(ds,smap_label,si,ei):
ps,ps_flag,attr = pfp_utils.GetSeriesasMA(ds,"ps",si=si,ei=ei)
ps = ps/float(1000)
ps = numpy.ma.filled(ps,float(-9999))
return ps,ps_flag
def smap_doshortwave(ds,smap_label,si,ei):
3
Source : pfp_io.py
with BSD 3-Clause "New" or "Revised" License
from OzFlux
with BSD 3-Clause "New" or "Revised" License
from OzFlux
def smap_doshortwave(ds,smap_label,si,ei):
ts = int(ds.globalattributes["time_step"])
Fsd,Fsd_flag,a = pfp_utils.GetSeriesasMA(ds,"Fsd",si=si,ei=ei)
Fsd = (Fsd*ts*60)/1E6
Fsd = numpy.ma.filled(Fsd,float(-9999))
return Fsd,Fsd_flag
def smap_parseformat(fmt):
3
Source : test_old_ma.py
with Apache License 2.0
from pierreant
with Apache License 2.0
from pierreant
def test_testBasic1d(self):
# Test of basic array creation and properties in 1 dimension.
(x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
self.assertFalse(isMaskedArray(x))
self.assertTrue(isMaskedArray(xm))
self.assertEqual(shape(xm), s)
self.assertEqual(xm.shape, s)
self.assertEqual(xm.dtype, x.dtype)
self.assertEqual(xm.size, reduce(lambda x, y:x * y, s))
self.assertEqual(count(xm), len(m1) - reduce(lambda x, y:x + y, m1))
self.assertTrue(eq(xm, xf))
self.assertTrue(eq(filled(xm, 1.e20), xf))
self.assertTrue(eq(x, xm))
def test_testBasic2d(self):
3
Source : test_old_ma.py
with Apache License 2.0
from pierreant
with Apache License 2.0
from pierreant
def test_testMasked(self):
# Test of masked element
xx = arange(6)
xx[1] = masked
self.assertTrue(str(masked) == '--')
self.assertTrue(xx[1] is masked)
self.assertEqual(filled(xx[1], 0), 0)
def test_testAverage1(self):
0
Source : mrecords.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def __setattr__(self, attr, val):
"""
Sets the attribute attr to the value val.
"""
# Should we call __setmask__ first ?
if attr in ['mask', 'fieldmask']:
self.__setmask__(val)
return
# Create a shortcut (so that we don't have to call getattr all the time)
_localdict = object.__getattribute__(self, '__dict__')
# Check whether we're creating a new field
newattr = attr not in _localdict
try:
# Is attr a generic attribute ?
ret = object.__setattr__(self, attr, val)
except Exception:
# Not a generic attribute: exit if it's not a valid field
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
optinfo = ndarray.__getattribute__(self, '_optinfo') or {}
if not (attr in fielddict or attr in optinfo):
exctype, value = sys.exc_info()[:2]
raise exctype(value)
else:
# Get the list of names
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
# Check the attribute
if attr not in fielddict:
return ret
if newattr:
# We just added this one or this setattr worked on an
# internal attribute.
try:
object.__delattr__(self, attr)
except Exception:
return ret
# Let's try to set the field
try:
res = fielddict[attr][:2]
except (TypeError, KeyError):
raise AttributeError("record array has no attribute %s" % attr)
if val is masked:
_fill_value = _localdict['_fill_value']
if _fill_value is not None:
dval = _localdict['_fill_value'][attr]
else:
dval = val
mval = True
else:
dval = filled(val)
mval = getmaskarray(val)
obj = ndarray.__getattribute__(self, '_data').setfield(dval, *res)
_localdict['_mask'].__setitem__(attr, mval)
return obj
def __getitem__(self, indx):
0
Source : test_old_ma.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def test_testBasic2d(self):
# Test of basic array creation and properties in 2 dimensions.
for s in [(4, 3), (6, 2)]:
(x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
x.shape = s
y.shape = s
xm.shape = s
ym.shape = s
xf.shape = s
assert_(not isMaskedArray(x))
assert_(isMaskedArray(xm))
assert_equal(shape(xm), s)
assert_equal(xm.shape, s)
assert_equal(xm.size, reduce(lambda x, y:x * y, s))
assert_equal(count(xm),
len(m1) - reduce(lambda x, y:x + y, m1))
assert_(eq(xm, xf))
assert_(eq(filled(xm, 1.e20), xf))
assert_(eq(x, xm))
self.setup()
def test_testArithmetic(self):
0
Source : test_old_ma.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def test_testAddSumProd(self):
# Test add, sum, product.
(x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
assert_(eq(np.add.reduce(x), add.reduce(x)))
assert_(eq(np.add.accumulate(x), add.accumulate(x)))
assert_(eq(4, sum(array(4), axis=0)))
assert_(eq(4, sum(array(4), axis=0)))
assert_(eq(np.sum(x, axis=0), sum(x, axis=0)))
assert_(eq(np.sum(filled(xm, 0), axis=0), sum(xm, axis=0)))
assert_(eq(np.sum(x, 0), sum(x, 0)))
assert_(eq(np.product(x, axis=0), product(x, axis=0)))
assert_(eq(np.product(x, 0), product(x, 0)))
assert_(eq(np.product(filled(xm, 1), axis=0),
product(xm, axis=0)))
if len(s) > 1:
assert_(eq(np.concatenate((x, y), 1),
concatenate((xm, ym), 1)))
assert_(eq(np.add.reduce(x, 1), add.reduce(x, 1)))
assert_(eq(np.sum(x, 1), sum(x, 1)))
assert_(eq(np.product(x, 1), product(x, 1)))
def test_testCI(self):
0
Source : test_old_ma.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def test_testCopySize(self):
# Tests of some subtle points of copying and sizing.
n = [0, 0, 1, 0, 0]
m = make_mask(n)
m2 = make_mask(m)
assert_(m is m2)
m3 = make_mask(m, copy=1)
assert_(m is not m3)
x1 = np.arange(5)
y1 = array(x1, mask=m)
assert_(y1._data is not x1)
assert_(allequal(x1, y1._data))
assert_(y1.mask is m)
y1a = array(y1, copy=0)
assert_(y1a.mask is y1.mask)
y2 = array(x1, mask=m3, copy=0)
assert_(y2.mask is m3)
assert_(y2[2] is masked)
y2[2] = 9
assert_(y2[2] is not masked)
assert_(y2.mask is m3)
assert_(allequal(y2.mask, 0))
y2a = array(x1, mask=m, copy=1)
assert_(y2a.mask is not m)
assert_(y2a[2] is masked)
y2a[2] = 9
assert_(y2a[2] is not masked)
assert_(y2a.mask is not m)
assert_(allequal(y2a.mask, 0))
y3 = array(x1 * 1.0, mask=m)
assert_(filled(y3).dtype is (x1 * 1.0).dtype)
x4 = arange(4)
x4[2] = masked
y4 = resize(x4, (8,))
assert_(eq(concatenate([x4, x4]), y4))
assert_(eq(getmask(y4), [0, 0, 1, 0, 0, 0, 1, 0]))
y5 = repeat(x4, (2, 2, 2, 2), axis=0)
assert_(eq(y5, [0, 0, 1, 1, 2, 2, 3, 3]))
y6 = repeat(x4, 2, axis=0)
assert_(eq(y5, y6))
def test_testPut(self):
0
Source : test_old_ma.py
with GNU General Public License v3.0
from adityaprakash-bobby
with GNU General Public License v3.0
from adityaprakash-bobby
def test_testTakeTransposeInnerOuter(self):
# Test of take, transpose, inner, outer products
x = arange(24)
y = np.arange(24)
x[5:6] = masked
x = x.reshape(2, 3, 4)
y = y.reshape(2, 3, 4)
assert_(eq(np.transpose(y, (2, 0, 1)), transpose(x, (2, 0, 1))))
assert_(eq(np.take(y, (2, 0, 1), 1), take(x, (2, 0, 1), 1)))
assert_(eq(np.inner(filled(x, 0), filled(y, 0)),
inner(x, y)))
assert_(eq(np.outer(filled(x, 0), filled(y, 0)),
outer(x, y)))
y = array(['abc', 1, 'def', 2, 3], object)
y[2] = masked
t = take(y, [0, 3, 4])
assert_(t[0] == 'abc')
assert_(t[1] == 2)
assert_(t[2] == 3)
def test_testInplace(self):
0
Source : NESOSIM.py
with GNU General Public License v3.0
from akpetty
with GNU General Public License v3.0
from akpetty
def loadData(yearT, dayT, precipVar, windVar, concVar, driftVar, dxStr, extraStr):
""" Load daily forcings
Temp added transpose to convert grid to proper row/column index.
"""
dayStr='%03d' %dayT
#------- Read in precipitation -----------
try:
print('Loading gridded snowfall forcing from:', forcingPath+'Precip/'+precipVar+'/'+str(yearT)+'/'+precipVar+'sf'+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr)
precipDayG=np.load(forcingPath+'Precip/'+precipVar+'/'+str(yearT)+'/'+precipVar+'sf'+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr, allow_pickle=True)
except:
if (dayStr=='365'):
precipDayG=np.load(forcingPath+'Precip/'+precipVar+'/sf/'+str(yearT)+'/'+precipVar+'sf'+dxStr+'-'+str(yearT)+'_d'+'364'+extraStr, allow_pickle=True)
print('no leap year data, used data from the previous day')
else:
print('No precip data so exiting!')
exit()
#------- Read in wind magnitude -----------
try:
print('Loading gridded wind forcing from:', forcingPath+'Winds/'+windVar+'/'+str(yearT)+'/'+windVar+'winds'+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr)
windDayG=np.load(forcingPath+'Winds/'+windVar+'/'+str(yearT)+'/'+windVar+'winds'+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr, allow_pickle=True)
except:
if (dayStr=='365'):
print('no leap year data, using data from the previous day')
windDayG=np.load(forcingPath+'Winds/'+windVar+'/'+str(yearT)+'/'+windVar+'winds'+dxStr+'-'+str(yearT)+'_d'+'364'+extraStr, allow_pickle=True)
else:
print('No wind data so exiting!')
exit()
#------- Read in ice concentration -----------
try:
print('Loading gridded ice conc forcing from:', forcingPath+'IceConc/'+concVar+'/'+str(yearT)+'/iceConcG_'+concVar+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr)
iceConcDayG=np.load(forcingPath+'IceConc/'+concVar+'/'+str(yearT)+'/iceConcG_'+concVar+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr, allow_pickle=True)
except:
if (dayStr=='365'):
print('no leap year data, using data from the previous day')
iceConcDayG=np.load(forcingPath+'IceConc/'+concVar+'/'+str(yearT)+'/iceConcG_'+concVar+dxStr+'-'+str(yearT)+'_d'+'364'+extraStr, allow_pickle=True)
else:
print('No ice conc data so exiting!')
exit()
# fill with zero
iceConcDayG[~np.isfinite(iceConcDayG)]=0.
#------- Read in ice drifts -----------
try:
print('Loading gridded ice drift forcing from:', forcingPath+'IceDrift/'+driftVar+'/'+str(yearT)+'/'+driftVar+'_driftG'+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr)
driftGdayG=np.load(forcingPath+'IceDrift/'+driftVar+'/'+str(yearT)+'/'+driftVar+'_driftG'+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr, allow_pickle=True)
except:
# if no drifts exist for that day then just set drifts to nan array (i.e. no drift).
print('No drift data')
driftGdayG = np.empty((2, iceConcDayG.shape[0], iceConcDayG.shape[1]))
driftGdayG[:] = np.nan
driftGdayG = ma.filled(driftGdayG, np.nan)
#print(driftGdayG)
#------- Read in temps (not currently used, placeholder) -----------
try:
tempDayG=np.load(forcingPath+'Temp/'+precipVar+'/t2m/'+str(yearT)+'/t2m'+dxStr+'-'+str(yearT)+'_d'+dayStr+extraStr, allow_pickle=True)
except:
# if no drifts exist for that day then just set drifts to masked array (i.e. no drift).
#print('No temp data')
tempDayG = np.empty((iceConcDayG.shape[0], iceConcDayG.shape[1]))
tempDayG[:] = np.nan
#tempDayG=ma.masked_all((iceConcDayG.shape[0], iceConcDayG.shape[1]))
return iceConcDayG, precipDayG, driftGdayG, windDayG, tempDayG
def densityCalc(snowDepthsT, iceConcDayT, region_maskT):
0
Source : transforms.py
with MIT License
from alvarobartt
with MIT License
from alvarobartt
def get_points(self):
if self._invalid:
p = self._bbox.get_points()
# Transform all four points, then make a new bounding box
# from the result, taking care to make the orientation the
# same.
points = self._transform.transform(
[[p[0, 0], p[0, 1]],
[p[1, 0], p[0, 1]],
[p[0, 0], p[1, 1]],
[p[1, 0], p[1, 1]]])
points = np.ma.filled(points, 0.0)
xs = min(points[:, 0]), max(points[:, 0])
if p[0, 0] > p[1, 0]:
xs = xs[::-1]
ys = min(points[:, 1]), max(points[:, 1])
if p[0, 1] > p[1, 1]:
ys = ys[::-1]
self._points = np.array([
[xs[0], ys[0]],
[xs[1], ys[1]]
])
self._invalid = 0
return self._points
get_points.__doc__ = Bbox.get_points.__doc__
0
Source : test_old_ma.py
with MIT License
from alvarobartt
with MIT License
from alvarobartt
def test_testCopySize(self):
# Tests of some subtle points of copying and sizing.
n = [0, 0, 1, 0, 0]
m = make_mask(n)
m2 = make_mask(m)
assert_(m is m2)
m3 = make_mask(m, copy=1)
assert_(m is not m3)
x1 = np.arange(5)
y1 = array(x1, mask=m)
assert_(y1._data is not x1)
assert_(allequal(x1, y1._data))
assert_(y1.mask is m)
y1a = array(y1, copy=0)
# For copy=False, one might expect that the array would just
# passed on, i.e., that it would be "is" instead of "==".
# See gh-4043 for discussion.
assert_(y1a._mask.__array_interface__ ==
y1._mask.__array_interface__)
y2 = array(x1, mask=m3, copy=0)
assert_(y2.mask is m3)
assert_(y2[2] is masked)
y2[2] = 9
assert_(y2[2] is not masked)
assert_(y2.mask is m3)
assert_(allequal(y2.mask, 0))
y2a = array(x1, mask=m, copy=1)
assert_(y2a.mask is not m)
assert_(y2a[2] is masked)
y2a[2] = 9
assert_(y2a[2] is not masked)
assert_(y2a.mask is not m)
assert_(allequal(y2a.mask, 0))
y3 = array(x1 * 1.0, mask=m)
assert_(filled(y3).dtype is (x1 * 1.0).dtype)
x4 = arange(4)
x4[2] = masked
y4 = resize(x4, (8,))
assert_(eq(concatenate([x4, x4]), y4))
assert_(eq(getmask(y4), [0, 0, 1, 0, 0, 0, 1, 0]))
y5 = repeat(x4, (2, 2, 2, 2), axis=0)
assert_(eq(y5, [0, 0, 1, 1, 2, 2, 3, 3]))
y6 = repeat(x4, 2, axis=0)
assert_(eq(y5, y6))
def test_testPut(self):
0
Source : mrecords.py
with GNU General Public License v3.0
from Artikash
with GNU General Public License v3.0
from Artikash
def __setattr__(self, attr, val):
"Sets the attribute attr to the value val."
# Should we call __setmask__ first ?
if attr in ['mask', 'fieldmask']:
self.__setmask__(val)
return
# Create a shortcut (so that we don't have to call getattr all the time)
_localdict = object.__getattribute__(self, '__dict__')
# Check whether we're creating a new field
newattr = attr not in _localdict
try:
# Is attr a generic attribute ?
ret = object.__setattr__(self, attr, val)
except:
# Not a generic attribute: exit if it's not a valid field
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
optinfo = ndarray.__getattribute__(self, '_optinfo') or {}
if not (attr in fielddict or attr in optinfo):
exctype, value = sys.exc_info()[:2]
raise exctype(value)
else:
# Get the list of names ......
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
# Check the attribute
if attr not in fielddict:
return ret
if newattr: # We just added this one
try: # or this setattr worked on an internal
# attribute.
object.__delattr__(self, attr)
except:
return ret
# Let's try to set the field
try:
res = fielddict[attr][:2]
except (TypeError, KeyError):
raise AttributeError("record array has no attribute %s" % attr)
#
if val is masked:
_fill_value = _localdict['_fill_value']
if _fill_value is not None:
dval = _localdict['_fill_value'][attr]
else:
dval = val
mval = True
else:
dval = filled(val)
mval = getmaskarray(val)
obj = ndarray.__getattribute__(self, '_data').setfield(dval, *res)
_localdict['_mask'].__setitem__(attr, mval)
return obj
def __getitem__(self, indx):
0
Source : stats_mstats_short.py
with MIT License
from birforce
with MIT License
from birforce
def quantiles(a, prob=list([.25,.5,.75]), alphap=.4, betap=.4, axis=None,
limit=(), masknan=False):
"""
Computes empirical quantiles for a data array.
Samples quantile are defined by :math:`Q(p) = (1-g).x[i] +g.x[i+1]`,
where :math:`x[j]` is the *j*th order statistic, and
`i = (floor(n*p+m))`, `m=alpha+p*(1-alpha-beta)` and `g = n*p + m - i`.
Typical values of (alpha,beta) are:
- (0,1) : *p(k) = k/n* : linear interpolation of cdf (R, type 4)
- (.5,.5) : *p(k) = (k+1/2.)/n* : piecewise linear
function (R, type 5)
- (0,0) : *p(k) = k/(n+1)* : (R type 6)
- (1,1) : *p(k) = (k-1)/(n-1)*. In this case, p(k) = mode[F(x[k])].
That's R default (R type 7)
- (1/3,1/3): *p(k) = (k-1/3)/(n+1/3)*. Then p(k) ~ median[F(x[k])].
The resulting quantile estimates are approximately median-unbiased
regardless of the distribution of x. (R type 8)
- (3/8,3/8): *p(k) = (k-3/8)/(n+1/4)*. Blom.
The resulting quantile estimates are approximately unbiased
if x is normally distributed (R type 9)
- (.4,.4) : approximately quantile unbiased (Cunnane)
- (.35,.35): APL, used with PWM ?? JP
- (0.35, 0.65): PWM ?? JP p(k) = (k-0.35)/n
Parameters
----------
a : array-like
Input data, as a sequence or array of dimension at most 2.
prob : array-like, optional
List of quantiles to compute.
alpha : float, optional
Plotting positions parameter, default is 0.4.
beta : float, optional
Plotting positions parameter, default is 0.4.
axis : int, optional
Axis along which to perform the trimming.
If None (default), the input array is first flattened.
limit : tuple
Tuple of (lower, upper) values.
Values of `a` outside this closed interval are ignored.
Returns
-------
quants : MaskedArray
An array containing the calculated quantiles.
Examples
--------
>>> from scipy.stats.mstats import mquantiles
>>> a = np.array([6., 47., 49., 15., 42., 41., 7., 39., 43., 40., 36.])
>>> mquantiles(a)
array([ 19.2, 40. , 42.8])
Using a 2D array, specifying axis and limit.
>>> data = np.array([[ 6., 7., 1.],
[ 47., 15., 2.],
[ 49., 36., 3.],
[ 15., 39., 4.],
[ 42., 40., -999.],
[ 41., 41., -999.],
[ 7., -999., -999.],
[ 39., -999., -999.],
[ 43., -999., -999.],
[ 40., -999., -999.],
[ 36., -999., -999.]])
>>> mquantiles(data, axis=0, limit=(0, 50))
array([[ 19.2 , 14.6 , 1.45],
[ 40. , 37.5 , 2.5 ],
[ 42.8 , 40.05, 3.55]])
>>> data[:, 2] = -999.
>>> mquantiles(data, axis=0, limit=(0, 50))
masked_array(data =
[[19.2 14.6 --]
[40.0 37.5 --]
[42.8 40.05 --]],
mask =
[[False False True]
[False False True]
[False False True]],
fill_value = 1e+20)
"""
if isinstance(a, np.ma.MaskedArray):
return stats.mstats.mquantiles(a, prob=prob, alphap=alphap, betap=alphap, axis=axis,
limit=limit)
if limit:
marr = stats.mstats.mquantiles(a, prob=prob, alphap=alphap, betap=alphap, axis=axis,
limit=limit)
return ma.filled(marr, fill_value=np.nan)
if masknan:
nanmask = np.isnan(a)
if nanmask.any():
marr = ma.array(a, mask=nanmask)
marr = stats.mstats.mquantiles(marr, prob=prob, alphap=alphap, betap=alphap,
axis=axis, limit=limit)
return ma.filled(marr, fill_value=np.nan)
# Initialization & checks ---------
data = np.asarray(a)
p = np.array(prob, copy=False, ndmin=1)
m = alphap + p*(1.-alphap-betap)
isrolled = False
#from _quantiles1d
if (axis is None):
data = data.ravel() #reshape(-1,1)
axis = 0
else:
axis = np.arange(data.ndim)[axis]
data = np.rollaxis(data, axis)
isrolled = True # keep track, maybe can be removed
x = np.sort(data, axis=0)
n = x.shape[0]
returnshape = list(data.shape)
returnshape[axis] = p
#TODO: check these
if n == 0:
return np.empty(len(p), dtype=float)
elif n == 1:
return np.resize(x, p.shape)
aleph = (n*p + m)
k = np.floor(aleph.clip(1, n-1)).astype(int)
ind = [None]*x.ndim
ind[0] = slice(None)
gamma = (aleph-k).clip(0,1)[ind]
q = (1.-gamma)*x[k-1] + gamma*x[k]
if isrolled:
return np.rollaxis(q, 0, axis+1)
else:
return q
def scoreatpercentile(data, per, limit=(), alphap=.4, betap=.4, axis=0, masknan=None):
0
Source : transforms.py
with MIT License
from buds-lab
with MIT License
from buds-lab
def get_points(self):
# docstring inherited
if self._invalid:
p = self._bbox.get_points()
# Transform all four points, then make a new bounding box
# from the result, taking care to make the orientation the
# same.
points = self._transform.transform(
[[p[0, 0], p[0, 1]],
[p[1, 0], p[0, 1]],
[p[0, 0], p[1, 1]],
[p[1, 0], p[1, 1]]])
points = np.ma.filled(points, 0.0)
xs = min(points[:, 0]), max(points[:, 0])
if p[0, 0] > p[1, 0]:
xs = xs[::-1]
ys = min(points[:, 1]), max(points[:, 1])
if p[0, 1] > p[1, 1]:
ys = ys[::-1]
self._points = np.array([
[xs[0], ys[0]],
[xs[1], ys[1]]
])
self._invalid = 0
return self._points
if DEBUG:
0
Source : mrecords.py
with Apache License 2.0
from dashanji
with Apache License 2.0
from dashanji
def __setattr__(self, attr, val):
"""
Sets the attribute attr to the value val.
"""
# Should we call __setmask__ first ?
if attr in ['mask', 'fieldmask']:
self.__setmask__(val)
return
# Create a shortcut (so that we don't have to call getattr all the time)
_localdict = object.__getattribute__(self, '__dict__')
# Check whether we're creating a new field
newattr = attr not in _localdict
try:
# Is attr a generic attribute ?
ret = object.__setattr__(self, attr, val)
except Exception:
# Not a generic attribute: exit if it's not a valid field
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
optinfo = ndarray.__getattribute__(self, '_optinfo') or {}
if not (attr in fielddict or attr in optinfo):
raise
else:
# Get the list of names
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
# Check the attribute
if attr not in fielddict:
return ret
if newattr:
# We just added this one or this setattr worked on an
# internal attribute.
try:
object.__delattr__(self, attr)
except Exception:
return ret
# Let's try to set the field
try:
res = fielddict[attr][:2]
except (TypeError, KeyError):
raise AttributeError("record array has no attribute %s" % attr)
if val is masked:
_fill_value = _localdict['_fill_value']
if _fill_value is not None:
dval = _localdict['_fill_value'][attr]
else:
dval = val
mval = True
else:
dval = filled(val)
mval = getmaskarray(val)
obj = ndarray.__getattribute__(self, '_data').setfield(dval, *res)
_localdict['_mask'].__setitem__(attr, mval)
return obj
def __getitem__(self, indx):
0
Source : test_old_ma.py
with Apache License 2.0
from dashanji
with Apache License 2.0
from dashanji
def test_testCopySize(self):
# Tests of some subtle points of copying and sizing.
n = [0, 0, 1, 0, 0]
m = make_mask(n)
m2 = make_mask(m)
assert_(m is m2)
m3 = make_mask(m, copy=True)
assert_(m is not m3)
x1 = np.arange(5)
y1 = array(x1, mask=m)
assert_(y1._data is not x1)
assert_(allequal(x1, y1._data))
assert_(y1._mask is m)
y1a = array(y1, copy=0)
# For copy=False, one might expect that the array would just
# passed on, i.e., that it would be "is" instead of "==".
# See gh-4043 for discussion.
assert_(y1a._mask.__array_interface__ ==
y1._mask.__array_interface__)
y2 = array(x1, mask=m3, copy=0)
assert_(y2._mask is m3)
assert_(y2[2] is masked)
y2[2] = 9
assert_(y2[2] is not masked)
assert_(y2._mask is m3)
assert_(allequal(y2.mask, 0))
y2a = array(x1, mask=m, copy=1)
assert_(y2a._mask is not m)
assert_(y2a[2] is masked)
y2a[2] = 9
assert_(y2a[2] is not masked)
assert_(y2a._mask is not m)
assert_(allequal(y2a.mask, 0))
y3 = array(x1 * 1.0, mask=m)
assert_(filled(y3).dtype is (x1 * 1.0).dtype)
x4 = arange(4)
x4[2] = masked
y4 = resize(x4, (8,))
assert_(eq(concatenate([x4, x4]), y4))
assert_(eq(getmask(y4), [0, 0, 1, 0, 0, 0, 1, 0]))
y5 = repeat(x4, (2, 2, 2, 2), axis=0)
assert_(eq(y5, [0, 0, 1, 1, 2, 2, 3, 3]))
y6 = repeat(x4, 2, axis=0)
assert_(eq(y5, y6))
def test_testPut(self):
0
Source : mrecords.py
with GNU General Public License v3.0
from dnn-security
with GNU General Public License v3.0
from dnn-security
def __setattr__(self, attr, val):
"""
Sets the attribute attr to the value val.
"""
# Should we call __setmask__ first ?
if attr in ['mask', 'fieldmask']:
self.__setmask__(val)
return
# Create a shortcut (so that we don't have to call getattr all the time)
_localdict = object.__getattribute__(self, '__dict__')
# Check whether we're creating a new field
newattr = attr not in _localdict
try:
# Is attr a generic attribute ?
ret = object.__setattr__(self, attr, val)
except Exception:
# Not a generic attribute: exit if it's not a valid field
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
optinfo = ndarray.__getattribute__(self, '_optinfo') or {}
if not (attr in fielddict or attr in optinfo):
raise
else:
# Get the list of names
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
# Check the attribute
if attr not in fielddict:
return ret
if newattr:
# We just added this one or this setattr worked on an
# internal attribute.
try:
object.__delattr__(self, attr)
except Exception:
return ret
# Let's try to set the field
try:
res = fielddict[attr][:2]
except (TypeError, KeyError) as e:
raise AttributeError(
f'record array has no attribute {attr}') from e
if val is masked:
_fill_value = _localdict['_fill_value']
if _fill_value is not None:
dval = _localdict['_fill_value'][attr]
else:
dval = val
mval = True
else:
dval = filled(val)
mval = getmaskarray(val)
obj = ndarray.__getattribute__(self, '_data').setfield(dval, *res)
_localdict['_mask'].__setitem__(attr, mval)
return obj
def __getitem__(self, indx):
0
Source : _gridprop_lowlevel.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def update_carray(self, undef=None, discrete=None, dtype=None, order="F"):
"""Copy (update) values from numpy to SWIG, 1D array, returns a pointer
to SWIG C array. If discrete is defined as True or False, force
the SWIG array to be of that kind.
Note that dtype will "override" current datatype if set. The resulting
carray will be in Fortran order, unless order is specified as 'C'
"""
dstatus = self._isdiscrete
if discrete is not None:
dstatus = bool(discrete)
if undef is None:
undef = xtgeo.UNDEF
if dstatus:
undef = xtgeo.UNDEF_INT
logger.debug("Entering conversion from numpy to C array ...")
values = self._values.copy()
if not dtype:
if dstatus:
values = values.astype(np.int32)
else:
values = values.astype(np.float64)
else:
values = values.astype(dtype)
values = ma.filled(values, undef)
if order == "F":
values = np.asfortranarray(values)
values1d = np.ravel(values, order="K")
if values1d.dtype == "float64" and dstatus and not dtype:
values1d = values1d.astype("int32")
logger.debug("Casting has been done")
if values1d.dtype == "float64":
logger.debug("Convert to carray (double)")
carray = _cxtgeo.new_doublearray(self.ntotal)
_cxtgeo.swig_numpy_to_carr_1d(values1d, carray)
elif values1d.dtype == "float32":
logger.debug("Convert to carray (float)")
carray = _cxtgeo.new_floatarray(self.ntotal)
_cxtgeo.swig_numpy_to_carr_f1d(values1d, carray)
elif values1d.dtype == "int32":
logger.debug("Convert to carray (int32)")
carray = _cxtgeo.new_intarray(self.ntotal)
_cxtgeo.swig_numpy_to_carr_i1d(values1d, carray)
else:
raise RuntimeError("Unsupported dtype, probable bug in {}".format(__name__))
return carray
def delete_carray(self, carray):
0
Source : regular_surface.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def get_values1d(
self, order="C", asmasked=False, fill_value=xtgeo.UNDEF, activeonly=False
):
"""Get a 1D numpy or masked array of the map values.
Args:
order (str): Flatteting is in C (default) or F order
asmasked (bool): If true, return as MaskedArray, other as standard
numpy ndarray with undef as np.nan or fill_value
fill_value (str): Relevent only if asmasked is False, this
will be the value of undef entries
activeonly (bool): If True, only active cells. Keys 'asmasked' and
'fill_value' are not revelant.
Returns:
A numpy 1D array or MaskedArray
"""
val = self.values.copy()
if order == "F":
val = ma.filled(val, fill_value=np.nan)
val = np.array(val, order="F")
val = ma.masked_invalid(val)
val = val.ravel(order="K")
if activeonly:
val = val[~val.mask]
if not asmasked and not activeonly:
val = ma.filled(val, fill_value=fill_value)
return val
def set_values1d(self, val, order="C"):
0
Source : _regsurf_export.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def export_petromod_binary(self, mfile, pmd_dataunits):
"""Export to petromod binary format."""
validunits = False
unitd = 15
unitz = 10
if isinstance(pmd_dataunits, tuple) and len(pmd_dataunits) == 2:
unitd, unitz = pmd_dataunits
if isinstance(unitd, int) and isinstance(unitz, int):
if unitd in PMD_DATAUNITDISTANCE.keys() and unitz in PMD_DATAUNITZ.keys():
validunits = True
if unitd < = 0 or unitz < = 0:
raise ValueError("Values for pmd_dataunits cannot be negative!")
if not validunits:
UserWarning(
"Format or values for pmd_dataunits out of range: Pair should be in ranges "
"{} and {}".format(PMD_DATAUNITDISTANCE, PMD_DATAUNITZ)
)
undef = 99999
dsc = "Content=Map,"
dsc += "DataUnitDistance={},".format(unitd)
dsc += "DataUnitZ={},".format(unitz)
dsc += "GridNoX={},".format(self.ncol)
dsc += "GridNoY={},".format(self.nrow)
dsc += "GridStepX={},".format(self.xinc)
dsc += "GridStepY={},".format(self.yinc)
dsc += "MapType=GridMap,"
dsc += "OriginX={},".format(self.xori)
dsc += "OriginY={},".format(self.yori)
dsc += "RotationAngle={},".format(self.rotation)
dsc += "RotationOriginX={},".format(self.xori)
dsc += "RotationOriginY={},".format(self.yori)
dsc += "Undefined={},".format(undef)
dsc += "Version=1.0"
values = np.ma.filled(self.values1d, fill_value=undef)
_cxtgeo.surf_export_petromod_bin(
mfile.get_cfhandle(),
dsc,
values,
)
mfile.cfclose()
def export_xtgregsurf(self, mfile):
0
Source : _regsurf_grid3d.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def from_grid3d(grid, template=None, where="top", mode="depth", rfactor=1):
"""Private function for deriving a surface from a 3D grid.
Note that rotated maps are currently not supported!
.. versionadded:: 2.1
"""
if where == "top":
klayer = 1
option = 0
elif where == "base":
klayer = grid.nlay
option = 1
else:
klayer, what = where.split("_")
klayer = int(klayer)
if grid.nlay < klayer < 0:
raise ValueError("Klayer out of range in where={}".format(where))
option = 0
if what == "base":
option = 1
if rfactor < 0.5:
raise KeyError("Refinefactor rfactor is too small, should be >= 0.5")
args = _update_regsurf(template, grid, rfactor=float(rfactor))
args["rotation"] = 0.0
# call C function to make a map
val = args["values"]
val = val.ravel(order="K")
val = ma.filled(val, fill_value=xtgeo.UNDEF)
svalues = val * 0.0 + xtgeo.UNDEF
ivalues = svalues.copy()
jvalues = svalues.copy()
grid._xtgformat1()
_cxtgeo.surf_sample_grd3d_lay(
grid.ncol,
grid.nrow,
grid.nlay,
grid._coordsv,
grid._zcornsv,
grid._actnumsv,
klayer,
args["ncol"],
args["nrow"],
args["xori"],
args["xinc"],
args["yori"],
args["yinc"],
args["rotation"],
svalues,
ivalues,
jvalues,
option,
)
logger.info("Extracted surfaces from 3D grid...")
svalues = np.ma.masked_greater(svalues, xtgeo.UNDEF_LIMIT)
ivalues = np.ma.masked_greater(ivalues, xtgeo.UNDEF_LIMIT)
jvalues = np.ma.masked_greater(jvalues, xtgeo.UNDEF_LIMIT)
if mode == "i":
ivalues = ivalues.reshape((args["ncol"], args["nrow"]))
ivalues = ma.masked_invalid(ivalues)
args["values"] = ivalues
return args, None, None
if mode == "j":
jvalues = jvalues.reshape((args["ncol"], args["nrow"]))
jvalues = ma.masked_invalid(jvalues)
args["values"] = jvalues
return args, None, None
svalues = svalues.reshape((args["ncol"], args["nrow"]))
svalues = ma.masked_invalid(svalues)
args["values"] = svalues
return args, ivalues, jvalues
def _update_regsurf(template, grid, rfactor=1.0):
0
Source : _regsurf_gridding.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def surf_fill(self, fill_value=None):
"""Replace the value of invalid 'data' cells (indicated by 'invalid')
by the value of the nearest valid data cell or a constant.
This is a quite fast method to fill undefined areas of the map.
The surface values are updated 'in-place'
.. versionadded:: 2.1
.. versionchanged:: 2.6 Added fill_value
"""
logger.info("Do fill...")
if fill_value is not None:
if np.isscalar(fill_value) and not isinstance(fill_value, str):
self.values = ma.filled(self.values, fill_value=float(fill_value))
else:
raise ValueError("Keyword fill_value must be int or float")
else:
invalid = ma.getmaskarray(self.values)
ind = scipy.ndimage.distance_transform_edt(
invalid, return_distances=False, return_indices=True
)
self._values = self._values[tuple(ind)]
logger.info("Do fill... DONE")
def smooth_median(self, iterations=1, width=1):
0
Source : _regsurf_oper.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def resample(self, other, mask=True):
"""Resample from other surface object to this surf."""
logger.info("Resampling...")
# a special case occur of the maps have same topology, but
# different masks
if self.compare_topology(other, strict=False):
self.values = other.values.copy()
return
svalues = np.ma.filled(self.values, fill_value=xtgeo.UNDEF)
ovalues = np.ma.filled(other.values, fill_value=xtgeo.UNDEF)
_cxtgeo.surf_resample(
other._ncol,
other._nrow,
other._xori,
other._xinc,
other._yori,
other._yinc,
other._yflip,
other._rotation,
ovalues,
self._ncol,
self._nrow,
self._xori,
self._xinc,
self._yori,
self._yinc,
self._yflip,
self._rotation,
svalues,
0 if not mask else 1,
)
self.values = np.ma.masked_greater(svalues, xtgeo.UNDEF_LIMIT)
self.set_values1d(svalues)
self._filesrc = "Resampled"
def distance_from_point(self, point=(0, 0), azimuth=0.0):
0
Source : _regsurf_oper.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def get_xy_values(self, order="C", asmasked=False):
"""Get X Y coordinate values as numpy 2D arrays."""
nno = self.ncol * self.nrow
ier, xvals, yvals = _cxtgeo.surf_xy_as_values(
self.xori,
self.xinc,
self.yori,
self.yinc * self.yflip,
self.ncol,
self.nrow,
self.rotation,
nno,
nno,
0,
)
if ier != 0:
raise XTGeoCLibError(f"Error in surf_xy_as_values, error code: {ier}")
# reshape
xvals = xvals.reshape((self.ncol, self.nrow))
yvals = yvals.reshape((self.ncol, self.nrow))
if order == "F":
xvals = np.array(xvals, order="F")
yvals = np.array(yvals, order="F")
if asmasked:
tmpv = ma.filled(self.values, fill_value=np.nan)
tmpv = np.array(tmpv, order=order)
tmpv = ma.masked_invalid(tmpv)
mymask = ma.getmaskarray(tmpv)
xvals = ma.array(xvals, mask=mymask, order=order)
yvals = ma.array(yvals, mask=mymask, order=order)
return xvals, yvals
def get_xy_values1d(self, order="C", activeonly=True):
0
Source : _regsurf_oper.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def operation_polygons(self, poly, value, opname="add", inside=True):
"""Operations restricted to polygons"""
if not isinstance(poly, Polygons):
raise ValueError("The poly input is not a Polygons instance")
if opname not in VALID_OPER_POLYS:
raise ValueError(f"Operation key opname has invalid value: {opname}")
# make a copy of the RegularSurface which is used a "filter" or "proxy"
# value will be 1 inside polygons, 0 outside. Undef cells are kept as is
proxy = self.copy()
proxy.values *= 0.0
vals = proxy.get_values1d(fill_value=xtgeo.UNDEF)
# value could be a scalar or another surface; if another surface,
# must ensure same topology
if isinstance(value, type(self)):
if not self.compare_topology(value):
raise ValueError("Input is RegularSurface, but not same map " "topology")
value = value.values.copy()
else:
# turn scalar value into numpy array
value = self.values.copy() * 0 + value
idgroups = poly.dataframe.groupby(poly.pname)
for _, grp in idgroups:
xcor = grp[poly.xname].values
ycor = grp[poly.yname].values
ier = _cxtgeo.surf_setval_poly(
proxy.xori,
proxy.xinc,
proxy.yori,
proxy.yinc,
proxy.ncol,
proxy.nrow,
proxy.yflip,
proxy.rotation,
vals,
xcor,
ycor,
1.0,
0,
)
if ier == -9:
xtg.warn("Polygon is not closed")
proxy.set_values1d(vals)
proxyv = proxy.values.astype(np.int8)
proxytarget = 1
if not inside:
proxytarget = 0
tmp = None
if opname == "add":
tmp = self.values.copy() + value
elif opname == "sub":
tmp = self.values.copy() - value
elif opname == "mul":
tmp = self.values.copy() * value
elif opname == "div":
# Dividing a map of zero is always a hazzle; try to obtain 0.0
# as result in these cases
if 0.0 in value:
xtg.warn(
"Dividing a surface with value or surface with zero "
"elements; may get unexpected results, try to "
"achieve zero values as result!"
)
with np.errstate(divide="ignore", invalid="ignore"):
this = ma.filled(self.values, fill_value=1.0)
that = ma.filled(value, fill_value=1.0)
mask = ma.getmaskarray(self.values)
tmp = np.true_divide(this, that)
tmp = np.where(np.isinf(tmp), 0, tmp)
tmp = np.nan_to_num(tmp)
tmp = ma.array(tmp, mask=mask)
elif opname == "set":
tmp = value
elif opname == "eli":
tmp = value * 0 + xtgeo.UNDEF
tmp = ma.masked_greater(tmp, xtgeo.UNDEF_LIMIT)
self.values[proxyv == proxytarget] = tmp[proxyv == proxytarget]
del tmp
0
Source : _blockedwell_roxapi.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def _roxapi_import_bwell(
self, rox, gname, bwname, wname, lognames, ijk, realisation
): # pylint: disable=too-many-statements
"""Private function for ROXAPI well import (get well from Roxar)"""
if gname in rox.project.grid_models:
gmodel = rox.project.grid_models[gname]
logger.info("RMS grid model < %s> OK", gname)
else:
raise ValueError("No such grid name present: {}".format(gname))
if bwname in gmodel.blocked_wells_set:
bwset = gmodel.blocked_wells_set[bwname]
logger.info("Blocked well set < %s> OK", bwname)
else:
raise ValueError("No such blocked well set: {}".format(bwname))
if wname in bwset.get_well_names():
self._wname = wname
else:
raise WellNotFoundError("No such well in blocked well set: {}".format(wname))
# pylint: disable=unnecessary-comprehension
bwprops = [item for item in bwset.properties]
bwnames = [item.name for item in bwset.properties]
bw_cellindices = bwset.get_cell_numbers()
dind = bwset.get_data_indices([wname])
cind = bw_cellindices[dind]
xyz = np.transpose(gmodel.get_grid(realisation=realisation).get_cell_centers(cind))
logs = OrderedDict()
logs["X_UTME"] = xyz[0].astype(np.float64)
logs["Y_UTMN"] = xyz[1].astype(np.float64)
logs["Z_TVDSS"] = xyz[2].astype(np.float64)
if ijk:
ijk = np.transpose(
gmodel.get_grid(realisation=realisation).grid_indexer.get_indices(cind)
)
logs["I_INDEX"] = ijk[0].astype(np.float64)
logs["J_INDEX"] = ijk[1].astype(np.float64)
logs["K_INDEX"] = ijk[2].astype(np.float64)
usenames = []
if lognames and lognames == "all":
usenames = bwnames
elif lognames:
usenames = lognames
for bwprop in bwprops:
lname = bwprop.name
if lname not in usenames:
continue
propvalues = bwprop.get_values(realisation=realisation)
tmplog = propvalues[dind].astype(np.float64)
tmplog = npma.filled(tmplog, fill_value=np.nan)
tmplog[tmplog == -999] = np.nan
if "discrete" in str(bwprop.type):
self._wlogtypes[lname] = "DISC"
self._wlogrecords[lname] = bwprop.code_names
else:
self._wlogtypes[lname] = "CONT"
self._wlogrecords[lname] = None
logs[lname] = tmplog
self._df = pd.DataFrame.from_dict(logs)
self._gname = gname
self._filesrc = None
# finally get some other metadata like RKB and topside X Y; as they
# seem to miss for the BW in RMS, try and get them from the
# well itself...
if wname in rox.project.wells:
self._rkb = rox.project.wells[wname].rkb
self._xpos, self._ypos = rox.project.wells[wname].wellhead
else:
self._rkb = None
self._xpos, self._ypos = self._df["X_UTME"][0], self._df["Y_UTMN"][0]
def _roxapi_export_bwell(
0
Source : test_etc_make_hcpv_maps.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def test_hcpvfz1(tmpdir, generate_plot):
"""HCPV thickness map."""
# It is important that input are pure numpies, not masked
logger.info("Name is %s", __name__)
logger.info("Import roff...")
grd = xtgeo.grid_from_file(ROFF1_GRID, fformat="roff")
st = xtgeo.gridproperty_from_file(ROFF1_PROPS, name="Oil_HCPV")
to = xtgeo.gridproperty_from_file(ROFF1_PROPS, name="Oil_bulk")
# get the dz and the coordinates, with no mask (ie get value for outside)
dz = grd.get_dz(asmasked=False)
xc, yc, _zc = grd.get_xyz(asmasked=False)
xcv = ma.filled(xc.values3d)
ycv = ma.filled(yc.values3d)
dzv = ma.filled(dz.values3d)
hcpfz = ma.filled(st.values3d, fill_value=0.0)
tov = ma.filled(to.values3d, fill_value=10)
tov[tov < 1.0e-32] = 1.0e-32
hcpfz = hcpfz * dzv / tov
# make a map... estimate from xc and yc
xmin = xcv.min()
xmax = xcv.max()
ymin = ycv.min()
ymax = ycv.max()
xinc = (xmax - xmin) / 50
yinc = (ymax - ymin) / 50
logger.debug(
"xmin xmax ymin ymax, xinc, yinc: %s %s %s %s %s %s",
xmin,
xmax,
ymin,
ymax,
xinc,
yinc,
)
hcmap = RegularSurface(
ncol=50,
nrow=50,
xinc=xinc,
yinc=yinc,
xori=xmin,
yori=ymin,
values=np.zeros((50, 50)),
)
hcmap2 = RegularSurface(
ncol=50,
nrow=50,
xinc=xinc,
yinc=yinc,
xori=xmin,
yori=ymin,
values=np.zeros((50, 50)),
)
zp = np.ones((grd.ncol, grd.nrow, grd.nlay))
# now make hcpf map
t1 = xtg.timer()
hcmap.hc_thickness_from_3dprops(
xprop=xcv,
yprop=ycv,
dzprop=dzv,
hcpfzprop=hcpfz,
zoneprop=zp,
zone_minmax=(1, 1),
)
assert hcmap.values.mean() == pytest.approx(1.447, abs=0.1)
t2 = xtg.timer(t1)
logger.info("Speed basic is %s", t2)
t1 = xtg.timer()
hcmap2.hc_thickness_from_3dprops(
xprop=xcv,
yprop=ycv,
dzprop=dzv,
hcpfzprop=hcpfz,
zoneprop=zp,
coarsen=2,
zone_avg=True,
zone_minmax=(1, 1),
mask_outside=True,
)
t2 = xtg.timer(t1)
logger.info("Speed zoneavg coarsen 2 is %s", t2)
if generate_plot:
hcmap.quickplot(filename=join(tmpdir, "quickplot_hcpv.png"))
hcmap2.quickplot(filename=join(tmpdir, "quickplot_hcpv_zavg_coarsen.png"))
0
Source : test_grid_property.py
with GNU Lesser General Public License v3.0
from equinor
with GNU Lesser General Public License v3.0
from equinor
def test_io_roff_discrete(tmpdir):
"""Import ROFF discrete property; then export to ROFF int."""
logger.info("Name is {}".format(__name__))
po = xtgeo.gridproperty_from_file(TESTFILE8, fformat="roff", name="Zone")
logger.info("\nCodes ({})\n{}".format(po.ncodes, po.codes))
# tests:
assert po.ncodes == 3
logger.debug(po.codes[3])
assert po.codes[3] == "Below_Low_reek"
# export discrete to ROFF ...TODO
po.to_file(
os.path.join(tmpdir, "reek_zone_export.roff"), name="Zone", fformat="roff"
)
# fix some zero values (will not be fixed properly as grid ACTNUM differs?)
val = po.values
val = npma.filled(val, fill_value=3) # trick
print(val.min(), val.max())
po.values = val
print(po.values.min(), po.values.max())
po.values[:, :, 13] = 1 # just for fun test
po.to_file(
os.path.join(tmpdir, "reek_zonefix_export.roff"), name="ZoneFix", fformat="roff"
)
@given(grid_properties())
0
Source : read_fes2012.py
with GNU General Public License v3.0
from hrmartens
with GNU General Public License v3.0
from hrmartens
def main(filename):
# Read In NetCDF File
f = netCDF4.Dataset(filename)
#print(f.dimensions['lon'])
#print(f.dimensions['lat'])
#print(f.variables)
lon = f.variables['lon'][:]
lat = f.variables['lat'][:]
amp = f.variables['Ha'][:]
pha = f.variables['Hg'][:]
f.close()
# Replace Masked Data with Fill Values
amp = np.ma.filled(amp,fill_value=0.)
pha = np.ma.filled(pha,fill_value=0.)
# Convert Amplitude from Centimeters to Meters
amp = np.divide(amp,100.)
# Write Phase to New Array to Make it Writeable
pha = np.divide(pha,1.)
# NaN 'FillValue' for FES2012 is 1.e+10 (in cm)
# Obtained from: amp_var = f.variables['Ha']; print(amp_var)
# Set NaN Values to Zero
# Note: This Should Already Be Taken Care of With the Unmasking Above
amp[amp > 10.**5.] = 0.
pha[pha > 10.**5.] = 0.
# Save Arrays in Original Format
lon1dseq = lon
lat1dseq = lat
amp2darr = amp
pha2darr = pha
# Reformat Load Points into 1D Vectors
grid_olon, grid_olat = sc.meshgrid(lon,lat)
olon = grid_olon.flatten()
olat = grid_olat.flatten()
amp = amp.flatten()
pha = pha.flatten()
# NaN 'FillValue' for FES2012 is 1.e+10 (in cm)
# Set NaN Values to Zero
nanidx = np.where(amp > 10**5)
amp[nanidx] = 0.
pha[nanidx] = 0.
# Return Parameters
return olat,olon,amp,pha,lat1dseq,lon1dseq,amp2darr,pha2darr
0
Source : read_fes2014.py
with GNU General Public License v3.0
from hrmartens
with GNU General Public License v3.0
from hrmartens
def main(filename):
# Read In NetCDF File
f = netCDF4.Dataset(filename)
#print(f.dimensions['lon'])
#print(f.dimensions['lat'])
#print(f.variables)
lon = f.variables['lon'][:]
lat = f.variables['lat'][:]
amp = f.variables['amplitude'][:]
pha = f.variables['phase'][:]
f.close()
# Replace Masked Data with Fill Values
amp = np.ma.filled(amp,fill_value=0.)
pha = np.ma.filled(pha,fill_value=0.)
# Convert Amplitude from Centimeters to Meters
amp = np.divide(amp,100.)
# Write Phase to New Array to Make it Writeable
pha = np.divide(pha,1.)
# NaN 'FillValue' for FES2014 is 1.84467e+19 (in cm)
# Obtained from: amp_var = f.variables['amplitude']; print(amp_var)
# Set NaN Values to Zero
# Note: This Should Already Be Taken Care of With the Unmasking Above
amp[amp > 10.**17.] = 0.
pha[pha > 10.**17.] = 0.
# Save Arrays in Original Format
lon1dseq = lon
lat1dseq = lat
amp2darr = amp
pha2darr = pha
# Reformat Load Points into 1D Vectors
grid_olon, grid_olat = sc.meshgrid(lon,lat)
olon = grid_olon.flatten()
olat = grid_olat.flatten()
amp = amp.flatten()
pha = pha.flatten()
# NaN 'FillValue' for FES2014 is 1.84467e+19 (in cm)
# Set NaN Values to Zero
nanidx = np.where(amp > 10**10)
amp[nanidx] = 0.
pha[nanidx] = 0.
# Return Parameters
return olat,olon,amp,pha,lat1dseq,lon1dseq,amp2darr,pha2darr
0
Source : read_hamtide11a.py
with GNU General Public License v3.0
from hrmartens
with GNU General Public License v3.0
from hrmartens
def main(filename):
# Read In NetCDF File
f = netCDF4.Dataset(filename)
#print f.variables
lon = f.variables['LON'][:]
lat = f.variables['LAT'][:]
real = f.variables['RE'][:]
imag = f.variables['IM'][:]
amp = f.variables['AMPL'][:]
pha = f.variables['PHAS'][:]
#real_var = f.variables['re']; print(real_var)
#amp = np.absolute(real + np.multiply(imag,1.j))
#pha = np.multiply(np.arctan2(imag,real),180./pi)
f.close()
# Replace Masked Data with Fill Values
amp = np.ma.filled(amp,fill_value=0.)
pha = np.ma.filled(pha,fill_value=0.)
# Convert Amplitude from Centimeters to Meters
amp = np.divide(amp,100.)
# Write Phase to New Array to Make it Writeable
pha = np.divide(pha,1.)
# Save Arrays in Original Format
lon1dseq = lon
lat1dseq = lat
amp2darr = amp
pha2darr = pha
# Reformat Load Points into 1D Vectors
grid_olon, grid_olat = sc.meshgrid(lon,lat)
olon = grid_olon.flatten()
olat = grid_olat.flatten()
amp = amp.flatten()
pha = pha.flatten()
# Return Parameters
return olat,olon,amp,pha,lat1dseq,lon1dseq,amp2darr,pha2darr
0
Source : mrecords.py
with BSD 3-Clause "New" or "Revised" License
from leobago
with BSD 3-Clause "New" or "Revised" License
from leobago
def __setattr__(self, attr, val):
"""
Sets the attribute attr to the value val.
"""
# Should we call __setmask__ first ?
if attr in ['mask', 'fieldmask']:
self.__setmask__(val)
return
# Create a shortcut (so that we don't have to call getattr all the time)
_localdict = object.__getattribute__(self, '__dict__')
# Check whether we're creating a new field
newattr = attr not in _localdict
try:
# Is attr a generic attribute ?
ret = object.__setattr__(self, attr, val)
except Exception:
# Not a generic attribute: exit if it's not a valid field
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
optinfo = ndarray.__getattribute__(self, '_optinfo') or {}
if not (attr in fielddict or attr in optinfo):
raise
else:
# Get the list of names
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
# Check the attribute
if attr not in fielddict:
return ret
if newattr:
# We just added this one or this setattr worked on an
# internal attribute.
try:
object.__delattr__(self, attr)
except Exception:
return ret
# Let's try to set the field
try:
res = fielddict[attr][:2]
except (TypeError, KeyError):
raise AttributeError(f'record array has no attribute {attr}')
if val is masked:
_fill_value = _localdict['_fill_value']
if _fill_value is not None:
dval = _localdict['_fill_value'][attr]
else:
dval = val
mval = True
else:
dval = filled(val)
mval = getmaskarray(val)
obj = ndarray.__getattribute__(self, '_data').setfield(dval, *res)
_localdict['_mask'].__setitem__(attr, mval)
return obj
def __getitem__(self, indx):
0
Source : mrecords.py
with GNU Affero General Public License v3.0
from nccgroup
with GNU Affero General Public License v3.0
from nccgroup
def __setattr__(self, attr, val):
"""
Sets the attribute attr to the value val.
"""
# Should we call __setmask__ first ?
if attr in ['mask', 'fieldmask']:
self.__setmask__(val)
return
# Create a shortcut (so that we don't have to call getattr all the time)
_localdict = object.__getattribute__(self, '__dict__')
# Check whether we're creating a new field
newattr = attr not in _localdict
try:
# Is attr a generic attribute ?
ret = object.__setattr__(self, attr, val)
except:
# Not a generic attribute: exit if it's not a valid field
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
optinfo = ndarray.__getattribute__(self, '_optinfo') or {}
if not (attr in fielddict or attr in optinfo):
exctype, value = sys.exc_info()[:2]
raise exctype(value)
else:
# Get the list of names
fielddict = ndarray.__getattribute__(self, 'dtype').fields or {}
# Check the attribute
if attr not in fielddict:
return ret
if newattr:
# We just added this one or this setattr worked on an
# internal attribute.
try:
object.__delattr__(self, attr)
except:
return ret
# Let's try to set the field
try:
res = fielddict[attr][:2]
except (TypeError, KeyError):
raise AttributeError("record array has no attribute %s" % attr)
if val is masked:
_fill_value = _localdict['_fill_value']
if _fill_value is not None:
dval = _localdict['_fill_value'][attr]
else:
dval = val
mval = True
else:
dval = filled(val)
mval = getmaskarray(val)
obj = ndarray.__getattribute__(self, '_data').setfield(dval, *res)
_localdict['_mask'].__setitem__(attr, mval)
return obj
def __getitem__(self, indx):
0
Source : raster.py
with GNU General Public License v2.0
from OSGeo
with GNU General Public License v2.0
from OSGeo
def _pred_fun(img, estimator):
"""Prediction function for classification or regression response
Parameters
----
img : numpy.ndarray
3d ndarray of raster data with the dimensions in order of
(band, rows, columns).
estimator : estimator object implementing 'fit'
The object to use to fit the data.
Returns
-------
numpy.ndarray
2d numpy array representing a single band raster containing the
classification or regression result.
"""
n_features, rows, cols = img.shape[0], img.shape[1], img.shape[2]
# reshape each image block matrix into a 2D matrix
# first reorder into rows,cols,bands(transpose)
# then resample into 2D array (rows=sample_n, cols=band_values)
n_samples = rows * cols
flat_pixels = img.transpose(1, 2, 0).reshape((n_samples, n_features))
# create mask for NaN values and replace with number
flat_pixels_mask = flat_pixels.mask.copy()
flat_pixels = np.ma.filled(flat_pixels, -99999)
# prediction
result = estimator.predict(flat_pixels)
# replace mask and fill masked values with nodata value
result = np.ma.masked_array(result, mask=flat_pixels_mask.any(axis=1))
# reshape the prediction from a 1D matrix/list
# back into the original format [band, row, col]
result = result.reshape((1, rows, cols))
return result
@staticmethod
0
Source : raster.py
with GNU General Public License v2.0
from OSGeo
with GNU General Public License v2.0
from OSGeo
def _prob_fun(img, estimator):
"""Class probabilities function
Parameters
----------
img : numpy.ndarray
3d numpy array of raster data with the dimensions in order of
(band, row, column).
estimator : estimator object implementing 'fit'
The object to use to fit the data.
Returns
-------
numpy.ndarray
Multi band raster as a 3d numpy array containing the probabilities
associated with each class. ndarray dimensions are in the order of
(class, row, column).
"""
n_features, rows, cols = img.shape[0], img.shape[1], img.shape[2]
# reshape each image block matrix into a 2D matrix
# first reorder into rows,cols,bands(transpose)
# then resample into 2D array (rows=sample_n, cols=band_values)
n_samples = rows * cols
flat_pixels = img.transpose(1, 2, 0).reshape((n_samples, n_features))
# create mask for NaN values and replace with number
flat_pixels_mask = flat_pixels.mask.copy()
flat_pixels = np.ma.filled(flat_pixels, -99999)
# predict probabilities
result = estimator.predict_proba(flat_pixels)
# reshape class probabilities back to 3D [iclass, rows, cols]
result = result.reshape((rows, cols, result.shape[1]))
flat_pixels_mask = flat_pixels_mask.reshape((rows, cols, n_features))
# flatten mask into 2d
mask2d = flat_pixels_mask.any(axis=2)
mask2d = np.where(mask2d != mask2d.min(), True, False)
mask2d = np.repeat(mask2d[:, :, np.newaxis], result.shape[2], axis=2)
# convert proba to masked array using mask2d
result = np.ma.masked_array(result, mask=mask2d, fill_value=np.nan)
# reshape band into raster format [band, row, col]
result = result.transpose(2, 0, 1)
return result
@staticmethod
See More Examples