[SciPy-User] GIS Raster Calculation: Combination of NumPy Arrays and SciPy Density function

josef.pktd at gmail.com josef.pktd at gmail.com
Fri Jan 7 08:18:43 EST 2011 

On Fri, Jan 7, 2011 at 7:08 AM,  <josef.pktd at gmail.com> wrote:
> On Fri, Jan 7, 2011 at 6:50 AM, Johannes Radinger <JRadinger at gmx.at> wrote:
>>
>> -------- Original-Nachricht --------
>>> Datum: Thu, 6 Jan 2011 11:50:20 -0500
>>> Von: josef.pktd at gmail.com
>>> An: SciPy Users List <scipy-user at scipy.org>
>>> Betreff: Re: [SciPy-User] GIS Raster Calculation: Combination of NumPy Arrays and SciPy Density function
>>
>>> On Thu, Jan 6, 2011 at 9:59 AM, Johannes Radinger <JRadinger at gmx.at>
>>> wrote:
>>> > Hello...
>>> >
>>> > I am working mostly in ArcGIS but as it's Raster Calculator resp. Map
>>> Algebra is limited I want to use NumPy Arrays and perform calcultations in
>>> SciPy.
>>> >
>>> > I can get my Rasterfiles in ArcGIS quite easily transformed/exported
>>> into a NumPy Array with the function:
>>> >
>>> > newArray = arcpy.RasterToNumPyArray(inRaster)
>>> >
>>> > This Raster/Array contains distance values (ranging from around -20000
>>> to + 20000 float).
>>> >
>>> > Some time before you helped with some distance-density functions which I
>>> want to apply now. New Values should be assigned to the cells/elements in
>>> the array by following function:
>>> >
>>> > DensityRaster = p * stats.norm.pdf(x, loc=m, scale=s)
>>> >
>>> > where p = 0.3, m=0 and s=200, the x-value is the distance value from the
>>> single cells in newArray. Can I just define x=newArray ? and use an array
>>> as input variable?
>>>
>>> Yes, just try it out. It's fully vectorized.
>>> loc=0 is the default, so in this case you wouldn't need to specify loc=m.
>>>
>>> If your array is very large, there is a shortcut to (maybe) avoid some
>>> intermediate arrays.
>>
>> Oh thank you, that works perfectly. What do you mean with a shortcut?

use the internal method _pdf . It avoids the generic overhead, but the
user is responsible that it makes sense

>>> scale = 2
>>> stats.norm.pdf(np.ones((3,3)), scale=scale)
array([[ 0.17603266,  0.17603266,  0.17603266],
       [ 0.17603266,  0.17603266,  0.17603266],
       [ 0.17603266,  0.17603266,  0.17603266]])

>>> stats.norm._pdf(np.ones((3,3))/scale)/scale
array([[ 0.17603266,  0.17603266,  0.17603266],
       [ 0.17603266,  0.17603266,  0.17603266],
       [ 0.17603266,  0.17603266,  0.17603266]])

It might save a bit of memory and time, but I never actually compared
the performance

Josef

