Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy

2014-10-14 Thread Ian Henriksen
On Mon, Oct 13, 2014 at 8:39 PM, Charles R Harris charlesr.har...@gmail.com
 wrote:



 On Mon, Oct 13, 2014 at 12:54 PM, Sebastian Berg 
 sebast...@sipsolutions.net wrote:

 On Mo, 2014-10-13 at 13:35 +0200, Daniele Nicolodi wrote:
  Hello,
 
  I have a C++ application that collects float, int or complex data in a
  possibly quite large std::vector. The application has some SWIG
  generated python wrappers that expose this vector to python. However,
  the standard way in which SWIG exposes the data is to create a touple
  and pass this to python, where it is very often converted to a numpy
  array for processing. Of course this is not efficient.
 
  I would like therefore to generate a smarter python interface. In python
  3 I would without doubt go for implementing the buffer protocol, which
  enables seamless interfacing with numpy. However, I need to support also
  python 2 and there the buffer protocol is not as nice.

 Isn't the new buffer protocol in python 2.6 or 2.7 already? There is at
 least a memoryview object in python 2, which maybe could be used to the
 same effect?


 No memoryview in python2.6, but the older buffer protocol it there. Is
 Boost Python an option?

 snip

 Chuck

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


Here's an idea on how to avoid the buffer protocol entirely.
It won't be perfectly optimal, but it is a simple solution that
may be good enough.

It could be that the shape and type inference and the Python
loops that are run when converting the large tuple to a NumPy
array are the main sticking points. In that case, the primary
bottleneck could be eliminated by looping and copying at the C level.
You could use Cython and just copy the needed information
manually into the array. Here's a rough outline of what the
Cython file might look like:

# distutils: language = c++
# distutils: libraries = (whatever shared library you need)
# distutils: library_dirs = (folders where the libraries are)
# distutils: include_dirs = (wherever the headers are)
# cython: boundscheck = False
# cython: wraparound = False

import numpy as np
from libcpp.vector cimport vector

# Here do the call to your C++ library
# This can be done by declaring the
# functions or whatever else you need in
# a cdef extern from header file block
# and then calling the functions to get
# the vector you need.

cdef int[::1] arr_from_vector(vector[int] v):
cdef:
int i
int[::1] a = np.empty(v.size(), int)
for i in xrange(v.size()):
a[i] = v[i]
return np.asarray(a)

def get_data():
# Run your C++ computation here.
return arr_from_vector(v)
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy

2014-10-14 Thread Stefan Seefeld
On Oct 14, 2014 4:40 AM, Charles R Harris charlesr.har...@gmail.com
wrote:



 On Mon, Oct 13, 2014 at 12:54 PM, Sebastian Berg 
sebast...@sipsolutions.net wrote:

 On Mo, 2014-10-13 at 13:35 +0200, Daniele Nicolodi wrote:
  Hello,
 
  I have a C++ application that collects float, int or complex data in a
  possibly quite large std::vector. The application has some SWIG
  generated python wrappers that expose this vector to python. However,
  the standard way in which SWIG exposes the data is to create a touple
  and pass this to python, where it is very often converted to a numpy
  array for processing. Of course this is not efficient.
 
  I would like therefore to generate a smarter python interface. In
python
  3 I would without doubt go for implementing the buffer protocol, which
  enables seamless interfacing with numpy. However, I need to support
also
  python 2 and there the buffer protocol is not as nice.

 Isn't the new buffer protocol in python 2.6 or 2.7 already? There is at
 least a memoryview object in python 2, which maybe could be used to the
 same effect?


 No memoryview in python2.6, but the older buffer protocol it there. Is
Boost Python an option?


Boost.Numpy may be a useful tool to use here.

Stefan

 snip

 Chuck

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy

2014-10-14 Thread Daniele Nicolodi
On 14/10/14 04:39, Charles R Harris wrote:
 On Mon, Oct 13, 2014 at 12:54 PM, Sebastian Berg
 sebast...@sipsolutions.net mailto:sebast...@sipsolutions.net wrote:
 
 On Mo, 2014-10-13 at 13:35 +0200, Daniele Nicolodi wrote:
  Hello,
 
  I have a C++ application that collects float, int or complex data in a
  possibly quite large std::vector. The application has some SWIG
  generated python wrappers that expose this vector to python. However,
  the standard way in which SWIG exposes the data is to create a touple
  and pass this to python, where it is very often converted to a numpy
  array for processing. Of course this is not efficient.
 
  I would like therefore to generate a smarter python interface. In python
  3 I would without doubt go for implementing the buffer protocol, which
  enables seamless interfacing with numpy. However, I need to support also
  python 2 and there the buffer protocol is not as nice.
 
 Isn't the new buffer protocol in python 2.6 or 2.7 already? There is at
 least a memoryview object in python 2, which maybe could be used to the
 same effect?
 
 No memoryview in python2.6, but the older buffer protocol it there. Is
 Boost Python an option?

