Mailman 3 April 2014 - SciPy-User

How to map spherical coordinate data to cartesian grid?
by Steve Miller April 24, 2014

April 24, 2014

Dear all. I was send 4 long 1D arrays (flattened) representing 3D data in spherical coordinates: rho: density array r: spherical coordinate r array th: spherical coordinate theta array phi: spherical coordinate phi array where the r, th, and phi arrays are the coordinates for the density array rho. These arrays are not uniformly spaces in r, th and phi. I would like to map rho to cartesian coordinates of a 3D cube. IE, define: X, Y, Z = mgrid[etc...] Is there a handy scipy routine that will map/interpolate these rho values in spherical coordinates across my new cartesian box? Naively I thought I could use something like map_coordinates but I am unclear how I would use this for this calculation. Could someone provide a simple example script for how I could use scipy to do this? Thank you.

2 2

numpy array assignment
by marco cammarata April 18, 2014

April 18, 2014

Hi, I have used many times syntax like import numpy as np def result_of_long_calculation_returning_the_right_shape(): return np.random.random( (2,3) ) N=10 out = np.empty( (N,2,3) ) for i in range(N): out[i] = result_of_long_calculation_returning_the_right_shape() if I have to loop trough the axis 1 I would do something like N=10 out = np.empty( (2,N,3) ) for i in range(N): out[:,i,:] = result_of_long_calculation_returning_the_right_shape() The problem is that - in the function I am writing - the number of dimensions and the axis to loop trough are not known in advance. Is there a smart what to do it ? I would expect that a function like np.assign(array,element=0,axis=0) that but I can't find such function. For the moment I ended up doing something really really inelegant: res = result_of_long_calculation_returning_the_right_shape() s = "out[" + (":,")*axis + "i" + (",:")*(nDim-axis-1) + "]=res" exec(s) Thanks a lot in advance, marco

3 3

scipy.signal vs Matlab: filtfilt and reflection
by John Krasting - NOAA Federal April 15, 2014

April 15, 2014

