![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
Changeset 2303 in lmdz_wrf for trunk/tools
- Timestamp:
- Jan 30, 2019, 2:56:09 PM (6 years ago)
- File:
-
- 1 edited
-
trunk/tools/nc_var_tools.py (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/nc_var_tools.py
r2300 r2303 28335 28335 28336 28336 if Nt > 0: 28337 print 'Lluis:', islc, ':', Nt 28338 gridsin[0,0:Nt,0,islc] = ainds[0:Nt,0] 28339 gridsin[1,0:Nt,0,islc] = ainds[0:Nt,1] 28337 gridsin[1,0:Nt,0,islc] = ainds[0:Nt,0] 28338 gridsin[0,0:Nt,0,islc] = ainds[0:Nt,1] 28340 28339 percens[0:Nt,0,islc] = 1. 28341 28340 for iv in range(Nt): … … 28343 28342 j = ainds[iv,1] 28344 28343 areas[j,i] = 1. 28345 if j == 94: print islc,'Lluis:', j,i,areas[j,i]28346 28344 28347 28345 # Values for file … … 28427 28425 anewvar.setncattr('coordinates',' '.join(varfinaldims[::-1])) 28428 28426 anewvar[:] = areas[:] 28429 if dn == 'HGT':28430 for i in range(areas.shape[1]):28431 if areas[94,i] == 1.: print 'Lluis 2:', 94,i,areas[94,i]28432 28427 28433 28428 aanewvar = onewnc.createVariable(dn+'area','f',tuple(Srgrid), \ … … 28533 28528 sliceinBt = sliceinB.transpose() 28534 28529 28535 print 'Lluis shapes sliceNAt:', sliceNAt.shape, 'sliceinAt:', sliceinAt.shape28536 print 'Lluis shapes sliceNBt:', sliceNBt.shape, 'sliceinBt:', sliceinBt.shape28537 28538 28530 NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \ 28539 28531 npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt, \ … … 28545 28537 28546 28538 # Remembering that it is python (C-like...) 28547 # inpointsA = inpointsA 28548 # inpointsB = inpointsB 28549 28550 jB = 0 28551 iB = 3 28552 jA = 10 28553 iA = 2 28554 print 'iA: _____ ' 28555 for iv in range(sliceNA[jA,iA]): 28556 print iv, sliceinA[:,iv,jA,iA] 28557 print 'iB: _____ ' 28558 for iv in range(sliceNB[jB,iB]): 28559 print iv, sliceinB[:,iv,jB,iB] 28560 for iv in range(sliceNA[jA,iA]): 28561 for iv2 in range(sliceNB[jB,iB]): 28562 if sliceinA[0,iv,jA,iA] == sliceinB[0,iv2,jB,iB] and sliceinA[1,iv,jA,iA] == sliceinB[1,iv2,jB,iB]: 28563 print iv,':', sliceinA[:,iv,jA,iA], ';', iv2, ':', sliceinB[:,iv2,jB,iB] 28564 28565 print 'Lluis jB iB jA iA _______',jB,iB,jA,iA,':',NpointsAB[jB,iB,jA,iA] 28566 print 'Lluis inpA shapes', inpointsA.shape, ':', inpointsA[:,jA,iA] 28567 print 'Lluis inpB shapes', inpointsB.shape, ':', inpointsB[:,jB,iB] 28568 for iv in range(NpointsAB[jB,iB,jA,iA]): 28569 print iv, 'NA:', sliceNA[jA,iA], 'pA:', inpointsA[iv,jA,iA], 'c:', sliceinA[:,inpointsA[iv,jA,iA],jA,iA], \ 28570 'NB:', sliceNB[jB,iB], 'pB:', inpointsB[iv,jB,iB], 'c:', sliceinB[:,inpointsB[iv,jB,iB],jB,iB], '=', \ 28571 pointsAB[:,iv,jB,iB,jA,iA] 28572 28573 quit() 28539 inpointsA = inpointsA-1 28540 inpointsB = inpointsB-1 28574 28541 28575 28542 maxNpointsAB = np.max(NpointsAB) … … 28615 28582 for iB in range(dxB): 28616 28583 Nin = NpointsAB[jB,iB,jA,iA] 28617 print 'jB iB jA iA:', jB,iB,jA,iA, ':', Nin28618 28584 innewvar[:,0:Nin,jB,iB,jA,iA] = pointsAB[:,0:Nin,jB,iB,jA,iA] 28585 aanewvar[jB,iB,jA,iA]= gen.fillValueF 28586 anewvar[:,jB,iB,jA,iA]= gen.fillValueF 28587 pnewvar[:,jB,iB,jA,iA]= gen.fillValueF 28619 28588 slicearea = 0. 28620 28589 for iv in range(Nin): 28621 pA = oslicepA[inpointsA[iv,j A,iA],jA,iA]28622 pB = oslicepB[inpointsB[iv,jB,iB ],jB,iB]28590 pA = oslicepA[inpointsA[iv,jB,iB,jA,iA],jA,iA] 28591 pB = oslicepB[inpointsB[iv,jB,iB,jA,iA],jB,iB] 28623 28592 pnewvar[iv,jB,iB,jA,iA]= pA*pB 28624 ixA = sliceinA[1,inpointsA[iv,j A,iA],jA,iA]28625 iyA = sliceinA[0,inpointsA[iv,j A,iA],jA,iA]28593 ixA = sliceinA[1,inpointsA[iv,jB,iB,jA,iA],jA,iA] 28594 iyA = sliceinA[0,inpointsA[iv,jB,iB,jA,iA],jA,iA] 28626 28595 aA = osliceaA[iyA,ixA] 28627 ixB = sliceinB[1,inpointsB[iv,jB,iB ],jB,iB]28628 iyB = sliceinB[0,inpointsB[iv,jB,iB ],jB,iB]28596 ixB = sliceinB[1,inpointsB[iv,jB,iB,jA,iA],jB,iB] 28597 iyB = sliceinB[0,inpointsB[iv,jB,iB,jA,iA],jB,iB] 28629 28598 aB = osliceaB[iyB,ixB] 28630 # I do not understand why this is needed !! 28631 if type(aA)!=type(gen.mamat[1]) and \ 28632 type(aB)!=type(gen.mamat[1]): 28633 anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB) 28634 slicearea = slicearea + (aA*pA)*(aB*pB) 28635 else: 28636 print type(aA) != type(gen.mamat[1]), type(aB) != type(gen.mamat[1]) 28637 print 'Lluis aA ________', iyA, ixA, 'shape osliceaA:', osliceaA.shape, ':', aA, gen.mamat[1] 28638 print osliceaA[iyA-2:iyA+3,ixA-2:ixA+3] 28639 print 'Lluis aB ________', iyB, ixB, 'shape osliceaB:', osliceaB.shape, ':', aB 28640 print osliceaB[iyB-2:iyB+3,ixB-2:ixB+3] 28641 #quit() 28642 aanewvar[jB,iB,jA,iA]= slicearea 28599 anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB) 28600 slicearea = slicearea + (aA*pA)*(aB*pB) 28601 aanewvar[jB,iB,jA,iA]= slicearea 28643 28602 28644 28603 onewnc.sync() … … 28734 28693 Nin = NpointsAB[jB,iB,jA,iA] 28735 28694 innewvar[:,0:Nin,jB,iB,jA,iA] = pointsAB[:,0:Nin,jB,iB,jA,iA] 28695 aanewvar[jB,iB,jA,iA]= gen.fillValueF 28696 anewvar[:,jB,iB,jA,iA]= gen.fillValueF 28697 pnewvar[:,jB,iB,jA,iA]= gen.fillValueF 28736 28698 slicearea = 0. 28737 28699 for iv in range(Nin): 28738 pA = oslicepA[inpointsA[iv,j A,iA],jA,iA]28739 pB = oslicepB[inpointsB[iv,jB,iB ],jB,iB]28700 pA = oslicepA[inpointsA[iv,jB,iB,jA,iA],jA,iA] 28701 pB = oslicepB[inpointsB[iv,jB,iB,jA,iA],jB,iB] 28740 28702 pnewvar[iv,jB,iB,jA,iA]= pA*pB 28741 ixA = sliceinA[1,inpointsA[iv,j A,iA],jA,iA]28742 iyA = sliceinA[0,inpointsA[iv,j A,iA],jA,iA]28703 ixA = sliceinA[1,inpointsA[iv,jB,iB,jA,iA],jA,iA] 28704 iyA = sliceinA[0,inpointsA[iv,jB,iB,jA,iA],jA,iA] 28743 28705 aA = osliceaA[iyA,ixA] 28744 ixB = sliceinB[1,inpointsB[iv,jB,iB ],jB,iB]28745 iyB = sliceinB[0,inpointsB[iv,jB,iB ],jB,iB]28706 ixB = sliceinB[1,inpointsB[iv,jB,iB,jA,iA],jB,iB] 28707 iyB = sliceinB[0,inpointsB[iv,jB,iB,jA,iA],jB,iB] 28746 28708 aB = osliceaB[iyB,ixB] 28747 # I do not understand why this is needed !! 28748 if type(aA)!=type(gen.mamat[1]) and \ 28749 type(aB)!=type(gen.mamat[1]): 28750 slicearea = slicearea + (aA*pA)*(aB*pB) 28751 aanewvar[jB,iB,jA,iA]= slicearea 28752 28709 anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB) 28710 slicearea = slicearea + (aA*pA)*(aB*pB) 28711 aanewvar[jB,iB,jA,iA]= slicearea 28753 28712 onewnc.sync() 28754 28713 … … 28768 28727 sga = oslicegridarea[:] 28769 28728 28729 # Percentages as the areal fraction of each grid with respect slice size 28730 sgpa = np.ones(tuple(oslicegridpercen.shape), dtype=np.float)*gen.fillValueF 28731 panewvar = onewnc.createVariable(dn+'areagridpercen','f',tuple(Spgrid), \ 28732 fill_value=gen.fillValueF) 28733 basicvardef(panewvar ,dn+'areagridpercen', "percentages of the " + \ 28734 "grids cells from " + newslcvarns[0] + " laying within " + newslcvarns[1] + \ 28735 'as grid-sliced area/slice-area', '1') 28736 pnewvar.setncattr('coordinates',' '.join(Spgrid[::-1])) 28737 28738 (dyA, dxA, dyB, dxB) = sN.shape 28739 for jA in range(dyA): 28740 for iA in range(dxA): 28741 for jB in range(dyB): 28742 for iB in range(dxB): 28743 Nin = sN[jA,iA,jB,iB] 28744 sgpa[0:Nin,jA,iA,jB,iB] = anewvar[0:Nin,jA,iA,jB,iB]/ \ 28745 aanewvar[jA,iA,jB,iB] 28746 panewvar[0:Nin,jA,iA,jB,iB] = sgpa[0:Nin,jA,iA,jB,iB] 28747 onewnc.sync() 28748 28770 28749 # removing monotones. Fortran only accepts rank <= 7 arrays !! 28771 28750 sN, Ndims = gen.remove_monotones(sN, oNslice.dimensions) … … 28773 28752 sa, adims = gen.remove_monotones(sa, oslicearea.dimensions) 28774 28753 sga, gadims = gen.remove_monotones(sga, oslicegridarea.dimensions) 28775 28776 # Percentages as the areal fraction of each grid with respect slice size 28777 sgp = sga/sa 28778 i=3 28779 j=10 28780 k=2 28781 print 'Lluis percentages ________', i, j,k 28782 Npt = sN[i,j,k] 28783 for iv in range(Npt): 28784 print 'Lluis:', iv,'pt:', sin[:,iv,i,j,k], 'a:', sa[i,j,k], 'ga:', sga[iv,i,j,k] 28785 quit() 28786 28754 sgpa, gpadims = gen.remove_monotones(sgpa, oslicegridpercen.dimensions) 28755 28787 28756 ssh = sN.shape 28788 28757 dimsslice = list(Ndims) … … 28791 28760 sNt = sN.transpose() 28792 28761 sint = sin.transpose() 28793 spt = sgp .transpose()28762 spt = sgpa.transpose() 28794 28763 28795 28764 for varn in varns: … … 28865 28834 varvt = varv.transpose() 28866 28835 28867 print 'Lluis vshape:', vshape, 'sN:', sN.shape28868 28869 28836 if len(vshape) == 3 and len(sN.shape) == 4: 28870 28837 voutt = fsci.module_scientific.multi_spacewightstats_in3drk3_slc3v4( \ … … 28888 28855 quit(-1) 28889 28856 28890 print 'Lluis shapes vout:', vout.shape28891 28892 28857 ist = 0 28893 28858 for statn in statns: … … 28899 28864 ist = ist + 1 28900 28865 28901 quit()28902 28903 # mask in var slice28904 maskvarslcs = {}28905 allNslices = []28906 slicedimn = []28907 for vslcn in slcvarns:28908 varslcv = slicesinf[vslcn]28909 Nslices = varslcv[0]28910 ovardims = varslcv[1]28911 slcvar = varslcv[2]28912 28913 NOTcoinc = list(set(ovardims) - set(vdimns))28914 if len(NOTcoinc) > 0:28915 notcoincdim = {}28916 for dn in NOTcoinc: notcoincdim[dn] = 028917 else:28918 notcoincdim = None28919 maskvarslice = []28920 for islc in range(Nslices):28921 varmask = gen.mat_dimreshape(list(ovardims), slcvar[islc,], vdimns, \28922 vshape, notcoincdim)28923 maskvarslice.append(varmask)28924 28925 maskvarslcs[vslcn] = maskvarslice28926 allNslices.append(Nslices)28927 slicedimn.append('slice_'+vslcn)28928 28929 newvarshape = allNslices + vshape28930 print ' Shape final sliced variable:', newvarshape28931 28932 # variables28933 28934 # This can get quite huge !!28935 #newvard = slicedimn + vdimns28936 #newvar = onewnc.createVariable(varn+'sliced', 'f', tuple(newvard), \28937 # fill_value=gen.fillValueF)28938 #basicvardef(newvar, varn+'sliced', varn + 'sliced by '+ slcvarnS, vu)28939 #add_varattrs(onc, onewnc, [varn], [varn+'sliced'])28940 28941 masknodims = []28942 vrank = np.arange(len(vdimns))28943 if slicestatsdim == 'any':28944 slcmin = np.ones(tuple(allNslices), dtype=ovar.dtype)*gen.fillValueF28945 slcmax = np.ones(tuple(allNslices), dtype=ovar.dtype)*gen.fillValueF28946 slcmean = np.ones(tuple(allNslices), dtype=ovar.dtype)*gen.fillValueF28947 slcstd= np.ones(tuple(allNslices), dtype=ovar.dtype)*gen.fillValueF28948 newvard = slicedimn28949 else:28950 nodims = gen.str_list(slicestatsdim, ',')28951 stvardims = vdimns + []28952 stvdn = []28953 stvid = []28954 idd = 028955 vrank = list(vrank)28956 for nod in nodims:28957 if not gen.searchInlist(onc.dimensions, nod):28958 print errormsg28959 print ' ' +fname+ ": file '" + ncfile + "' does not have " + \28960 "dimension: '" + dn + "' !!"28961 dimns = list(onc.dimensions)28962 dimns.sort()28963 print ' available ones:', dimns28964 quit(-1)28965 28966 if gen.searchInlist(vdimns, nod):28967 stvdn.append(nod)28968 stvid.append(vshape[idd])28969 stvardims.remove(nod)28970 vrank.pop(idd)28971 masknodims.append(slice(0,vshape[idd]))28972 idd = idd + 128973 vrank = np.array(vrank)28974 28975 slcmin = np.ones(tuple(allNslices+stvid), dtype=ovar.dtype)*gen.fillValueF28976 slcmax = np.ones(tuple(allNslices+stvid), dtype=ovar.dtype)*gen.fillValueF28977 slcmean = np.ones(tuple(allNslices+stvid), dtype=ovar.dtype)*gen.fillValueF28978 slcstd= np.ones(tuple(allNslices+stvid), dtype=ovar.dtype)*gen.fillValueF28979 newvard = slicedimn + stvdn28980 28981 newvarN = onewnc.createVariable(varn+'Nsliced', 'f', tuple(newvard), \28982 fill_value=gen.fillValueF)28983 basicvardef(newvarN, varn+'Nsliced','number of values in slice of ' + varn + \28984 'sliced by '+ slcvarnS, vu)28985 add_varattrs(onc, onewnc, [varn], [varn+'Nsliced'])28986 28987 newvarn = onewnc.createVariable(varn+'minsliced', 'f', tuple(newvard), \28988 fill_value=gen.fillValueF)28989 basicvardef(newvarn, varn+'minsliced', 'minimum value of '+varn+'sliced by '+\28990 slcvarnS, vu)28991 add_varattrs(onc, onewnc, [varn], [varn+'minsliced'])28992 28993 newvarx = onewnc.createVariable(varn+'maxsliced', 'f', tuple(newvard), \28994 fill_value=gen.fillValueF)28995 basicvardef(newvarx, varn+'maxsliced', 'maximum value of '+varn+'sliced by '+\28996 slcvarnS, vu)28997 add_varattrs(onc, onewnc, [varn], [varn+'maxsliced'])28998 28999 newvarm = onewnc.createVariable(varn+'meansliced', 'f', tuple(newvard), \29000 fill_value=gen.fillValueF)29001 basicvardef(newvarm, varn+'meansliced', 'mean value of '+varn+'sliced by ' + \29002 slcvarnS, vu)29003 add_varattrs(onc, onewnc, [varn], [varn+'meansliced'])29004 29005 newvars = onewnc.createVariable(varn+'stdsliced', 'f', tuple(newvard), \29006 fill_value=gen.fillValueF)29007 basicvardef(newvars, varn+'stdsliced', 'standard deviation of '+varn+ \29008 'sliced by ' + slcvarnS, vu)29009 add_varattrs(onc, onewnc, [varn], [varn+'stdsliced'])29010 29011 # Obviously we get memory problems... proceed slide by slide to fill the file29012 #newvarNmasked = np.ones(tuple(newvarshape), dtype=ovar.dtype)*gen.fillValueF29013 newvarNmasked = np.ones(tuple(vshape), dtype=ovar.dtype)*gen.fillValueF29014 29015 # Masking!29016 allslices = gen.all_consecutive_combs(allNslices)29017 Nallslices = len(allslices)29018 29019 for iallslc in range(Nallslices):29020 islcN = allslices[iallslc]29021 slcnewvar = []29022 iv = 029023 for vslcn in slcvarns:29024 masksv = maskvarslcs[vslcn]29025 slcnewvar.append(islcN[iv])29026 if iv == 0: newmask = masksv[islcN[iv]]29027 else: newmask = newmask+masksv[islcN[iv]]29028 iv = iv + 129029 29030 for iid in vshape:29031 slcnewvar.append(slice(0,iid))29032 29033 mavals = ma.array(varv, mask=newmask)29034 #newvarNmasked[tuple(slcnewvar)] = mavals.filled(gen.fillValueF)29035 newvarNmasked[:] = mavals.filled(gen.fillValueF)29036 29037 # This can get quite huge !!!29038 #newvar[tuple(slcnewvar)] = newvarNmasked[:]29039 29040 manewvar = ma.masked_equal(newvarNmasked, gen.fillValueF)29041 29042 if len(masknodims) > 0:29043 islcst = list(islcN) + list(masknodims)29044 if iallslc == 0:29045 print ' variable slice statisitcs not using:', nodims29046 print ' size statistics variables:', allNslices+stvid29047 print ' slice statistics variables:', islcst29048 print ' running shapes of the variable:', vrank29049 else: islcst = islcN29050 29051 newvarN[tuple(islcst)]=np.sum(~manewvar.mask, axis=tuple(vrank))29052 29053 slcmin[tuple(islcst)] = np.min(manewvar, axis=tuple(vrank))29054 slcmax[tuple(islcst)] = np.max(manewvar, axis=tuple(vrank))29055 slcmean[tuple(islcst)] = np.mean(manewvar, axis=tuple(vrank))29056 slcstd[tuple(islcst)] = np.std(manewvar, axis=tuple(vrank))29057 29058 newvarn[:] = slcmin29059 newvarx[:] = slcmax29060 newvarm[:] = slcmean29061 newvars[:] = slcstd29062 29063 28866 onewnc.sync() 29064 # Adding dimension-variables29065 for dn in vdimns:29066 if not dimvars.has_key(dn):29067 print errormsg29068 print ' ' + fname + ": no dimension-variable provided for " + \29069 "dimension '" + dn + "' !!"29070 print ' provided ones _______'29071 gen.printing_dictionary(dimvars)29072 quit(-1)29073 if not onewnc.variables.has_key(dimvars[dn]):29074 add_vars(onc, onewnc, [dimvars[dn]])29075 onewnc.sync()29076 28867 29077 28868 # Add global attributes
Note: See TracChangeset
for help on using the changeset viewer.