The old buffer protocol is an option, but it is much less nice than the
new one, as it requires to use numpy.frombuffer() with an exlicit dtype
instead of the siumpler numpy.asarray() sufficient in python 3.

Boost Python may be an option as the codebase already depends on Boost,
but probably not yet on Boost Python. Can you point me to the relevant
documentation, and maybe to an example? One of the problems I have is
that the current wrapping is done auto-magically with SWIG and I would
like to deviate the less possible from that patter.

Thank you!

Cheers,
Daniele

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] numpy.mean slicing in a netCDF file

2014-10-14 Thread Fadzil Mnor
Hi all,
I wrote a script and plot monthly mean zonal wind (from a netcdf file names
uwnd.mon.mean.nc) and I'm not sure if I'm doing it correctly. What I have:


*#this part calculates mean values for january only, from 1980-2010; thus
the index looks like this 396:757:12*

def ujan():
f = nc.Dataset('~/data/ncep/uwnd.mon.mean.nc')
u10_1 = f.variables['uwnd']
u10_2 = np.mean(u10_1[396:757:12,:,38,39:43],axis=0)
return u10_2

uJan = ujan()* #calling function*

*#this part is only to define lon, lat and level *
q = nc.Dataset('~/data/ncep/uwnd.mon.mean.nc')
lon=q.variables['lon']
lat=q.variables['lat']
lev=q.variables['level']

*#for some reason I need to define this unless it gave error length of x
must be number of column in z*

lon=lon[39:43]

*#begin plotting*

clevs=np.arange(-10.,10.,0.5)
fig = plt.figure(figsize=(11, 8))
fig.clf()
ax = fig.add_subplot(111)
ax.axis([97.5, 105., 1000., 10.])
ax.tick_params(direction='out', which='both')
ax.set_xlabel('Lon (degrees)')
ax.set_ylabel('Pressure (mb)')
ax.set_xticks(np.arange(97.5, 105., .5))
ax.set_yticks([1000, 700, 500, 300, 100, 10])
cs=ax.contourf(lon, lev, uJan, clevs, extend='both',cmap='seismic')
plt.title('Zonal winds average (Jan, 1981-2010)')
cax = fig.add_axes([0.99, 0.1, 0.03, 0.8])
aa=fig.colorbar(cs,cax=cax,orientation='vertical')
aa.set_label('m/s')
plt.savefig('~/uwind-crossection-test.png', bbox_inches='tight')
***

the result is attached.
I have no idea how to confirm the result (at least until this email is
written) , but I believe the lower altitude values should be mostly
negative, because over this region, the zonal wind are usually easterly
(thus,negative values), but I got positive values.

Put the information above aside, *I just want to know if my slicing in the
ujan() function is correct*. If it is, then, there must be nothing
wrong(except my above mentioned assumption).
The data file dimension is:
*[time,level,latitude,longitude]*

This part:
*u10_2 = np.mean(u10_1[396:757:12,:,38,39:43],axis=0)*
The line above will calculate the mean of zonal wind (uwnd) in a range of
time index 396 to 757 for each year (january only), for all vertical level,
at latitude index 38 (5 N) and in between longitude index 39 to 43
(97.5E-105E).
I assume it will calculate a 30-year average of zonal wind for january only.
Is this correct?

Thank you.

Fadzil
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy

2014-10-14 Thread Nathaniel Smith
If the goal is to have something that works kind of like the new buffer
protocol but with a wider variety of python versions, then you might find
the old array interface useful:
http://docs.scipy.org/doc/numpy/reference/arrays.interface.html

I always get confused by the history here but I believe that that's the
numpy-only interface that later got cleaned up and generalized to become
the new buffer interface. Numpy itself still supports it.