>>
>> And just another short question about the stats.norm.cdf function: That is calculating the cumulative probabilty from the smallest value (most negative x) to the point x, but is there also a function which gives the cumulative sum away from 0 to point x in both directions (negativ and positive? Like an integration of the probabilty function from 0 to (-)x?
>
> cdf is the standard definition of the cdf, integration from lower
> support limit, attribute .a, for normal it's -np.inf, to point x.
>
> there is also sf(x) = 1-cdf(x)
>
> in general
> Prob(0<y<=x) =  cdf(x) - cdf(0)
>
> for symmetric centered  Prob(0<y<=x) = (1-2*cdf(-x))/2 , but I suspect
> for the normal you don't gain anything compared to cdf(x) - 0.5
>
> Josef
>>
>> thank you
>> /j
>>>
>>> Josef
>>>
>>>
>>> >
>>> > thank you
>>> > /j
>>> >
>>> >
>>> > -------- Original-Nachricht --------
>>> >> Datum: Thu, 6 Jan 2011 07:10:58 -0500
>>> >> Von: josef.pktd at gmail.com
>>> >> An: SciPy Users List <scipy-user at scipy.org>
>>> >> Betreff: Re: [SciPy-User] solving integration, density function
>>> >
>>> >> On Thu, Jan 6, 2011 at 7:01 AM, Johannes Radinger <JRadinger at gmx.at>
>>> >> wrote:
>>> >> > Thank you for the simplification of the formula,
>>> >> > but I still get a different result in the case
>>> >> > when x<s1 (eg. x=2, s1=3).
>>> >> >
>>> >> > here a code to try:
>>> >> >
>>> >> > ********************
>>> >> > import math
>>> >> > from scipy import stats
>>> >> >
>>> >> > s1 = 3
>>> >> > m = 0
>>> >> > p = 1
>>> >> > x = 2
>>> >> >
>>> >> > func = stats.norm.pdf(x, loc=m, scale=(s1))
>>> >> > func2 = (1/(s1*math.sqrt(2*math.pi)) * math.exp(-0.5*((x-m)/s1)**2))
>>> >> >
>>> >> > print func
>>> >> > print func2
>>> >> > ********************************
>>> >>
>>> >> use floats, I think you just run into integer division
>>> >>
>>> >> (x-m)/s1
>>> >>
>>> >> Josef
>>> >>
>>> >> >
>>> >> > /j
>>> >> >
>>> >> > -------- Original-Nachricht --------
>>> >> >> Datum: Thu, 6 Jan 2011 05:48:25 -0600
>>> >> >> Von: Warren Weckesser <warren.weckesser at enthought.com>
>>> >> >> An: SciPy Users List <scipy-user at scipy.org>
>>> >> >> Betreff: Re: [SciPy-User] solving integration, density function
>>> >> >
>>> >> >> On Thu, Jan 6, 2011 at 5:27 AM, Johannes Radinger <JRadinger at gmx.at>
>>> >> >> wrote:
>>> >> >>
>>> >> >> > Hey
>>> >> >> >
>>> >> >> > Last time you helped me a lot with my normal
>>> >> >> > probabilty density function. My problem now is
>>> >> >> > quite simple, I think it's just a problem with
>>> >> >> > the syntax (brackets):
>>> >> >> >
>>> >> >> > There are two ways to calculate the pdf, with the
>>> >> >> > stats-function and with pure mathematically, but
>>> >> >> > the give different results and I can't find the
>>> >> >> > where I make the mistake:
>>> >> >> >
>>> >> >> >
>>> >> >> > func1 = stats.norm.pdf(x, loc=m, scale=(s1))
>>> >> >> > func2 =
>>> >> >> >
>>> 1/((s1)*(math.sqrt(2*math.pi))))*(math.exp(((-0.5)*((x-m)/(s1)))**2)
>>> >> >> >
>>> >> >> > Where is the problem
>>> >> >> >
>>> >> >>
>>> >> >>
>>> >> >> func2 = 1/(s1*math.sqrt(2*math.pi)) * math.exp(-0.5*((x-m)/s1)**2)
>>> >> >>
>>> >> >>
>>> >> >> Warren
>>> >> >>
>>> >> >>
>>> >> >>
>>> >> >> > thank you...
>>> >> >> >
>>> >> >> > /j
>>> >> >> >
>>> >> >> > -------- Original-Nachricht --------
>>> >> >> > > Datum: Tue, 21 Dec 2010 09:18:15 -0500
>>> >> >> > > Von: Skipper Seabold <jsseabold at gmail.com>
>>> >> >> > > An: SciPy Users List <scipy-user at scipy.org>
>>> >> >> > > Betreff: Re: [SciPy-User] solving integration, density function
>>> >> >> >
>>> >> >> > > On Tue, Dec 21, 2010 at 7:48 AM, Johannes Radinger
>>> >> <JRadinger at gmx.at>
>>> >> >> > > wrote:
>>> >> >> > > >
>>> >> >> > > > -------- Original-Nachricht --------
>>> >> >> > > >> Datum: Tue, 21 Dec 2010 13:20:47 +0100
>>> >> >> > > >> Von: Gregor Thalhammer <Gregor.Thalhammer at gmail.com>
>>> >> >> > > >> An: SciPy Users List <scipy-user at scipy.org>
>>> >> >> > > >> Betreff: Re: [SciPy-User] solving integration, density
>>> function
>>> >> >> > > >
>>> >> >> > > >>
>>> >> >> > > >> Am 21.12.2010 um 12:06 schrieb Johannes Radinger:
>>> >> >> > > >>
>>> >> >> > > >> > Hello,
>>> >> >> > > >> >
>>> >> >> > > >> > I am really new to python and Scipy.
>>> >> >> > > >> > I want to solve a integrated function with a python script
>>> >> >> > > >> > and I think Scipy should do that :)
>>> >> >> > > >> >
>>> >> >> > > >> > My task:
>>> >> >> > > >> >
>>> >> >> > > >> > I do have some variables (s, m, K,) which are now
>>> absolutely
>>> >> set,
>>> >> >> > but
>>> >> >> > > in
>>> >> >> > > >> future I'll get the values via another process of pyhton.
>>> >> >> > > >> >
>>> >> >> > > >> > s = 400
>>> >> >> > > >> > m = 0
>>> >> >> > > >> > K = 1
>>> >> >> > > >> >
>>> >> >> > > >> > And have have following function:
>>> >> >> > > >> > (1/((s*K)*sqrt(2*pi)))*exp((-1/2*(((x-m)/s*K))^2) which is
>>> the
>>> >> >> > > density
>>> >> >> > > >> function of the normal distribution a symetrical curve with
>>> the
>>> >> >> mean
>>> >> >> > > (m) of
>>> >> >> > > >> 0.
>>> >> >> > > >> >
>>> >> >> > > >> > The total area under the curve is 1 (100%) which is for an
>>> >> >> > > integration
>>> >> >> > > >> from -inf to +inf.
>>> >> >> > > >> > I want to know x in the case of 99%: meaning that the
>>> integral
>>> >> >> (-x
>>> >> >> > to
>>> >> >> > > >> +x) of the function is 0.99. Due to the symetry of the curve
>>> you
>>> >> >> can
>>> >> >> > > also set
>>> >> >> > > >> the integral from 0 to +x equal to (0.99/2):
>>> >> >> > > >> >
>>> >> >> > > >> > 0.99 =
>>> >> >> integral((1/((s*K)*sqrt(2*pi)))*exp((-1/2*(((x-m)/s*K))^2)),
>>> >> >> > > -x,
>>> >> >> > > >> x)
>>> >> >> > > >> > resp.
>>> >> >> > > >> > (0.99/2) =
>>> >> >> > > integral((1/((s*K)*sqrt(2*pi)))*exp((-1/2*(((x-m)/s*K))^2)),
>>> >> >> > > >> 0, x)
>>> >> >> > > >> >
>>> >> >> > > >> > How can I solve that question in Scipy/python
>>> >> >> > > >> > so that I get x in the end. I don't know how to write
>>> >> >> > > >> > the code...
>>> >> >> > > >>
>>> >> >> > > >>
>>> >> >> > > >> --->
>>> >> >> > > >> erf(x[, out])
>>> >> >> > > >>
>>> >> >> > > >>     y=erf(z) returns the error function of complex argument
>>> >> defined
>>> >> >> > > as
>>> >> >> > > >>     as 2/sqrt(pi)*integral(exp(-t**2),t=0..z)
>>> >> >> > > >> ---
>>> >> >> > > >>
>>> >> >> > > >> from scipy.special import erf, erfinv
>>> >> >> > > >> erfinv(0.99)*sqrt(2)
>>> >> >> > > >>
>>> >> >> > > >>
>>> >> >> > > >> Gregor
>>> >> >> > > >>
>>> >> >> > > >
>>> >> >> > > >
>>> >> >> > > > Thank you Gregor,
>>> >> >> > > > I only understand a part of your answer... I know that the
>>> >> integral
>>> >> >> of
>>> >> >> > > the density function is a error function and I know that the
>>> >> argument
>>> >> >> > "from
>>> >> >> > > scipy.special import erf, erfinv" is to load the module.
>>> >> >> > > >
>>> >> >> > > > But how do I write the code including my orignial function so
>>> >> that I
>>> >> >> > can
>>> >> >> > > modify it (I have also another function I want to integrate).
>>> how
>>> >> do i
>>> >> >> > > start? I want to save the whole code to a python-script I can
>>> then
>>> >> >> load
>>> >> >> > e.g.
>>> >> >> > > into ArcGIS where I want to use the value of x for further
>>> >> >> calculations.
>>> >> >> > > >
>>> >> >> > >
>>> >> >> > > Are you always integrating densities?  If so, you don't want to
>>> >> use
>>> >> >> > > integrals probably, but you could use scipy.stats
>>> >> >> > >
>>> >> >> > > erfinv(.99)*np.sqrt(2)
>>> >> >> > > 2.5758293035489004
>>> >> >> > >
>>> >> >> > > from scipy import stats
>>> >> >> > >
>>> >> >> > > stats.norm.ppf(.995)
>>> >> >> > > 2.5758293035489004
>>> >> >> > >
>>> >> >> > > Skipper
>>> >> >> > > _______________________________________________
>>> >> >> > > SciPy-User mailing list
>>> >> >> > > SciPy-User at scipy.org
>>> >> >> > > http://mail.scipy.org/mailman/listinfo/scipy-user
>>> >> >> >
>>> >> >> > --
>>> >> >> > NEU: FreePhone - kostenlos mobil telefonieren und surfen!
>>> >> >> > Jetzt informieren: http://www.gmx.net/de/go/freephone
>>> >> >> > _______________________________________________
>>> >> >> > SciPy-User mailing list
>>> >> >> > SciPy-User at scipy.org
>>> >> >> > http://mail.scipy.org/mailman/listinfo/scipy-user
>>> >> >> >
>>> >> >
>>> >> > --
>>> >> > NEU: FreePhone - kostenlos mobil telefonieren und surfen!
>>> >> > Jetzt informieren: http://www.gmx.net/de/go/freephone
>>> >> > _______________________________________________
>>> >> > SciPy-User mailing list
>>> >> > SciPy-User at scipy.org
>>> >> > http://mail.scipy.org/mailman/listinfo/scipy-user
>>> >> >
>>> >> _______________________________________________
>>> >> SciPy-User mailing list
>>> >> SciPy-User at scipy.org
>>> >> http://mail.scipy.org/mailman/listinfo/scipy-user
>>> >
>>> > --
>>> > Sicherer, schneller und einfacher. Die aktuellen Internet-Browser -
>>> > jetzt kostenlos herunterladen! http://portal.gmx.net/de/go/atbrowser
>>> > _______________________________________________
>>> > SciPy-User mailing list
>>> > SciPy-User at scipy.org
>>> > http://mail.scipy.org/mailman/listinfo/scipy-user
>>> >
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>> --
>> NEU: FreePhone - kostenlos mobil telefonieren und surfen!
>> Jetzt informieren: http://www.gmx.net/de/go/freephone
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>

More information about the SciPy-User mailing list