[ncl-talk] Discrepancies between Python and NCL results

[Guido Cioni]

As a matter of fact it was, yes.

Slicing in Python works differently so that using 1:150 would yield an array with vertical dimension 149 and not 150. :) 

> Il giorno 21 feb 2018, alle ore 19:25, Bill Ladwig <ladwig at ucar.edu> ha scritto:
> 
> Hi Guido,
> 
> Just at a quick glance, I see this in your NCL script:
> 
> w=data->w(:,1:150,ind_dry)
> 
> In Python, you have this:
> 
>  w=np.array(nc2.variables["w"][:][:,0:150,ind_dry])
> 
> You've got difference slices going (both NCL and Python use 0-based indexes), was that intentional?
> 
> Bill
>  
> 
> On Wed, Feb 21, 2018 at 2:58 AM, Guido Cioni <guidocioni at gmail.com <mailto:guidocioni at gmail.com>> wrote:
> Hey all,
> I have a weird question addressed to experts of both NCL and Python. 
> 
> In my effort to speed up a NCL script which computes the Moist Static Energy of updraft points and of the environment in a LES model simulation of diurnal convection, I've tried to translate the entire script into Python. I was interested to see whether Python was faster than NCL: long story short, it's not :) The problem is that my data is really big (72 x 150 x 1280000 points) so that a single assignment a = file->var takes approx. 40 minutes. A single plot takes about 3 hours to be produced. And I'm using no single loop! :) 
> 
> But let's come to the problem. I have to select updrafts, which I defined as grid points having a vertical velocity of at least 1 m/s and a cloud water content of at least 1E-4 kg/kg. For these points I want to construct a vertical profile of moist static energy. Thus, it seems to me that using the where function was going to be my best shot. This is what I tried to do to select the grid points:
> 
>    w_cloud=where(qc.gt.th_qc .and. w.ge.1, mse, mse at _FillValue)
>    w_envir=where(qc.le.th_qc .or. w.lt.1, mse, mse at _FillValue)
> 
> And then averaged over time and space to get a single vertical profile (I know that NCL will take care of the missing values using these operations). The result is not bad..
> 
> 
> The red line shows the difference w_cloud-w_envir which is, as expected, positive meaning that updrafts are moister and warmer than the surrounding. 
> 
> Coming to Python I decided to use a similar approach with numpy masked arrays
> 
> mse_cloud_3D=np.ma.masked_where( ((w<1) | (qc<1E-04)), mse)
> mse_cloud=np.ma.mean(mse_cloud_3D, axis=(0,2))
> 
> mse_envir_3D=np.ma.masked_where( ((w>=1) & (qc>=1E-04)), mse)
> mse_envir=np.ma.mean(mse_envir_3D, axis=(0,2))
> 
> For the average MSE the profile looks almost exactly the same (good!) but the difference instead is completely different, also considering the order of magnitude.
> 
> So my question is...what I'm doing wrong? How could the average MSE be fine and the difference between the updraft and the environment not? I would tend to trust more NCL as the masked arrays of numpy are behaving pretty weird sometimes.
> 
> I'm attaching the scripts cleaned as much as possible, although it is impossible (for obvious reasons) to send the data along :) 
> 
> If anyone has opinion/comments...
> 
> Cheers 
> 
> 
> 
> Guido Cioni
> http://guidocioni.altervista <http://guidocioni.altervista/>.org
> 
> 
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk <http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
> 
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180221/63b961df/attachment.html>