numpy.ma.filled

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 7

3 Source : recfunctions.py
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

    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

    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

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

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

    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

    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

    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

    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

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

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

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

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

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

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

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

    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

    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

    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

    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

    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

    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

    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

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

    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

    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

    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

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

    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

    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

    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

    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

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

    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

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

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

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

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

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

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

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

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

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

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

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

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

    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

    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

    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

    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