Hi Scipy Users - Am I correct in reading that filtfilt in scipy.signal (v. 13.0) does not extrapolate data at the beginning and the end of a time series when using the filtfilt function? Based on the Matlab code ( http://chmielowski.eu/POLITECHNIKA/Dydaktyka/AUTOMATYKA/AutoLab/Matlab/TOOL…), I see that there is a block of code that does the reflection: % Extrapolate beginning and end of data sequence using a "reflection % method". Slopes of original and extrapolated sequences match at % the end points. % This reduces end effects. y = [2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)]; I did not see this in the source for scipy's filtfilt function. Are there any plans to add this as a possible option to the function? Thanks, John K.

2 5

Re: [SciPy-User] Can create, but not access, netCDF variables with no dimension
by Mannucci, Anthony J (335G) April 14, 2014

April 14, 2014

Thanks for your help using getValue. That should help with the example I provided. However, the variable does not show up in the keys(). Here is an ncdump excerpt. More information: I'm using scipy 0.10. Perhaps this issue has been resolved in a later version? -Tony netcdf apr2011storm_s008D110429_0100_110501_0000 { dimensions: time = UNLIMITED ; // (48 currently) lon = 72 ; lat = 36 ; lev = 29 ; ilev = 29 ; mlon = 81 ; mlat = 97 ; mlev = 32 ; imlev = 32 ; mtimedim = 3 ; latlon = 2 ; dtidedim = 2 ; sdtidedim = 10 ; datelen = 24 ; filelen = 1024 ; variables: double time(time) ; time:long_name = "time" ; time:units = "minutes since 2011-3-21 0:0:0" ; time:initial_year = 2011 ; time:initial_day = 80 ; time:initial_mtime = 80, 0, 0 ; double lon(lon) ; lon:long_name = "geographic longitude (-west, +east)" ; lon:units = "degrees_east" ; double lat(lat) ; lat:long_name = "geographic latitude (-south, +north)" ; lat:units = "degrees_north" ; double lev(lev) ; lev:long_name = "midpoint levels" ; lev:short_name = "ln(p0/p)" ; lev:units = "" ; lev:positive = "up" ; lev:standard_name = "atmosphere_ln_pressure_coordinate" ; lev:formula_terms = "p0: p0 lev: lev" ; lev:formula = "p(k) = p0 * exp(-lev(k))" ; double ilev(ilev) ; ilev:long_name = "interface levels" ; ilev:short_name = "ln(p0/p)" ; ilev:units = "" ; ilev:positive = "up" ; ilev:standard_name = "atmosphere_ln_pressure_coordinate" ; ilev:formula_terms = "p0: p0 lev: ilev" ; ilev:formula = "p(k) = p0 * exp(-ilev(k))" ; <<< More stuff. The following don't show up>>> double p0 ; p0:long_name = "Reference pressure" ; p0:units = "millibars" ; double p0_model ; p0_model:long_name = "Reference pressure (as used by the model)" ; p0_model:units = "microbars" ; -- Tony Mannucci Supervisor, Ionospheric and Atmospheric Remote Sensing Group Mail-Stop 138-308, Tel > (818) 354-1699 Jet Propulsion Laboratory, Fax > (818) 393-5115 California Institute of Technology, Email > Tony.Mannucci(a)jpl.nasa.gov 4800 Oak Grove Drive, http://scienceandtechnology.jpl.nasa.gov/people/a_mannucci/ Pasadena, CA 91109

2 1

Saving into file.txt
by Antonelli Maria Rosaria April 14, 2014

April 14, 2014

Hi, I need to save same data arrays into a file.txt. I would like to save them as different columns. At the moment I can manage to save them as different lines. Can you please give me some hints ? Best regards, Rosa

3 3

create sparse matrix from diagonal arrays
by Darcoux Christine April 14, 2014

April 14, 2014

I would like to solve a PDE using an implicit finite volume discretization. Since the problem is 3D, the matrix should be septadiagonal. I have computed the 7 diagonals (for top, bottom, east, west, north south and central point) in 7 different arrays. Now, I would like to create a sparse matrix object from these arrays, but I feel a little lost with the documentation of scipy.sparse.spdiags. Help would be greatly appreciated. Christine

3 3

Re: [SciPy-User] defining new distributions with, rv_continuous
by Francesco De Gasperin April 14, 2014

April 14, 2014

Many thanks for the quick help. I see that the overhead is quite large but I don't look for high performances now, I need to generate only 1000 numbers. For interested people, I attach here a more flexible code which should work as it numerically normalizes the pdf given "a" and "b". ##################################################### from scipy.stats import distributions from scipy import integrate class test(distributions.rv_continuous): def __init__(self,a=-np.inf,b=np.inf,name=''): super().__init__(a=a,b=b,name=name) self.integral = integrate.quad(self.f, self.a, self.b)[0] def f(self, x): # in this case is a gaussian, but it could be a more complex pdf sigma = 1. c = 0. return np.exp( -(x-c)**2 / (2*sigma**2) ) def _pdf(self, x): return self.f(x)/self.integral t_obj = test(a=0,name="test") ###################################################### thanks, Francesco

1 0

defining new distributions with rv_continuous
by Francesco De Gasperin April 14, 2014

April 14, 2014

Hi, I have some problems in defining a new distribution using rv_continuous. As an example I tried to create a gaussian defined only for positive values of the random variable. This is the code: ########################################### from scipy.stats import distributions class gauss_pos(distributions.rv_continuous): def _pdf(self, x): if x <= 0: return 0 sigma = 1. c = 0. return 1/(np.sqrt(2*np.pi)*sigma) * np.exp( -(x-c)**2 / (2*sigma**2) ) t_obj = gauss_pos(name="gauss_pos") print (t_obj.rvs(size=10)) ############################################# but when I run it I always end up with: OverflowError: (34, 'Numerical result out of range') Then, I removed the "if x <= 0: return 0" line and I tried to define the "a" and "b" parameters when creating the t_obj. Defining only the "a=0" produces the same error, while using something as: t_obj = gauss_pos(a=0,b=2) instead produces a result of this type: [ 2. 2. 1.5460072 1.09621782 1.62184086 0.42107306 2. 1.09281742 0.82957403 2. ] which I still do not understand... it seems that "b" is returned every time the constraint that the random var must be between a and b is violated... any help is appreciated! cheers, Francesco

3 3

Re: [SciPy-User] defining new distributions with, rv_continuous
by Francesco De Gasperin April 14, 2014

April 14, 2014

Thanks for the tip. I am actually working with a quite more complex distribution which is not a simple truncated normal distribution. My example was to simplify the code which generates the problems. > In this specific case, have a look at truncnorm: > > http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.stats.trun… > > import scipy.stats as stats > import numpy as np > stats.truncnorm.pdf(0.1, 0, np.inf) > 0.79390509495402362 > stats.truncnorm.pdf(-0.1, 0, np.inf) > 0.0

1 0

Can create, but not access, netCDF variables with no dimension
by Mannucci, Anthony J (335G) April 14, 2014

April 14, 2014

I am trying to read a netcdf file that has variables created without dimension (they are scalars). Since the file is obtained externally, and written without python, I cannot modify the file. So, the variable I am trying to access is not accessible using scipy.io.netcdf. (I can see the variable with ncdump). I cannot fully reproduce the problem here, although I can reproduce something similar. Here is some code for discussion: --- import scipy.io.netcdf as CDF f = CDF.netcdf_file('simple.nc', 'w') f.createDimension('time', 10) time = f.createVariable('time', 'i', ('time',)) time[:] = np.arange(10) nodimvar = f.createVariable('nodimvar','f8',()) # Guessing this creates a scalar variable. nodimvar = 8.0 # nodimvar[:] fails with an error. Scalar not indexable. f.close() f = CDF.netcdf_file('simple.nc', 'r') print f.variables.keys() print f.variables['time'][:] print f.variables['nodimvar'] # Does not print variable value, just the object. print f.variables['nodimvar'][:] # Error, scalar is not indexable. f.close() ----- The variable name shows up in f.variables.keys(). However, I do not know how to access the variable. Appending [:] to gives an error "IndexError: 0-d arrays can't be indexed.". In the file I was provided, not included here because it's very large, the variable with no dimensions does not even show with up in keys(). I note that the package netCDF4 seems to handle this case correctly. Even zero-dimension variables are accessed via [:]. No error is generated and the value of the scalar is returned. Thanks! -Tony -- Tony Mannucci Supervisor, Ionospheric and Atmospheric Remote Sensing Group Mail-Stop 138-308, Tel > (818) 354-1699 Jet Propulsion Laboratory, Fax > (818) 393-5115 California Institute of Technology, Email > Tony.Mannucci(a)jpl.nasa.gov 4800 Oak Grove Drive, http://scienceandtechnology.jpl.nasa.gov/people/a_mannucci/ Pasadena, CA 91109

2 2