-n
On 14 Oct 2014 11:19, Daniele Nicolodi dani...@grinta.net wrote:

 On 14/10/14 04:39, Charles R Harris wrote:
  On Mon, Oct 13, 2014 at 12:54 PM, Sebastian Berg
  sebast...@sipsolutions.net mailto:sebast...@sipsolutions.net wrote:
 
  On Mo, 2014-10-13 at 13:35 +0200, Daniele Nicolodi wrote:
   Hello,
  
   I have a C++ application that collects float, int or complex data
 in a
   possibly quite large std::vector. The application has some SWIG
   generated python wrappers that expose this vector to python.
 However,
   the standard way in which SWIG exposes the data is to create a
 touple
   and pass this to python, where it is very often converted to a
 numpy
   array for processing. Of course this is not efficient.
  
   I would like therefore to generate a smarter python interface. In
 python
   3 I would without doubt go for implementing the buffer protocol,
 which
   enables seamless interfacing with numpy. However, I need to
 support also
   python 2 and there the buffer protocol is not as nice.
 
  Isn't the new buffer protocol in python 2.6 or 2.7 already? There is
 at
  least a memoryview object in python 2, which maybe could be used to
 the
  same effect?
 
  No memoryview in python2.6, but the older buffer protocol it there. Is
  Boost Python an option?

 The old buffer protocol is an option, but it is much less nice than the
 new one, as it requires to use numpy.frombuffer() with an exlicit dtype
 instead of the siumpler numpy.asarray() sufficient in python 3.

 Boost Python may be an option as the codebase already depends on Boost,
 but probably not yet on Boost Python. Can you point me to the relevant
 documentation, and maybe to an example? One of the problems I have is
 that the current wrapping is done auto-magically with SWIG and I would
 like to deviate the less possible from that patter.

 Thank you!

 Cheers,
 Daniele

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy

2014-10-14 Thread Daniele Nicolodi
On 14/10/14 13:39, Nathaniel Smith wrote:
 If the goal is to have something that works kind of like the new buffer
 protocol but with a wider variety of python versions, then you might
 find the old array interface useful:
 http://docs.scipy.org/doc/numpy/reference/arrays.interface.html
 
 I always get confused by the history here but I believe that that's the
 numpy-only interface that later got cleaned up and generalized to become
 the new buffer interface. Numpy itself still supports it.

Hello Nathaniel,

thanks for the pointer, that's what I need.

However, reading the documentation you linked, it looks like the new
buffer protocol is also available in Python 2.6, and I don't have the
need to support older versions, and this is confirmed by the python
documentation:

https://docs.python.org/2.7/c-api/buffer.html

So the issue is just on how to tell SWIG to implement this interface in
the wrappers it generates. Does anyone have a pointer to some relevant
documentation or examples?

Thanks!

Cheers,
Daniele

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy

2014-10-14 Thread John Zwinck
On Tue, Oct 14, 2014 at 6:19 PM, Daniele Nicolodi dani...@grinta.net wrote:
 On Mo, 2014-10-13 at 13:35 +0200, Daniele Nicolodi wrote:
  I have a C++ application that collects float, int or complex data in a
  possibly quite large std::vector. The application has some SWIG
  generated python wrappers that expose this vector to python. However,
  the standard way in which SWIG exposes the data is to create a touple
  and pass this to python, where it is very often converted to a numpy
  array for processing. Of course this is not efficient.

 Boost Python may be an option as the codebase already depends on Boost,
 but probably not yet on Boost Python. Can you point me to the relevant
 documentation, and maybe to an example? One of the problems I have is
 that the current wrapping is done auto-magically with SWIG and I would
 like to deviate the less possible from that patter.

Some time ago I needed to do something similar.  I fused the NumPy C
API and Boost.Python with a small bit of code which I then
open-sourced as part of a slightly larger library.  The most relevant
part for you is here:
https://github.com/jzwinck/pccl/blob/master/NumPyArray.hpp

In particular, it offers this function:

  // use already-allocated storage for the array
  // (it will not be initialized, and its format must match the given dtype)
  boost::python::object makeNumPyArrayWithData(
boost::python::list const dtype, unsigned count, void* data);

What is dtype?  For that you can use another small widget in my
library: https://github.com/jzwinck/pccl/blob/master/NumPyDataType.hpp

So the outline for your use case: you have a C array or C++ vector of
contiguous plain old data.  You create a NumPyDataType(), call
append() on it with your data type, then pass all of the above to
makeNumPyArrayWithData().  The result is a boost::python::object which
you can return from C++ to Python with zero copies made of your data.

