![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
Changeset 2319 in lmdz_wrf for trunk/tools
- Timestamp:
- Feb 5, 2019, 5:04:56 PM (6 years ago)
- Location:
- trunk/tools
- Files:
-
- 2 edited
-
module_scientific.f90 (modified) (3 diffs)
-
nc_var_tools.py (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/module_scientific.f90
r2318 r2319 5663 5663 INTEGER :: ii3, ss1, ss2, ss3 5664 5664 INTEGER :: Ncounts, Nin 5665 INTEGER, DIMENSION(1) :: dmaxvarin 5665 5666 CHARACTER(len=3) :: val1S, val2S 5666 5667 CHARACTER(len=30) :: val3S … … 5699 5700 ss3 = 3 + 1 5700 5701 ii3 = 1 + 1 5702 5703 dmaxvarin = UBOUND(varin) 5701 5704 5702 5705 ! Let's be efficient? … … 5749 5752 ELSE 5750 5753 i1 = gridsin(s1,s2,s3,1,idv) 5751 varout(s1,s2,s3,1) = varin(i1) 5752 varout(s1,s2,s3,2) = varin(i1) 5753 varout(s1,s2,s3,3) = varin(i1) 5754 varout(s1,s2,s3,4) = varin(i1)*varin(i1) 5755 varout(s1,s2,s3,5) = zeroRK 5756 varout(s1,s2,s3,6) = varin(i1) 5757 varout(s1,s2,s3,7) = Nin*1. 5754 IF (i1 > 0 .AND. i1 <= dmaxvarin(1)) THEN 5755 varout(s1,s2,s3,1) = varin(i1) 5756 varout(s1,s2,s3,2) = varin(i1) 5757 varout(s1,s2,s3,3) = varin(i1) 5758 varout(s1,s2,s3,4) = varin(i1)*varin(i1) 5759 varout(s1,s2,s3,5) = zeroRK 5760 varout(s1,s2,s3,6) = varin(i1) 5761 varout(s1,s2,s3,7) = Nin*1. 5762 END IF 5758 5763 END IF 5759 5764 END DO -
trunk/tools/nc_var_tools.py
r2315 r2319 27873 27873 newlonbndsv = np.zeros((4,sdimy,sdimx), dtype=np.float) 27874 27874 newlatbndsv = np.zeros((4,sdimy,sdimx), dtype=np.float) 27875 print 'Lluis shapes olonbnds1D:', olonbnds1D.shape27876 27875 for j in range(sdimy): 27877 27876 for i in range(sdimx): … … 28404 28403 getvarybndst = ydimvarbnds2D.transpose() 28405 28404 28406 Ngridsint, gridsint, areas t, percenst =\28405 Ngridsint, gridsint, areas2Dt, areast, percenst = \ 28407 28406 fsci.module_scientific.grid_spacepercen(xcavals=reflont, \ 28408 28407 ycavals=reflatt, xbavals=refvarxbndst, ybavals=refvarybndst, \ … … 28420 28419 gridsin = gridsin - 1 28421 28420 28421 areas2D = areas2Dt.transpose() 28422 28422 areas = areast.transpose() 28423 28423 percens = percenst.transpose() … … 28497 28497 28498 28498 gridsin = np.zeros((2,Ngridsinmax,1,Nslices), dtype=int) 28499 areas = np.ones( tuple(varfinalshape), dtype=np.float)*gen.fillValueF28499 areas = np.ones((Ngridsinmax,1,Nslices), dtype=np.float)*gen.fillValueF 28500 28500 percens = np.zeros((Ngridsinmax,1,Nslices), dtype=np.float) 28501 28501 … … 28538 28538 gridsin[1,0:Nt,0,islc] = ainds[0:Nt,0] 28539 28539 gridsin[0,0:Nt,0,islc] = ainds[0:Nt,1] 28540 areas[0:Nt,0,islc] = 1. 28540 28541 percens[0:Nt,0,islc] = 1. 28541 for iv in range(Nt):28542 i = ainds[iv,0]28543 j = ainds[iv,1]28544 areas[j,i] = 1.28545 28542 28546 28543 # Values for file … … 28632 28629 innewvar.setncattr('coordinates',' '.join(Scgrid[::-1])) 28633 28630 28634 anewvar = onewnc.createVariable(dn+'gridarea','f',tuple( varfinaldims),\28631 anewvar = onewnc.createVariable(dn+'gridarea','f',tuple(Spgrid), \ 28635 28632 fill_value = gen.fillValueF) 28636 28633 basicvardef(anewvar ,dn+'gridarea', "area of the grids cells from " + \ 28637 28634 dn, vu) 28638 28635 anewvar.setncattr('coordinates',' '.join(varfinaldims[::-1])) 28639 anewvar[:] = areas[ :]28636 anewvar[:] = areas[0:Ngridsinmax,...] 28640 28637 28641 28638 aanewvar = onewnc.createVariable(dn+'area','f',tuple(Srgrid), \ … … 28654 28651 j = 0 28655 28652 for i in range(Nslices): 28653 slicearea = 0. 28656 28654 innewvar[:,0:Ngridsin[j,i],i] = gridsin[:,0:Ngridsin[j,i],j,i] 28657 28655 pnewvar[0:Ngridsin[j,i],i]= percens[0:Ngridsin[j,i],j,i] 28658 slicearea = 0. 28659 for iv in range(Ngridsin[j,i]): 28660 ix = gridsin[1,iv,j,i] 28661 iy = gridsin[0,iv,j,i] 28662 slicearea = slicearea + areas[iy,ix]*percens[iv,j,i] 28656 slicearea = np.sum(areas[:,j,i]*percens[:,j,i]) 28663 28657 aanewvar[i]= slicearea 28664 28658 elif len(dv.shape) == 2: … … 28689 28683 ag = np.abs(fsci.module_scientific.shoelace_area_polygon(nvertex=4, poly=apolygon)) 28690 28684 slicearea = slicearea + ag*percens[iv,j,i] 28691 aa = aa + areas[i y,ix]28685 aa = aa + areas[iv,j,i] 28692 28686 pp = pp + percens[iv,j,i] 28693 28687 if np.isnan(slicearea): … … 28699 28693 ix = gridsin[1,iv,j,i] 28700 28694 iy = gridsin[0,iv,j,i] 28701 print iv, iy, ix, areas[i y,ix], percens[j,i]28695 print iv, iy, ix, areas[iv,j,i], percens[iv,j,i] 28702 28696 quit(-1) 28703 28697 aanewvar[j,i]= slicearea … … 28835 28829 for iv in range(Nin): 28836 28830 pA = oslicepA[inpointsA[iv,jB,iB,jA,iA],jA,iA] 28831 aA = oslicepA[inpointsA[iv,jB,iB,jA,iA],jA,iA] 28837 28832 pB = oslicepB[inpointsB[iv,jB,iB,jA,iA],jB,iB] 28833 aB = oslicepB[inpointsB[iv,jB,iB,jA,iA],jB,iB] 28838 28834 pnewvar[iv,jB,iB,jA,iA]= pA*pB 28839 28835 ixA = sliceinA[1,inpointsA[iv,jB,iB,jA,iA],jA,iA] 28840 28836 iyA = sliceinA[0,inpointsA[iv,jB,iB,jA,iA],jA,iA] 28841 aA = osliceaA[iyA,ixA]28842 28837 ixB = sliceinB[1,inpointsB[iv,jB,iB,jA,iA],jB,iB] 28843 28838 iyB = sliceinB[0,inpointsB[iv,jB,iB,jA,iA],jB,iB] 28844 aB = osliceaB[iyB,ixB]28845 28839 anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB) 28846 28840 slicearea = slicearea + (aA*pA)*(aB*pB) … … 28918 28912 28919 28913 # Getting respective A & B 28920 for jc in range(dyC):28921 for ic in range(dxC):28922 print jc, ic, 'Lluis shapes NpointsABt[jc,ic,:,:]', NpointsABt[jc,ic,:,:].shape, \28923 'pointsABt[jc,ic,:,:]', NpointsABt[jc,ic,:,:].shape, 'sliceNAt', sliceNAt.shape, \28924 'sliceinAt', sliceinAt.shape, 'A dx dy dxy:', dxA, dyA, dxyAB, \28925 'B dx dy dxy:', dxA, dyA, dxA28926 28927 NptsAt, inptsAt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \28928 npointsa=NpointsABt[jc,ic,:,:], pointsa=pointsABt[jc,ic,:,:], npointsb=sliceNAt,\28929 pointsb=sliceinAt, dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxA, dyb=dyA, dxyb=dxyA)28930 NpointsA[jc,ic,j,i,:,:] = NptsAt.transpose()28931 inpointsA[:,jc,ic,j,i,:,:] = inptsAt.transpose()28932 28933 NptsBt, inptsBt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \28934 npointsa=NpointsABt[jc,ic,:,:], pointsa=pointsABt[jc,ic,:,:], npointsb=sliceNBt,\28935 pointsb=sliceinBt, dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxB, dyb=dyB, dxyb=dxyB)28936 NpointsB[jc,ic,j,i,:,:] = NptsAt.transpose()28937 inpointsB[:,jc,ic,j,i,:,:] = inptsAt.transpose()28914 #for jc in range(dyC): 28915 # for ic in range(dxC): 28916 # print jc, ic, 'Lluis shapes NpointsABt[jc,ic,:,:]', NpointsABt[jc,ic,:,:].shape, \ 28917 # 'pointsABt[jc,ic,:,:]', NpointsABt[jc,ic,:,:].shape, 'sliceNAt', sliceNAt.shape, \ 28918 # 'sliceinAt', sliceinAt.shape, 'A dx dy dxy:', dxA, dyA, dxyAB, \ 28919 # 'B dx dy dxy:', dxA, dyA, dxA 28920 # 28921 # NptsAt, inptsAt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \ 28922 # npointsa=NpointsABt[jc,ic,:,:], pointsa=pointsABt[jc,ic,:,:], npointsb=sliceNAt,\ 28923 # pointsb=sliceinAt, dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxA, dyb=dyA, dxyb=dxyA) 28924 # NpointsA[jc,ic,j,i,:,:] = NptsAt.transpose() 28925 # inpointsA[:,jc,ic,j,i,:,:] = inptsAt.transpose() 28926 # 28927 # NptsBt, inptsBt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \ 28928 # npointsa=NpointsABt[jc,ic,:,:], pointsa=pointsABt[jc,ic,:,:], npointsb=sliceNBt,\ 28929 # pointsb=sliceinBt, dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxB, dyb=dyB, dxyb=dxyB) 28930 # NpointsB[jc,ic,j,i,:,:] = NptsAt.transpose() 28931 # inpointsB[:,jc,ic,j,i,:,:] = inptsAt.transpose() 28938 28932 28939 28933 # Remembering that it is python (C-like...) 28940 inpointsAB = inpointsAB-128941 inpointsC = inpointsC-128934 inpointsAB[:,:,:,j,i,:,:] = inpointsAB[:,:,:,j,i,:,:]-1 28935 inpointsC[:,:,:,j,i,:,:] = inpointsC[:,:,:,j,i,:,:]-1 28942 28936 28943 28937 maxNpointsABC = np.max(NpointsABC) … … 28950 28944 28951 28945 newvar = onewnc.createVariable(dn+'Ngrid','i', tuple(Srgrid)) 28952 print 'Lluis shapes: newvar', newvar.shape, 'NpointsABC:', NpointsABC.shape, 'Srgrid:', Srgrid28953 28946 newvar[:] = NpointsABC[:] 28954 28947 basicvardef(newvar, dn+'Ngrid', "number of grids cells from slice " + \ … … 28979 28972 "grids cells from " + newslcvarns[0] + " laying within " + newslcvarns[1], '1') 28980 28973 pnewvar.setncattr('coordinates',' '.join(Spgrid[::-1])) 28981 28982 print 'Lluis shapes NpointsANC:', NpointsABC.shape, 'pointsABC:', pointsABC.shape, \28983 'inpoints AB:', inpointsAB.shape, 'sliceinAB:', sliceinAB.shape, 'oslicepAB:', \28984 oslicepAB.shape, 'osliceaAB:', osliceaAB.shape, 'inpoints C:', inpointsC.shape, \28985 'sliceinC:', sliceinC.shape, 'oslicepC:', oslicepC.shape, 'osliceaC:', osliceaC.shape28986 28974 ixA = -1 28987 28975 iyA = -1 … … 28997 28985 for iC in range(dxC): 28998 28986 Nin = NpointsABC[jC,iC,jB,iB,jA,iA] 28999 innewvar[:,0:Nin,jC,iC,jB,iB,jA,iA] = pointsABC[:,0:Nin,jC,iC,jB,iB,jA,iA] 28987 innewvar[:,0:Nin,jC,iC,jB,iB,jA,iA] = \ 28988 pointsABC[:,0:Nin,jC,iC,jB,iB,jA,iA] 29000 28989 aanewvar[jC,iC,jB,iB,jA,iA]= gen.fillValueF 29001 28990 anewvar[:,jC,iC,jB,iB,jA,iA]= gen.fillValueF … … 29003 28992 slicearea = 0. 29004 28993 for iv in range(Nin): 29005 pA = oslicepAB[inpointsAB[iv,jC,iC,jB,iB,jA,iA],jB,iB,jA,iA] 29006 pB = oslicepC[inpointsC[iv,jC,iC,jB,iB,jA,iA],jC,iC] 28994 ivAB = inpointsAB[iv,jC,iC,jB,iB,jA,iA] 28995 ivC = inpointsC[iv,jC,iC,jB,iB,jA,iA] 28996 pA = oslicepAB[ivAB,jB,iB,jA,iA] 28997 pB = oslicepC[ivC,jC,iC] 28998 aA = osliceaAB[ivAB,jB,iB,jA,iA] 28999 aB = osliceaC[ivC,jC,iC] 29000 ixAB = sliceinAB[1,ivAB,jB,iB,jA,iA] 29001 iyAB = sliceinAB[0,ivAB,jB,iB,jA,iA] 29002 ixC = sliceinC[1,ivC,jC,iC] 29003 iyC = sliceinC[0,ivC,jC,iC] 29004 anewvar[iv,jC,iC,jB,iB,jA,iA]= (aA*pA)*(aB*pB) 29007 29005 pnewvar[iv,jC,iC,jB,iB,jA,iA]= pA*pB 29008 ixAB = sliceinAB[1,inpointsAB[iv,jC,iC,jB,iB,jA,iA],jB,iB,jA,iA]29009 iyAB = sliceinAB[0,inpointsAB[iv,jC,iC,jB,iB,jA,iA],jB,iB,jA,iA]29010 ixA = inpointsA[1,NpointsA[iv,jC,iC,jB,iB,jA,iA],jC,iC,jB,iB,jA,iA]29011 iyA = inpointsA[0,NpointsA[iv,jC,iC,jB,iB,jA,iA],jC,iC,jB,iB,jA,iA]29012 ixB = inpointsB[1,NpointsB[iv,jC,iC,jB,iB,jA,iA],jC,iC,jB,iB,jA,iA]29013 iyB = inpointsB[0,NpointsB[iv,jC,iC,jB,iB,jA,iA],jC,iC,jB,iB,jA,iA]29014 ixC = sliceinC[1,inpointsC[iv,jC,iC,jB,iB,jA,iA],jC,iC]29015 iyC = sliceinC[0,inpointsC[iv,jC,iC,jB,iB,jA,iA],jC,iC]29016 print jA,iA,jB,iB,jC,iC,':',iyA,ixA,iyB,ixB,iyC,ixC29017 aA = osliceaAB[iyB,ixB,iyA,ixA]29018 aB = osliceaC[iyC,ixC]29019 anewvar[iv,jC,iC,jB,iB,jA,iA]= (aA*pA)*(aB*pB)29020 29006 slicearea = slicearea + (aA*pA)*(aB*pB) 29007 if np.isnan(slicearea): 29008 print errormsg 29009 print ' ' + fname + ': Nan slice area for slice:', \ 29010 jC, iC, jB, iB, jA, iA, 'N grids:', Nin,' !!' 29011 print ' #grid AB #grid grid_area ' + \ 29012 'grid_percen C # grid grid_area grid_percen____' 29013 for iv in range(Nin): 29014 print iv, inpointsAB[iv,jC,iC,jB,iB,jA,iA], \ 29015 aA, pA, inpointsC[iv,jC,iC,jB,iB,jA,iA], \ 29016 aB, pB 29017 quit(-1) 29021 29018 aanewvar[jC,iC,jB,iB,jA,iA]= slicearea 29022 29019 onewnc.sync() … … 29046 29043 pnewvar.setncattr('coordinates',' '.join(Spgrid[::-1])) 29047 29044 29048 (dyA, dxA, dyB, dxB) = sN.shape 29049 for jA in range(dyA): 29050 for iA in range(dxA): 29051 for jB in range(dyB): 29052 for iB in range(dxB): 29053 Nin = sN[jA,iA,jB,iB] 29054 if Nin > 0 and aanewvar[jA,iA,jB,iB] != 0.: 29055 sgpa[0:Nin,jA,iA,jB,iB] = anewvar[0:Nin,jA,iA,jB,iB]/ \ 29056 aanewvar[jA,iA,jB,iB] 29057 panewvar[0:Nin,jA,iA,jB,iB] = sgpa[0:Nin,jA,iA,jB,iB] 29058 # Let's check 29059 psum = np.sum(sgpa[0:Nin,jA,iA,jB,iB]) 29060 if not np.isnan(psum): 29061 if not gen.almost_zero(psum, 1., 4): 29062 print warnmsg 29063 print ' ' + fname + 'at slide:', jA,iA,jB,iB, \ 29064 'percentages do NOT add one down to 0.0001 !!' 29065 print ' Ngrids:', Nin, 'slice area:', \ 29066 aanewvar[jA,iA,jB,iB], 'total sum:', psum 29067 else: 29068 print errormsg 29069 print ' ' + fname + ': grid Nan slice area for final '+ \ 29070 ' slice:', jA, iA, jB, iB, 'N grids:', Nin,' !!' 29071 print ' slice total area:', aanewvar[jA,iA,jB,iB] 29072 print ' #grid grid_area slice_area grid_percen ______' 29073 for iv in range(Nin): 29074 print iv, anewvar[iv,jA,iA,jB,iB], \ 29075 aanewvar[jA,iA,jB,iB], sgpa[iv,jA,iA,jB,iB] 29076 quit(-1) 29045 if len(sN.shape)/2 == 2: 29046 (dyA, dxA, dyB, dxB) = sN.shape 29047 for jA in range(dyA): 29048 for iA in range(dxA): 29049 for jB in range(dyB): 29050 for iB in range(dxB): 29051 Nin = sN[jA,iA,jB,iB] 29052 if Nin > 0 and aanewvar[jA,iA,jB,iB] != 0.: 29053 sgpa[0:Nin,jA,iA,jB,iB] = anewvar[0:Nin,jA,iA,jB,iB]/ \ 29054 aanewvar[jA,iA,jB,iB] 29055 panewvar[0:Nin,jA,iA,jB,iB] = sgpa[0:Nin,jA,iA,jB,iB] 29056 # Let's check 29057 psum = np.sum(sgpa[0:Nin,jA,iA,jB,iB]) 29058 if not np.isnan(psum): 29059 if not gen.almost_zero(psum, 1., 4): 29060 print warnmsg 29061 print ' ' + fname + 'at slide:', jA,iA,jB,iB, \ 29062 'percentages do NOT add one down to 0.0001 !!' 29063 print ' Ngrids:', Nin, 'slice area:', \ 29064 aanewvar[jA,iA,jB,iB], 'total sum:', psum 29065 else: 29066 print errormsg 29067 print ' ' + fname + ': grid Nan slice area for final '+ \ 29068 ' slice:', jA, iA, jB, iB, 'N grids:', Nin,' !!' 29069 print ' slice total area:', aanewvar[jA,iA,jB,iB] 29070 print ' #grid grid_area slice_area grid_percen ______' 29071 for iv in range(Nin): 29072 print iv, anewvar[iv,jA,iA,jB,iB], \ 29073 aanewvar[jA,iA,jB,iB], sgpa[iv,jA,iA,jB,iB] 29074 quit(-1) 29075 elif len(sN.shape)/2 == 3: 29076 (dyA, dxA, dyB, dxB, dyC, dxC) = sN.shape 29077 for jA in range(dyA): 29078 for iA in range(dxA): 29079 for jB in range(dyB): 29080 for iB in range(dxB): 29081 for jC in range(dyC): 29082 for iC in range(dxC): 29083 Nin = sN[jA,iA,jB,iB,jC,iC] 29084 if Nin > 0 and aanewvar[jA,iA,jB,iB,jC,iC] != 0.: 29085 sgpa[0:Nin,jA,iA,jB,iB,jC,iC] = \ 29086 anewvar[0:Nin,jA,iA,jB,iB,jC,iC]/ \ 29087 aanewvar[jA,iA,jB,iB,jC,iC] 29088 panewvar[0:Nin,jA,iA,jB,iB,jC,iC] = \ 29089 sgpa[0:Nin,jA,iA,jB,iB,jC,iC] 29090 # Let's check 29091 psum = np.sum(sgpa[0:Nin,jA,iA,jB,iB,jC,iC]) 29092 if not np.isnan(psum): 29093 if not gen.almost_zero(psum, 1., 4): 29094 print warnmsg 29095 print ' ' + fname + 'at slide:', \ 29096 jA,iA,jB,iB,jC,iC, 'percentages do ' + \ 29097 'NOT add one down to 0.0001 !!' 29098 print ' Ngrids:', Nin, 'slice area:', \ 29099 aanewvar[jA,iA,jB,iB,jC,iC], \ 29100 'total sum:', psum 29101 else: 29102 print errormsg 29103 print ' ' + fname + ': grid Nan slice ' + \ 29104 'area for final slice:', jA, iA, jB, iB, \ 29105 jC, iC, 'N grids:', Nin,' !!' 29106 print ' slice total area:', \ 29107 aanewvar[jA,iA,jB,iB,jC,iC] 29108 print ' #grid grid_area slice_area ' + \ 29109 'grid_percen ______' 29110 for iv in range(Nin): 29111 print iv, anewvar[iv,jA,iA,jB,iB,jC,iC], \ 29112 aanewvar[jA,iA,jB,iB,jC,iC], \ 29113 sgpa[iv,jA,iA,jB,iB,jC,iC] 29114 quit(-1) 29115 else: 29116 print errormsg 29117 print ' ' + fname + ": rank ", len(sN.shape)/2,' of final slice not ready !!' 29118 quit(-1) 29119 29077 29120 onewnc.sync() 29078 29121 … … 29102 29145 29103 29146 print ' ' + infmsg 29104 print ' final slice:', Ndims, 's ape:', sN.shape29147 print ' final slice:', Ndims, 'shape:', sN.shape 29105 29148 29106 29149 # Checking 29107 i = 329150 i = 2 29108 29151 j = 5 29109 29152 k = 8 … … 29157 29200 if not gen.searchInlist(onewnc.dimensions,dnn): 29158 29201 add_dims(onc, onewnc, [dnn]) 29159 vardims.append(dnn)29160 varshape.append(len(onc.dimensions[dnn]))29161 29202 if not gen.searchInlist(onewnc.variables,vdimn): add_vars(onc, \ 29162 29203 onewnc, [vdimn]) … … 29166 29207 newvarshape = varshape + shapeslice 29167 29208 29168 # newvarslice = [] 29169 # for dimn in 29209 # Generic newvar slice 29210 newvarslice = [] 29211 for idd in range(len(newvarshape)): 29212 newvarslice.append(slice(0,newvarshape[idd])) 29170 29213 29171 29214 onewvars = {} … … 29193 29236 # Checking for size... 29194 29237 Nvdims = len(vshape) 29238 newvarslc = {} 29195 29239 varslc = {} 29196 29240 if Nvdims == 3: 29197 vSIZE = np.prod( varv.shape[:])*8.29241 vSIZE = np.prod(ovar.shape[:])*8. 29198 29242 if vSIZE > gen.oneGB*3.: 29199 29243 print ' ' + infmsg 29200 29244 print ' variable size: ', vSIZE / gen.oneGB,' larger than 3GB !!' 29201 29245 print ' splitting time evolution to compute' 29202 vHSIZE = np.prod( varv.shape[1:3])*8.29203 print ' var. horizontal size:', vHSIZE, 'dimt:', varv.shape[0]29246 vHSIZE = np.prod(ovar.shape[1:3])*8. 29247 print ' var. horizontal size:', vHSIZE, 'dimt:', ovar.shape[0] 29204 29248 tsteps = gen.oneGB*3. / vHSIZE 29205 tslices = range(0, varv.shape[2],tsteps)29206 print ' time-slices:', t klices29249 tslices = range(0,ovar.shape[0],tsteps) 29250 print ' time-slices:', tslices 29207 29251 for itt in tslcices[1,len(tslices)]: 29208 29252 iit = tslices[itt-1] … … 29210 29254 varslc[itt-1] = [slice(iit,eet), slice(0,vshape[1]), \ 29211 29255 slice(0,vbshape[2])] 29256 newvarslc[itt-1] = [slice(iit,eet)] + newvarslice[1:3] 29212 29257 else: 29213 tslices = range(0, varv.shape[0])29258 tslices = range(0,ovar.shape[0]) 29214 29259 varslc[0] = [slice(0,vshape[0]), slice(0,vshape[1]), \ 29215 29260 slice(0,vshape[2])] 29261 newvarslc[0] = newvarslice 29262 else: 29263 slicevv = [] 29264 for iid in range(Nvdims): 29265 slicevv.append(slice(0,vshape[0])) 29266 varslc[0] = slicevv 29267 newvarslc[0] = newvarslice 29268 29269 if Nvdims == 3: 29270 d0 = ovar.shape[0] 29271 d1 = ovar.shape[1] 29272 d2 = ovar.shape[2] 29273 elif Nvdims == 2: 29274 d0 = ovar.shape[0] 29275 d1 = ovar.shape[1] 29276 elif Nvdims == 1: 29277 d0 = ovar.shape[0] 29278 if slicespacedim.count(vdimns[0]) != 1: 29279 print errormsg 29280 print ' ' + fname + ": there is no way to be able to compute " + \ 29281 "sliced statistics for variable '" + varn + "' because it " + \ 29282 "does not have any horizontal spatial dimension:", slicespacedim, \ 29283 "!!" 29284 quit(-1) 29285 ind = slicespacedim.index(vdimns[0]) 29286 else: 29287 print errormsg 29288 print ' ' + fame + ": rank ", Nvdims, " of variable '" + varn + \ 29289 "' not ready !!" 29290 quit(-1) 29291 29292 print ' Lluis computing statistics ...' 29216 29293 29217 29294 for itt in varslc.keys(): 29295 nslcv = newvarslc[itt] 29218 29296 slcv = varslc[itt] 29219 29297 varv = ovar[tuple(slcv)] 29220 d0 = varv.shape[0] 29221 d1 = varv.shape[1] 29222 d2 = varv.shape[2] 29298 29223 29299 varvt = varv.transpose() 29224 29300 if Nvdims == 3 and len(sN.shape) == 4: … … 29232 29308 di1=d2, di2=d1, di3=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0], \ 29233 29309 maxngridsin=sin.shape[1]) 29310 elif Nvdims == 2 and len(sN.shape) == 3: 29311 voutt=fsci.module_scientific.multi_spaceweightstats_in2drkno_slc3v3( \ 29312 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt, \ 29313 di1=d1, di2=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0], \ 29314 maxngridsin=sin.shape[1]) 29315 elif Nvdims == 1 and len(sN.shape) == 3: 29316 voutt=fsci.module_scientific.multi_spaceweightstats_in1drkno_slc3v3( \ 29317 varin=varvt, idv=ind, ngridsin=sNt, gridsin=sint, percentages=spt, \ 29318 di1=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0], maxngridsin=sin.shape[1]) 29234 29319 else: 29235 29320 print errormsg … … 29237 29322 ' slice rank:', len(sN.shape), 'not ready!!' 29238 29323 print ' available ones _______' 29324 print ' var_rank: 3 slice_rank: 4' 29239 29325 print ' var_rank: 3 slice_rank: 3' 29240 print ' var_rank: 3 slice_rank: 4' 29326 print ' var_rank: 2 slice_rank: 3' 29327 print ' var_rank: 1 slice_rank: 3' 29241 29328 quit(-1) 29242 29329 29243 29330 vout = voutt.transpose() 29244 29331 ist = 0 29332 print ' Lluis filling file with statistics ...' 29245 29333 for statn in statns: 29246 29334 newvar = onewvars[statn] 29247 29335 rightvals = np.where(vout[6] == 0, gen.fillValueF, vout[ist]) 29248 newvar[tuple( slcv)] = rightvals29336 newvar[tuple(nslcv)] = rightvals 29249 29337 onewnc.sync() 29250 29338 ist = ist + 1
Note: See TracChangeset
for help on using the changeset viewer.
