Re: [Numpy-discussion] Best way to expose std::vector to be used with numpy
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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