Even if you don't decide to use Boost.Python, maybe you will find the
implementation instructive.  The functions I described take only a few
lines.

John Zwinck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy

2014-10-14 Thread Daniele Nicolodi
On 14/10/14 14:11, Daniele Nicolodi wrote:
 On 14/10/14 13:39, Nathaniel Smith wrote:
 If the goal is to have something that works kind of like the new buffer
 protocol but with a wider variety of python versions, then you might
 find the old array interface useful:
 http://docs.scipy.org/doc/numpy/reference/arrays.interface.html

 I always get confused by the history here but I believe that that's the
 numpy-only interface that later got cleaned up and generalized to become
 the new buffer interface. Numpy itself still supports it.
 
 Hello Nathaniel,
 
 thanks for the pointer, that's what I need.
 
 However, reading the documentation you linked, it looks like the new
 buffer protocol is also available in Python 2.6, and I don't have the
 need to support older versions, and this is confirmed by the python
 documentation:
 
 https://docs.python.org/2.7/c-api/buffer.html

I found one more problem: one of the data types I would like to expose
is a complex float array equivalent to numpy.complex64 dtype. However,
there is not such type format string defined in the python struct module,

Experimenting with numpy it seems like it extends the struct format
string definition with a new modifier 'Z' for complex:

 a = np.ones(10, dtype=np.complex64)
 m = memoryview(a)
 m.format
'Zf'

How standard is that?  Should I do the same in my extension?

Thank you.

Cheers,
Daniele

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy

2014-10-14 Thread Toby St Clere Smithe
John Zwinck jzwi...@gmail.com writes:
 Some time ago I needed to do something similar.  I fused the NumPy C
 API and Boost.Python with a small bit of code which I then
 open-sourced as part of a slightly larger library.  The most relevant
 part for you is here:
 https://github.com/jzwinck/pccl/blob/master/NumPyArray.hpp

There is also a 'Boost.NumPy', which is quite nice, though perhaps a bit
heavy-duty:

https://github.com/ndarray/Boost.NumPy/

Cheers,


Toby

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Changed behavior of np.gradient

2014-10-14 Thread Charles R Harris
On Sat, Oct 4, 2014 at 3:16 PM, Stéfan van der Walt ste...@sun.ac.za
wrote:

 On Oct 4, 2014 10:14 PM, Derek Homeier 
 de...@astro.physik.uni-goettingen.de wrote:
 
  +1 for an order=2 or maxorder=2 flag

 If you parameterize that flag, users will want to change its value (above
 two). Perhaps rather use a boolean flag such as second_order or
 high_order, unless it seems feasible to include additional orders in the
 future.

How about 'matlab={True, False}'. There is an open issue
https://github.com/numpy/numpy/issues/5184 for this and it would be good
to decide before 1.9.1 comes out.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Changed behavior of np.gradient

2014-10-14 Thread Nathaniel Smith
On 4 Oct 2014 22:17, Stéfan van der Walt ste...@sun.ac.za wrote:

 On Oct 4, 2014 10:14 PM, Derek Homeier 
de...@astro.physik.uni-goettingen.de wrote:
 
  +1 for an order=2 or maxorder=2 flag

 If you parameterize that flag, users will want to change its value (above
two). Perhaps rather use a boolean flag such as second_order or
high_order, unless it seems feasible to include additional orders in the
future.

Predicting the future is hard :-). And in particular high_order= would
create all kinds of confusion if in the future we added 3rd order
approximations but high_order=True continued to mean 2nd order because of
compatibility. I like maxorder (or max_order would be more pep8ish I guess)
because it leaves our options open. (Similar to how it's often better to
have a kwarg that can take two possible string values than to have a
boolean kwarg. It makes current code more explicit and makes future
enhancements easier.)

-n
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Changed behavior of np.gradient

2014-10-14 Thread Charles R Harris
On Tue, Oct 14, 2014 at 10:57 AM, Nathaniel Smith n...@pobox.com wrote:

 On 4 Oct 2014 22:17, Stéfan van der Walt ste...@sun.ac.za wrote:
 
  On Oct 4, 2014 10:14 PM, Derek Homeier 
 de...@astro.physik.uni-goettingen.de wrote:
  
   +1 for an order=2 or maxorder=2 flag
 
  If you parameterize that flag, users will want to change its value
 (above two). Perhaps rather use a boolean flag such as second_order or
 high_order, unless it seems feasible to include additional orders in the
 future.

 Predicting the future is hard :-). And in particular high_order= would
 create all kinds of confusion if in the future we added 3rd order
 approximations but high_order=True continued to mean 2nd order because of
 compatibility. I like maxorder (or max_order would be more pep8ish I guess)
 because it leaves our options open. (Similar to how it's often better to
 have a kwarg that can take two possible string values than to have a
 boolean kwarg. It makes current code more explicit and makes future
 enhancements easier.)


I think maxorder is a bit misleading. The both versions are second order in
the interior while at the ends the old is first order and the new is second
order. Maybe edge_order?

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Changed behavior of np.gradient

2014-10-14 Thread Nathaniel Smith
On 14 Oct 2014 18:29, Charles R Harris charlesr.har...@gmail.com wrote:



 On Tue, Oct 14, 2014 at 10:57 AM, Nathaniel Smith n...@pobox.com wrote:

 On 4 Oct 2014 22:17, Stéfan van der Walt ste...@sun.ac.za wrote:
 
  On Oct 4, 2014 10:14 PM, Derek Homeier 
de...@astro.physik.uni-goettingen.de wrote:
  
   +1 for an order=2 or maxorder=2 flag
 
  If you parameterize that flag, users will want to change its value
(above two). Perhaps rather use a boolean flag such as second_order or
high_order, unless it seems feasible to include additional orders in the
future.

 Predicting the future is hard :-). And in particular high_order= would
create all kinds of confusion if in the future we added 3rd order
approximations but high_order=True continued to mean 2nd order because of
compatibility. I like maxorder (or max_order would be more pep8ish I guess)
because it leaves our options open. (Similar to how it's often better to
have a kwarg that can take two possible string values than to have a
boolean kwarg. It makes current code more explicit and makes future
enhancements easier.)


 I think maxorder is a bit misleading. The both versions are second order
in the interior while at the ends the old is first order and the new is
second order. Maybe edge_order?

Ah, that makes sense. edge_order makes more sense to me too then - and we
can always add interior_order to complement it later, if appropriate.

The other thing to decide on is the default. Is the 2nd order version
generally preferred (modulo compatibility)? If so then it might make sense
to keep it the default, given that there are already numpy's in the wild
with that version, so we can't fully guarantee compatibility even if we
wanted to. But what do others think?

-n
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] numpy.mean slicing in a netCDF file

2014-10-14 Thread Stephan Hoyer
Hi Fadzil,

My strong recommendation is that you don't just use numpy/netCDF4 to
process your data, but rather use one of a multitude of packages that have
been developed specifically to facilitate working with labeled data from
netCDF files:
- Iris: http://scitools.org.uk/iris/
- CDAT: http://uvcdat.llnl.gov/
- xray (my project): http://xray.readthedocs.org

I can't answer your specific question without taking a careful look at your
data, but in very general terms, your code will have fewer bugs if you can
use meaningful labels to refer to your data rather than numeric ranges like
396:757:12.

Best,
Stephan


On Tue, Oct 14, 2014 at 3:50 AM, Fadzil Mnor fadzilmno...@gmail.com wrote:

 Hi all,
 I wrote a script and plot monthly mean zonal wind (from a netcdf file
 names uwnd.mon.mean.nc) and I'm not sure if I'm doing it correctly. What
 I have:


 
 *#this part calculates mean values for january only, from 1980-2010; thus
 the index looks like this 396:757:12*

 def ujan():
 f = nc.Dataset('~/data/ncep/uwnd.mon.mean.nc')
 u10_1 = f.variables['uwnd']
 u10_2 = np.mean(u10_1[396:757:12,:,38,39:43],axis=0)
 return u10_2

 uJan = ujan()* #calling function*

 *#this part is only to define lon, lat and level *
 q = nc.Dataset('~/data/ncep/uwnd.mon.mean.nc')
 lon=q.variables['lon']
 lat=q.variables['lat']
 lev=q.variables['level']

 *#for some reason I need to define this unless it gave error length of x
 must be number of column in z*

 lon=lon[39:43]

 *#begin plotting*

 clevs=np.arange(-10.,10.,0.5)
 fig = plt.figure(figsize=(11, 8))
 fig.clf()
 ax = fig.add_subplot(111)
 ax.axis([97.5, 105., 1000., 10.])
 ax.tick_params(direction='out', which='both')
 ax.set_xlabel('Lon (degrees)')
 ax.set_ylabel('Pressure (mb)')
 ax.set_xticks(np.arange(97.5, 105., .5))
 ax.set_yticks([1000, 700, 500, 300, 100, 10])
 cs=ax.contourf(lon, lev, uJan, clevs, extend='both',cmap='seismic')
 plt.title('Zonal winds average (Jan, 1981-2010)')
 cax = fig.add_axes([0.99, 0.1, 0.03, 0.8])
 aa=fig.colorbar(cs,cax=cax,orientation='vertical')
 aa.set_label('m/s')
 plt.savefig('~/uwind-crossection-test.png', bbox_inches='tight')

 ***

 the result is attached.
 I have no idea how to confirm the result (at least until this email is
 written) , but I believe the lower altitude values should be mostly
 negative, because over this region, the zonal wind are usually easterly
 (thus,negative values), but I got positive values.

 Put the information above aside, *I just want to know if my slicing in
 the ujan() function is correct*. If it is, then, there must be nothing
 wrong(except my above mentioned assumption).
 The data file dimension is:
 *[time,level,latitude,longitude]*

 This part:
 *u10_2 = np.mean(u10_1[396:757:12,:,38,39:43],axis=0)*
 The line above will calculate the mean of zonal wind (uwnd) in a range of
 time index 396 to 757 for each year (january only), for all vertical level,
 at latitude index 38 (5 N) and in between longitude index 39 to 43
 (97.5E-105E).
 I assume it will calculate a 30-year average of zonal wind for january
 only.
 Is this correct?

 Thank you.

 Fadzil

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] numpy.mean slicing in a netCDF file

2014-10-14 Thread Fadzil Mnor
Thank you Stephan,
I've been trying to install IRIS on my laptop (OS X) for months. Errors
everywhere.
I'll look at that IRIS again, and other links.

Cheers,


Fadzil

On Tue, Oct 14, 2014 at 7:09 PM, Stephan Hoyer sho...@gmail.com wrote:

 Hi Fadzil,

 My strong recommendation is that you don't just use numpy/netCDF4 to
 process your data, but rather use one of a multitude of packages that have
 been developed specifically to facilitate working with labeled data from
 netCDF files:
 - Iris: http://scitools.org.uk/iris/
 - CDAT: http://uvcdat.llnl.gov/
 - xray (my project): http://xray.readthedocs.org

 I can't answer your specific question without taking a careful look at
 your data, but in very general terms, your code will have fewer bugs if you
 can use meaningful labels to refer to your data rather than numeric ranges
 like 396:757:12.

 Best,
 Stephan


 On Tue, Oct 14, 2014 at 3:50 AM, Fadzil Mnor fadzilmno...@gmail.com
 wrote:

 Hi all,
 I wrote a script and plot monthly mean zonal wind (from a netcdf file
 names uwnd.mon.mean.nc) and I'm not sure if I'm doing it correctly. What
 I have:


 
 *#this part calculates mean values for january only, from 1980-2010; thus
 the index looks like this 396:757:12*

 def ujan():
 f = nc.Dataset('~/data/ncep/uwnd.mon.mean.nc')
 u10_1 = f.variables['uwnd']
 u10_2 = np.mean(u10_1[396:757:12,:,38,39:43],axis=0)
 return u10_2

 uJan = ujan()* #calling function*

 *#this part is only to define lon, lat and level *
 q = nc.Dataset('~/data/ncep/uwnd.mon.mean.nc')
 lon=q.variables['lon']
 lat=q.variables['lat']
 lev=q.variables['level']

 *#for some reason I need to define this unless it gave error length of x
 must be number of column in z*

 lon=lon[39:43]

 *#begin plotting*

 clevs=np.arange(-10.,10.,0.5)
 fig = plt.figure(figsize=(11, 8))
 fig.clf()
 ax = fig.add_subplot(111)
 ax.axis([97.5, 105., 1000., 10.])
 ax.tick_params(direction='out', which='both')
 ax.set_xlabel('Lon (degrees)')
 ax.set_ylabel('Pressure (mb)')
 ax.set_xticks(np.arange(97.5, 105., .5))
 ax.set_yticks([1000, 700, 500, 300, 100, 10])
 cs=ax.contourf(lon, lev, uJan, clevs, extend='both',cmap='seismic')
 plt.title('Zonal winds average (Jan, 1981-2010)')
 cax = fig.add_axes([0.99, 0.1, 0.03, 0.8])
 aa=fig.colorbar(cs,cax=cax,orientation='vertical')
 aa.set_label('m/s')
 plt.savefig('~/uwind-crossection-test.png', bbox_inches='tight')

 ***

 the result is attached.
 I have no idea how to confirm the result (at least until this email is
 written) , but I believe the lower altitude values should be mostly
 negative, because over this region, the zonal wind are usually easterly
 (thus,negative values), but I got positive values.

 Put the information above aside, *I just want to know if my slicing in
 the ujan() function is correct*. If it is, then, there must be nothing
 wrong(except my above mentioned assumption).
 The data file dimension is:
 *[time,level,latitude,longitude]*

 This part:
 *u10_2 = np.mean(u10_1[396:757:12,:,38,39:43],axis=0)*
 The line above will calculate the mean of zonal wind (uwnd) in a range of
 time index 396 to 757 for each year (january only), for all vertical level,
 at latitude index 38 (5 N) and in between longitude index 39 to 43
 (97.5E-105E).
 I assume it will calculate a 30-year average of zonal wind for january
 only.
 Is this correct?

 Thank you.

 Fadzil

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion



 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] array indexing question

2014-10-14 Thread Neal Becker
I'm using np.nonzero to construct the tuple:
(array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2]), array([1, 3, 5, 7, 2, 3, 6, 7, 4, 
5, 6, 7]))

Now what I want is the 2-D index array:

[1,3,5,7,
2,3,6,7,
4,5,6,7]

Any ideas?

-- 
-- Those who don't understand recursion are doomed to repeat it

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] array indexing question

2014-10-14 Thread Nathaniel Smith
For this to work at all you have to know a priori that there are the same
number of non-zero entries in each row of your mask. Given that you know
that, isn't it just a matter of calling reshape on the second array?
On 14 Oct 2014 20:37, Neal Becker ndbeck...@gmail.com wrote:

 I'm using np.nonzero to construct the tuple:
 (array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2]), array([1, 3, 5, 7, 2, 3, 6,
 7, 4,
 5, 6, 7]))

 Now what I want is the 2-D index array:

 [1,3,5,7,
 2,3,6,7,
 4,5,6,7]

 Any ideas?

 --
 -- Those who don't understand recursion are doomed to repeat it

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Changed behavior of np.gradient

2014-10-14 Thread Charles R Harris
On Tue, Oct 14, 2014 at 11:50 AM, Nathaniel Smith n...@pobox.com wrote:

 On 14 Oct 2014 18:29, Charles R Harris charlesr.har...@gmail.com
 wrote:
 
 
 
  On Tue, Oct 14, 2014 at 10:57 AM, Nathaniel Smith n...@pobox.com wrote:
 
  On 4 Oct 2014 22:17, Stéfan van der Walt ste...@sun.ac.za wrote:
  
   On Oct 4, 2014 10:14 PM, Derek Homeier 
 de...@astro.physik.uni-goettingen.de wrote:
   
+1 for an order=2 or maxorder=2 flag
  
   If you parameterize that flag, users will want to change its value
 (above two). Perhaps rather use a boolean flag such as second_order or
 high_order, unless it seems feasible to include additional orders in the
 future.
 
  Predicting the future is hard :-). And in particular high_order= would
 create all kinds of confusion if in the future we added 3rd order
 approximations but high_order=True continued to mean 2nd order because of
 compatibility. I like maxorder (or max_order would be more pep8ish I guess)
 because it leaves our options open. (Similar to how it's often better to
 have a kwarg that can take two possible string values than to have a
 boolean kwarg. It makes current code more explicit and makes future
 enhancements easier.)
 
 
  I think maxorder is a bit misleading. The both versions are second order
 in the interior while at the ends the old is first order and the new is
 second order. Maybe edge_order?

 Ah, that makes sense. edge_order makes more sense to me too then - and we
 can always add interior_order to complement it later, if appropriate.

 The other thing to decide on is the default. Is the 2nd order version
 generally preferred (modulo compatibility)? If so then it might make sense
 to keep it the default, given that there are already numpy's in the wild
 with that version, so we can't fully guarantee compatibility even if we
 wanted to. But what do others think?

I'd be inclined to keep the older as the default and regard adding the
keyword as a bugfix. I should have caught the incompatibility in review.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion