Changeset 706 in lmdz_wrf for trunk

Timestamp:

Apr 6, 2016, 4:34:35 PM (9 years ago)

Author:

lfita

Message:

Adding CoarselonlatFind' lonlatFind' functions to find values on lon, lat 2D matrices

• `CoarselonlatFind?': uses a first guess based on 'coarser' (sub sample) search from the original lon, lat 2D matrices (speed up, but only for continuos projections! :( )

Location:

trunk/tools

Files:

• nc_var.py (modified) (1 diff)
• nc_var_tools.py (modified) (18 diffs)

trunk/tools/nc_var.py

@@ -10 +10 @@
 ## e.g. # nc_var.py -o filter_2dim -S '80,y,x,lon,lat' -f 'tahiti_5m_ll.grd' -v 'z'
 ## e.g. # nc_var.py -o selvar -f /home/lluis/PY/met_em.d01.1979-01-01_00:00:00.nc -S 'west_east@XLONG_M,south_north@XLAT_M,num_metgrid_levels@int,Time@Times' -v TT,UU,VV,SKINTEMP
+## e.g. # nc_var.py -o 'Partialmap_Entiremap' -f carteveg5km.nc -S 'longitude,latitude,std,5000.,Goode,Goode_5km.nc' -v vegetation_map
 
 from optparse import OptionParser

trunk/tools/nc_var_tools.py

@@ -3 +3 @@
 import os
 import re
+import numpy.ma as ma
 
 main = 'nc_var_tools.py'
@@ -242 +243 @@
 
     return newstring
+
+def period_information(idate, edate, totunits):
+    """ Function to provide the information of a given period idate, edate (dates in [YYYY][MM][DD][HH][MI][SS] format)
+      [idate]= initial date of the period (in [YYYY][MM][DD][HH][MI][SS] format)
+      [edate]= end date of the period (in [YYYY][MM][DD][HH][MI][SS] format)
+      [totunits]= latest number of entire temporal units to compute: 'year', 'month', 'day', 'hour', 'minute', 'second'
+
+    resultant information given as ',' list of:
+      [dT]: sign of the temporal increment ('pos', 'neg')
+      [idate]: initial date as [YYYY]:[MM]:[DD]:[HH]:[MI]:[SS]
+      [edate]: end date as [YYYY]:[MM]:[DD]:[HH]:[MI]:[SS]
+      [Nsecs]: total seconds between dates
+      [Nyears]: Number of years taken as entire units ([iYYYY]0101000000 to [eYYYY+1]0101000000)
+      [ExactYears]: Exact year-dates of the period as a '@' separated list of  [YYYY]0101000000
+      [Nmonths]: Number of months taken as entire units ([iYYYY][iMM]01000000 to [eYYYY][eMM+1]01000000)
+      [ExactMonths]: Exact mon-dates of the period as  a '@' separated list of [YYYY][MM]01000000
+      [Ndays]: Number of days taken as entire units ([iYYYY][iMM][iDD]000000 to [eYYYY][eMM][eDD+1]000000)
+      [ExactDays]: Exact day-dates of the period as  a '@' separated list of [YYYY][MM][DD]000000
+      [Nhours]: Number of hours taken as entire units ([iYYYY][iMM][iDD][iHH]0000 to [eYYYY][eMM][eDD][eHH+1]0000)
+      [ExactHours]: Exact hour-dates of the period as a '@' separated list of [YYYY][MM][DD][HH]0000
+      [Nminutes]: Number of minutes taken as entire units ([iYYYY][iMM][iDD][iHH][iMI]00 to [eYYYY][eMM][eDD][eHH][eMI+1]00)
+      [ExactMinutes]: Exact minutes-dates of the period as a '@' separated list of [YYYY][MM][DD][HH][MI]00
+      [Nsecs]: Number of minutes taken as entire units ([iYYYY][iMM][iDD][iHH][iMI][iS] to [eYYYY][eMM][eDD][eHH][eMI][eSE+1])
+      [ExactSeconds]: Exact seconds-dates of the period as a '@' separated list of [YYYY][MM][DD][HH][MI][SS]
+
+    >>> period_information( '19760217083110', '19770828000000', 'month')
+    pos,1976:02:17:08:31:10,1977:08:28:00:00:00,48180530.0,19,19760201000000@19760301000000@19760401000000@
+    19760501000000@19760601000000@19760701000000@19760801000000@19760901000000@19761001000000@19761101000000@
+    19761201000000@19770101000000@19770201000000@19770301000000@19770401000000@19770501000000@19770601000000@
+    19770701000000@19770801000000@19770901000000
+    """
+    import datetime as dt
+
+    fname = 'period_information'
+
+    if idate == 'h':
+        print fname + '_____________________________________________________________'
+        print variables_values.__doc__
+        quit()
+
+    expectargs = '[idate],[edate],[totunits]'
+ 
+    check_arguments(fname,str(idate)+','+str(edate)+','+str(totunits),expectargs,',')
+
+# Checking length
+    Lidate = len(idate)
+    Ledate = len(edate)
+
+    if Lidate != Ledate:
+        print errormsg
+        print '   ' + fname + ': string dates of the period have different length !!'
+        print "    idate: '" + idate + "' (L=", Lidate, ") edate: '" + edate +       \
+          "' (L=", Ledate, ')'
+        quit(-1)
+    else:
+        if Lidate != 14:
+            print errormsg
+            print '  ' + fname + ": wrong initial date: '" + idate + "' must " +     \
+              'have 14 characters!!'
+            print '    and it has:', Lidate
+            quit(-1)
+        if Ledate != 14:
+            print errormsg
+            print '  ' + fname + ": wrong end date: '" + edate + "' must " +         \
+              'have 14 characters!!'
+            print '    and it has:', Ledate
+            quit(-1)
+
+# Checking for right order of dates
+    if int(idate) > int(edate):
+        print warnmsg
+        print '  ' + fname + ": initial date '" + idate + "' after end date '" +     \
+          edate + "' !!"
+        print "  re-sorting them!"
+        i1 = edate
+        e1 = idate
+        idate = i1 
+        edate = e1
+        oinf = 'neg'
+    else:
+        oinf = 'pos'
+
+# Year, month, day, hour, minute, second initial date
+    iyrS = idate[0:4]
+    imoS = idate[4:6]
+    idaS = idate[6:8]
+    ihoS = idate[8:10]
+    imiS = idate[10:12]
+    iseS = idate[12:14]
+
+    oinf = oinf + ',' + iyrS+':'+imoS+':'+idaS+':'+ihoS+':'+imiS+':'+iseS
+
+# Year, month, day, hour, minute, second end date
+    eyrS = edate[0:4]
+    emoS = edate[4:6]
+    edaS = edate[6:8]
+    ehoS = edate[8:10]
+    emiS = edate[10:12]
+    eseS = edate[12:14]
+
+    oinf = oinf + ',' + eyrS+':'+emoS+':'+edaS+':'+ehoS+':'+emiS+':'+eseS
+
+    idateT = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),int(iseS))
+    edateT = dt.datetime(int(eyrS),int(emoS),int(edaS),int(ehoS),int(emiS),int(eseS))
+
+# Seconds between dates
+    DT = edateT - idateT
+    Nsecs = DT.total_seconds()
+
+    oinf = oinf + ',' + str(Nsecs)
+
+# Looking for number of years tacking exact beginning of the units [iYYYY]0101000000, [eYYYY+1]0101000000
+##
+    if totunits == 'year':
+        itime = int(iyrS)
+        if edate[4:15] != '0101000000':
+            etime = int(eyrS)+1
+        else:
+            etime = int(eyrS)
+
+        itimeT = dt.datetime(itime,1,1,0,0,0)
+
+        Nyears = 1
+        ExactYears = itimeT.strftime("%Y%m%d%H%M%S")
+        it = int(iyrS)
+        while it + 1 <= etime:
+            it = it + 1
+            Nyears = Nyears + 1
+            itT = dt.datetime(it,1,1,0,0,0)
+            ExactYears = ExactYears + '@' + itT.strftime("%Y%m%d%H%M%S")
+
+        oinf = oinf + ',' + str(Nyears) + ',' + ExactYears
+
+# Looking for number of months tacking exact beginning of the units [iYYYY][iMM]01000000, [eYYYY][eMM+1]01000000
+##
+    if totunits == 'month':
+        itime = int(iyrS)*100 + int(imoS)
+        if edate[6:15] != '01000000':
+            emo1 = int(emoS)+1
+            if emo1 > 13:
+                eyr1 = str(int(eyrS) + 1)
+                emo1 = 01
+            else:
+                eyr1 = eyrS
+            
+        else:
+            eyr1 = eyrS
+            emo1 = emoS
+    
+        etime = int(eyr1)*100 + int(emo1) 
+        itimeT = dt.datetime(int(iyrS),int(imoS),1,0,0,0)
+        
+        Nmonths = 1
+        ExactMonths = itimeT.strftime("%Y%m%d%H%M%S")
+        it = itime
+        while it + 1 <= etime:
+            it = it + 1
+            Nmonths = Nmonths + 1
+            ityr = it/100
+            itmo = it - ityr*100
+            if itmo > 12:
+                ityr = ityr + 1
+                itmo = 1
+                it = ityr * 100 + itmo
+    
+            itT = dt.datetime(ityr,itmo,1,0,0,0)
+            ExactMonths = ExactMonths + '@' + itT.strftime("%Y%m%d%H%M%S")
+    
+        oinf = oinf + ',' + str(Nmonths) + ',' + ExactMonths
+
+# Looking for number of days tacking exact beginning of the units [iYYYY][iMM][iDD]000000, [eYYYY][eMM][eDD+1]000000
+##
+    if totunits == 'day':
+        itime = dt.datetime(int(iyrS),int(imoS),int(idaS),0,0,0)
+        if edate[8:15] != '000000':
+            etime = dt.datetime(int(eyrS), int(emoS), int(edaS)+1,0,0,0)
+        else:
+            etime = edateT
+    
+        Ndays = 1
+        ExactDays = itime.strftime("%Y%m%d%H%M%S")
+        it = itime
+        while it + dt.timedelta(days=1) <= etime:
+            it = it + dt.timedelta(days=1)
+            Ndays = Ndays + 1
+            ExactDays = ExactDays + '@' + it.strftime("%Y%m%d%H%M%S")
+    
+        oinf = oinf + ',' + str(Ndays) + ',' + ExactDays
+
+# Looking for number of hours tacking exact beginning of the units [iYYYY][iMM][iDD][iHH]0000, [eYYYY][eMM][eDD][iHH+1]0000
+##
+    if totunits == 'hour':
+        itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),0,0)
+        if edate[10:15] != '0000':
+            etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS)+1,0,0)
+        else:
+            etime = edateT
+    
+        Nhours = 1
+        ExactHours = itime.strftime("%Y%m%d%H%M%S")
+        it = itime
+        while it  + dt.timedelta(hours=1) <= etime:
+            it = it + dt.timedelta(hours=1)
+            Nhours = Nhours + 1
+            ExactHours = ExactHours + '@' + it.strftime("%Y%m%d%H%M%S")
+    
+        oinf = oinf + ',' + str(Nhours) + ',' + ExactHours
+
+# Looking for number of minutes tacking exact beginning of the units [iYYYY][iMM][iDD][iHH][iMI]00, [eYYYY][eMM][eDD][iHH][iMI+1]00
+##
+    if totunits == 'minute':
+        itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),0)
+        if edate[12:15] != '00':
+            etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS), int(emiS)+1,0)
+        else:
+            etime = edateT
+
+        Nminutes = 1
+        ExactMinutes = itime.strftime("%Y%m%d%H%M%S")
+        it = itime
+        while it + dt.timedelta(minutes=1) <= etime:
+            it = it + dt.timedelta(minutes=1)
+            Nminutes = Nminutes + 1
+            ExactMinutes = ExactMinutes + '@' + it.strftime("%Y%m%d%H%M%S")
+    
+        oinf = oinf + ',' + str(Nminutes) + ',' + ExactMinutes
+
+# Looking for number of seconds tacking exact beginning of the units [iYYYY][iMM][iDD][iHH][iMI][iSE], 
+#   [eYYYY][eMM][eDD][iHH][iMI][iSE+1]
+##
+    if totunits == 'second':
+        itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),int(iseS))
+        print 'Lluis ', edate[12:15], '00'
+        if edate[12:15] != '00':
+            etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS), int(emiS), int(eseS)+1)
+        else:
+            etime = edateT
+
+        Nseconds = 1
+        ExactSeconds = itime.strftime("%Y%m%d%H%M%S")
+        it = itime
+        while it + dt.timedelta(seconds=1) < etime:
+            it = it + dt.timedelta(seconds=1)
+            Nseconds = Nseconds + 1
+            ExactSeconds = ExactSeconds + '@' + it.strftime("%Y%m%d%H%M%S")
+    
+        oinf = oinf + ',' + str(Nseconds) + ',' + ExactSeconds
+
+    return oinf
+
+# LluisWORKING
 
 ####### ###### ##### #### ### ## #
@@ -6018 +6270 @@
                 quit(-1)
 
-
 class cls_seasonal_means(object):
     """ Class to compute the seasonal mean of a time series of values
@@ -12727 +12978 @@
 
     return origvar, dimsnewvar
+
+def var_3desc(ovar):
+    """ Function to provide std_name, long_name and units from an object variable
+      ovar= object variable
+    """
+    fname = 'var_desc'
+    varattrs = ovar.ncattrs()
+
+    if searchInlist(varattrs,'std_name'):
+        stdn = ovar.getncattr('std_name')
+    else:
+        vvalues = variables_values(ovar._name)
+        stdn = vvalues[1]
+
+    if searchInlist(varattrs,'long_name'):
+        lonn = ovar.getncattr('long_name')
+    else:
+        lonn = vvalues[4].replace('|',' ')
+
+    if searchInlist(varattrs,'units'):
+        un = ovar.getncattr('units')
+    else:
+        un = vvalues[5]
+
+    return stdn, lonn, un
 
 def file_oper_alongdims(values, ncfile, varn):
@@ -17773 +18049 @@
 #SpatialWeightedMean('reglonlat,XLONG,XLAT,west_east,south_north','wrfout_d01_2001-11-11_00:00:00','HGT')
 
-def fill_value(vartype, fillval0):
+def fillvalue_kind(vartype, fillval0):
     """ Function to provide the fiven fill_Value for a kind of variable
       [vartype]= type of the variable
@@ -17784 +18060 @@
            Integer32 = -99999
     """
-    fname = 'fill_value'
-
-    if vartype == '|S':
+    fname = 'fillvalue_kind'
+
+    if vartype == type('c'):
         if fillval0 == 'std':
             fillval = fillvalueC
@@ -17798 +18074 @@
     elif vartype == type(np.int16(1)):
         if fillval0 == 'std':
-            fillval = np.int16(fillValueI)
+            fillval = np.int(9999999)
         else:
             fillval = np.int16(fillval0)
@@ -17828 +18104 @@
     return fillval
 
+def nctype(vartype):
+    """ Function to provide the string for the variable creation in a netCDF file
+      [vartype]= type of the variable
+    """
+    fname = 'nctype'
+
+    if vartype == type('c'):
+         ncs = 'c'
+    elif vartype == type(int(1)):
+         ncs = 'i'
+    elif vartype == type(np.int16(1)):
+         ncs = 'i'
+    elif vartype == type(np.int32(1)):
+         ncs = 'i8'
+    elif vartype == type(np.float(1.)):
+         ncs = 'f'
+    elif vartype == type(np.float32(1.)):
+         ncs = 'f4'
+    elif vartype == type(np.float64(1.)):
+         ncs = 'f8'
+    else:
+        print errormsg
+        print '  ' + fname + ': variable type', vartype, 'not ready !!'
+        quit(-1)
+
+    return ncs
+
 def VarVal_FillValue(values, filen, varn):
     """ Function to transform a given value from a given variable to _FillValue in a netCDF file
@@ -18070 +18373 @@
 #quit()
 
+def CoarselonlatFind(ilon, ilat, lonv, latv, per):
+    """ Function to search a given value from a coarser version of the data
+      ilon, ilat= original 2D matrices with the longitudes and the latitudes
+      lonv, latv= longitude and latitude to find
+      per= period (as fraction over 1) of the fractions of the original grid to use to explore
+      >>> npt=100
+      >>> lonval0 = np.arange(0.,npt*1.,1.)
+      >>> latval0 = np.arange(0.,npt*1.,1.)
+      >>> lonval = np.zeros((npt,npt), dtype=np.float)
+      >>> latval = np.zeros((npt,npt), dtype=np.float)
+
+      >>> for i in range(npt):
+      >>>     lonval[:,i] = lonval0[i]
+      >>>     latval[i,:] = latval0[i]
+      >>> CoarselonlatFind(lonval, latval, 5.25, 5.25, .1)
+      (array([5, 5]), 0.35355339059327379)
+    """
+    fname = 'CoarselonlatFind'
+
+    dx = ilon.shape[1]
+    dy = ilon.shape[0]
+
+    nlon = np.min(ilon)
+    xlon = np.max(ilon)
+    nlat = np.min(ilat)
+    xlat = np.max(ilat)
+
+    if lonv < nlon or lonv > xlon:
+        print errormsg
+        print '  ' + fname + ': longitude outside data range!!'
+        print '    given value:', lonv,'outside [',nlon,',',xlon,']'
+        quit(-1) 
+    if latv < nlat or latv > xlat:
+        print errormsg
+        print '  ' + fname + ': latitude outside data range!!'
+        print '    given value:', latv,'outside [',nlat,',',xlat,']'
+        quit(-1) 
+
+    fracx = int(dx*per)
+    fracy = int(dy*per)
+
+    fraclon = ilon[::fracy,::fracx]
+    fraclat = ilat[::fracy,::fracx]
+    fracdx = fraclon.shape[1]
+    fracdy = fraclon.shape[0]
+
+#    print 'fraclon _______'
+#    print fraclon
+
+#    print 'fraclat _______'
+#    print fraclat
+
+# Fraction point
+    difffraclonlat = np.sqrt((fraclon-lonv)**2. + (fraclat-latv)**2.)
+    mindifffracLl = np.min(difffraclonlat)
+    ilatlonfrac = index_mat(difffraclonlat, mindifffracLl)
+
+#    print 'mindifffracLl:', mindifffracLl, 'ilatlonfrac:', ilatlonfrac
+#    print 'frac lon, lat:', fraclon[ilatlonfrac[0],ilatlonfrac[1]], ',',             \
+#      fraclat[ilatlonfrac[0],ilatlonfrac[1]]
+#    ilatlon = [ilatlonfrac[0]*fracy, ilatlonfrac[1]*fracx]
+#    print 'input index:',ilatlon, 'lon lat:', ilon[ilatlon[0],ilatlon[1]], ',',      \
+#      ilat[ilatlon[0],ilatlon[1]]
+      
+# Providing fraction range
+    if fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                            \
+      fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
+#        ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]-1]
+#        ifracend = [ilatlonfrac[0], ilatlonfrac[1]]
+        if ilatlonfrac[0] > 0: 
+            iybeg = (ilatlonfrac[0]-1)*fracy
+            iyend = ilatlonfrac[0]*fracy+1
+        else:
+            iybeg = 0
+            iyend = fracy+1
+        if ilatlonfrac[1] > 0: 
+            ixbeg = (ilatlonfrac[1]-1)*fracx
+            ixend = ilatlonfrac[1]*fracx+1
+        else:
+            ixbeg = 0
+            ixend = fracx+1
+    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                           \
+      fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
+#        ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]]
+#        ifracend = [ilatlonfrac[0], ilatlonfrac[1]+1]
+        if ilatlonfrac[0] > 0: 
+            iybeg = (ilatlonfrac[0]-1)*fracy
+            iyend = ilatlonfrac[0]*fracy+1
+        else:
+            iybeg = 0
+            iyend = fracy+1
+        if ilatlonfrac[1] < fracdx: 
+            ixbeg = (ilatlonfrac[1]-1)*fracx
+            ixend = ilatlonfrac[1]*fracx+1
+        else:
+            ixbeg = fracdx*fracx
+            ixend = dx+1
+    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                           \
+      fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
+#        ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]]
+#        ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]+1]
+        if ilatlonfrac[0] < fracdy: 
+            iybeg = (ilatlonfrac[0]-1)*fracy
+            iyend = ilatlonfrac[0]*fracy+1
+        else:
+            iybeg = fracdy*fracy
+            iyend = dy+1
+        if ilatlonfrac[1] < fracdx: 
+            ixbeg = (ilatlonfrac[1]-1)*fracx
+            ixend = ilatlonfrac[1]*fracx+1
+        else:
+            ixbeg = fracdx*fracx
+            ixend = dx+1
+    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                          \
+      fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
+#        ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]-1]
+#        ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]]
+        if ilatlonfrac[0] < fracdy: 
+            iybeg = (ilatlonfrac[0]-1)*fracy
+            iyend = ilatlonfrac[0]*fracy+1
+        else:
+            iybeg = fracdy*fracy
+            iyend = dy+1
+        if ilatlonfrac[1] > 0: 
+            ixbeg = (ilatlonfrac[1]-1)*fracx
+            ixend = ilatlonfrac[1]*fracx+1
+        else:
+            ixbeg = 0
+            ixend = fracx+1
+
+#    ibeg = [iybeg, ixbeg]
+#    iend = [iyend, ixend]
+#    print 'find window; beginning', ibeg, 'end:', iend
+    lon = ilon[iybeg:iyend,ixbeg:ixend]
+    lat = ilat[iybeg:iyend,ixbeg:ixend]
+
+#    print 'lon _______'
+#    print lon
+#    print 'lat _______'
+#    print lat
+
+# Find point
+    difflonlat = np.sqrt((lon-lonv)**2. + (lat-latv)**2.)
+    mindiffLl = np.min(difflonlat)
+    ilatlon = index_mat(difflonlat, mindiffLl)
+
+    ilatlon[0] = ilatlon[0] + iybeg
+    ilatlon[1] = ilatlon[1] + ixbeg
+
+#    print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
+#    print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]
+#    quit()
+
+    return ilatlon, mindiffLl
+
+def lonlatFind(ilon, ilat, lonv, latv):
+    """ Function to search a given value from a lon, lat 2D data
+      ilon, ilat= original 2D matrices with the longitudes and the latitudes
+      lonv, latv= longitude and latitude to find
+      >>> npt=100
+      >>> lonval0 = np.arange(0.,npt*1.,1.)
+      >>> latval0 = np.arange(0.,npt*1.,1.)
+      >>> lonval = np.zeros((npt,npt), dtype=np.float)
+      >>> latval = np.zeros((npt,npt), dtype=np.float)
+
+      >>> for i in range(npt):
+      >>>     lonval[:,i] = lonval0[i]
+      >>>     latval[i,:] = latval0[i]
+      >>> lonlatFind(lonval, latval, 5.25, 5.25)
+      (array([5, 5]), 0.35355339059327379)
+    """
+    fname = 'lonlatFind'
+
+    dx = ilon.shape[1]
+    dy = ilon.shape[0]
+
+    nlon = np.min(ilon)
+    xlon = np.max(ilon)
+    nlat = np.min(ilat)
+    xlat = np.max(ilat)
+
+    if lonv < nlon or lonv > xlon:
+        print errormsg
+        print '  ' + fname + ': longitude outside data range!!'
+        print '    given value:', lonv,'outside [',nlon,',',xlon,']'
+        quit(-1) 
+    if latv < nlat or latv > xlat:
+        print errormsg
+        print '  ' + fname + ': latitude outside data range!!'
+        print '    given value:', latv,'outside [',nlat,',',xlat,']'
+        quit(-1) 
+
+# Find point
+    difflonlat = np.sqrt((ilon-lonv)**2. + (ilat-latv)**2.)
+    mindiffLl = np.min(difflonlat)
+    ilatlon = index_mat(difflonlat, mindiffLl)
+
+#    print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
+#    print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]
+
+    return ilatlon, mindiffLl
+
 def Partialmap_Entiremap(values, filen, varn):
     """ Function to transform from a partial global map (e.g.: only land points) to an entire one
-      values= [lonmame],[latname],[fillVal],[resolution],[kind]
+      values= [lonmame],[latname],[fillVal],[resolution],[kind],[lonlatProjfile]
         [lonname]: name of the longitude variable
         [latname]: name of the latitude variable
@@ -18085 +18590 @@
          'lonlat': standard lon/lat projection tacking number of points from [dx/dy]res at the Equator
          'lonlat_dxdyFix': projection with a fixed grid distance
+         'Goode': Goode projection
+        [lonlatProjfile]: file with the lon,lat of the desired projection. 'None' to be computed and written on fly
       filen= name of the netCDF file
       varn= name of the variable
@@ -18093 +18600 @@
     fname = 'Partialmap_Entiremap'
 
-    arguments = '[lonmame],[latname],[fillVal],[resolution],[kind]'
-    check_arguments(fname, values, arguments, ',')
-
     if values == 'h':
         print fname + '_____________________________________________________________'
         print Partialmap_Entiremap.__doc__
         quit()
+
+    arguments = '[lonmame],[latname],[fillVal],[resolution],[kind],[lonlatProjfile]'
+    check_arguments(fname, values, arguments, ',')
 
     ofile = 'EntireGlobalMap.nc'
@@ -18115 +18622 @@
     resolution = np.float(values.split(',')[3])
     kind = values.split(',')[4]
-
-# Input reference lonlat map
-    lonlatProjfile = filen.split('.')[0] + '_' + kind + '.nc'
+    Projfile = values.split(',')[5]
+    if Projfile == 'None':
+        lonlatProjfile = None
+    else:
+        lonlatProjfile = Projfile
 
     if not onc.variables.has_key(lonname):
@@ -18136 +18645 @@
     latvs = olat[:]
 
+    Ninpts = len(lonvs)
+
     nlat = np.min(latvs)
     nlon = np.min(lonvs)
@@ -18188 +18699 @@
         quit(-1)
 
-    if not os.path.isfile(lonlatProjfile):
+    if lonlatProjfile is None:
+        lonlatProjfile = kind + '.nc'
         print warnmsg
-        print '  ' + fname + ": no reference map '" + lonlatProjfile + "' found !!"
-        print '  creating it via: lonlatProj(', resolution, ',', resolution,         \
+        print '  ' + fname + ": no reference map !!"
+        print "    creation of '" + lonlatProjfile + "'"
+        print '    creating it via: lonlatProj(', resolution, ',', resolution,       \
           ', ' + kind + ', True)'
         lonmap, latmap = lonlatProj(resolution, resolution, kind, True)
@@ -18197 +18710 @@
 
     oproj = NetCDFFile(lonlatProjfile, 'r')
-    if kind == 'lonlat':
+    if kind == 'Goode':
+        olon = oproj.variables['lon']
+        olat = oproj.variables['lat']
+        lonmap = olon[:]
+        latmap = olat[:]
+        projlon = olon[:]
+        projlat = olat[:]
+        omapvals = oproj.variables['ptid']
+        Ntotvals = omapvals.shape[0] * omapvals.shape[1]
+        Nmaplon = olon.shape[1]
+        Nmaplat = olat.shape[0]
+
+    elif kind == 'lonlat':
         olon = oproj.variables['lon']
         olat = oproj.variables['lat']
@@ -18233 +18758 @@
 ##
     ovar = onc.variables[varn]
-    vartype = ovar.dtype
+    vartype = type(ovar[0])
     varlongname = ovar.long_name
     varunits = ovar.units
 
-    fval = fill_value(vartype, fval0)
+    fval = fillvalue_kind(type(ovar[0]), fval0)
 
 # Final map creation
 ##
-    newnc = NetCDFFile(ofile, 'w')
-
+    if not os.path.isfile(ofile):
+        print '  ' + fname + "File '" + ofile + "' does not exist !!"
+        print '    cration of output file'
+        newnc = NetCDFFile(ofile, 'w')
 # dimensions
-    newdim = newnc.createDimension('lon',Nmaplon)
-    newdim = newnc.createDimension('lat',Nmaplat)
-    newdim = newnc.createDimension('Npts', Ntotvals)
+        newdim = newnc.createDimension('lon',Nmaplon)
+        newdim = newnc.createDimension('lat',Nmaplat)
+        newdim = newnc.createDimension('inpts', Ninpts)
+        newdim = newnc.createDimension('pts', Ntotvals)
 
 # Variables
-    newvar = newnc.createVariable('lon','f8',('lon'))
-    basicvardef(newvar, 'lon', 'Longitudes','degrees East')
-    newvar[:] = lonmap
-    newvar.setncattr('axis', 'X')
-    newvar.setncattr('_CoordinateAxisType', 'Lon')   
-
-    newvar = newnc.createVariable('lat','f8',('lat'))
-    basicvardef(newvar, 'lat', 'Latitudes','degrees North')
-    newvar[:] = latmap
-    newvar.setncattr('axis', 'Y')
-    newvar.setncattr('_CoordinateAxisType', 'Lat')   
+        if kind == 'Goode':
+            newvar = newnc.createVariable('lon','f8',('lat', 'lon'))
+        else:
+            newvar = newnc.createVariable('lon','f8',('lon'))
+        basicvardef(newvar, 'lon', 'Longitudes','degrees East')
+        newvar[:] = lonmap
+        newvar.setncattr('axis', 'X')
+        newvar.setncattr('_CoordinateAxisType', 'Lon')   
+
+        if kind == 'Goode':
+            newvar = newnc.createVariable('lat','f8',('lat','lon'))
+        else:
+            newvar = newnc.createVariable('lat','f8',('lat'))
+        basicvardef(newvar, 'lat', 'Latitudes','degrees North')
+        newvar[:] = latmap
+        newvar.setncattr('axis', 'Y')
+        newvar.setncattr('_CoordinateAxisType', 'Lat')   
+
+        newvarinpt = newnc.createVariable('locinpt','i',('inpts'))
+        basicvardef(newvarinpt, 'locinpt', 'input point located: 0: no, 1: yes','-')
+
+        newvarindiff = newnc.createVariable('locindiff','f4',('inpts'))
+        basicvardef(newvarindiff, 'locindiff', 'distance between input point and its final location','degree')
 
 # map variable
-    Lmapvalsshape = len(omapvals.shape)
-    if Lmapvalsshape == 1:
-        newvar = newnc.createVariable(varn,vartype,('Npts'), fill_value=fval)    
+        Lmapvalsshape = len(omapvals.shape)
+        if Lmapvalsshape == 1:
+            newvar = newnc.createVariable(varn, nctype(vartype), ('pts'), fill_value=fval)
+        else:
+            newvar = newnc.createVariable(varn, nctype(vartype), ('lat','lon'), fill_value=fval)
+
+        basicvardef(newvar, varn, varlongname, varunits)
+        newvar.setncattr('coordinates', 'lon lat')
+
+        if Lmapvalsshape == 1:
+            newvarin = newnc.createVariable('inpts', 'i4', ('pts'))    
+        else:
+            newvarin = newnc.createVariable('inpts', 'i4', ('lat','lon'))
+
+        basicvardef(newvarin, 'inpts', 'Equivalent point from the input source', '-')
+        newvar.setncattr('coordinates', 'lon lat')
+
     else:
-        newvar = newnc.createVariable(varn,vartype,('lat','lon'), fill_value=fval)
-
-    basicvardef(newvar, varn, varlongname, varunits)
-    newvar.setncattr('coordinates', 'lon lat')
-
-    print '  ' + fname + ': re-locating:',Ntotvals,'points...'
-    if kind == 'lonlat':
+        print '  ' + fname + "File '" + ofile + "' exists !!"
+        print '    reading values from file'
+        newnc = NetCDFFile(ofile, 'a')
+        newvar = newnc.variables[varn]
+        newvarin = newnc.variables['inpts']
+        newvarinpt = newnc.variables['locinpt']
+        newvarindiff = newnc.variables['locindiff']
+
+    amsk = np.arange(3)
+    amsk = ma.masked_equal(amsk, 0)
+
+#    fraclon = projlon[::Nmaplon*0.1]
+#    fraclat = projlat[::Nmaplat*0.1]
+
+#    print 'fraclon________________', fraclon.shape
+#    print fraclon
+
+#    print 'fraclat________________', fraclat.shape
+#    print fraclat
+
+    print '  ' + fname + ': re-locating:',Ninpts,'points...'
+    if kind == 'Goode':
+        for iv in range(Ninpts):
+            if newvarinpt[iv] == 0:
+
+                difflonlat = np.sqrt((projlon-lonvs[iv])**2. + (projlat-latvs[iv])**2.)
+                mindiffLl = np.min(difflonlat)
+                ilatlon = index_mat(difflonlat, mindiffLl)
+#                ilatlon2, mindiffLl2 = CoarselonlatFind(projlon, projlat, lonvs[iv], latvs[iv], .1)
+
+#                if ilatlon[0] != ilatlon2[0] or ilatlon[1] != ilatlon2[1]:
+#                    print ilatlon, '<', ilatlon2
+#                    print 'diff:', mindiffLl, '<', mindiffLl2
+#                    quit(-1)
+
+                if type(newvar[ilatlon[0],ilatlon[1]]) != type(amsk): 
+                    print errormsg
+                    print '  ' + main + ': point', projlon[ilatlon[0],ilatlon[1]],  \
+                      ',', projlat[ilatlon[0],ilatlon[1]], 'already filled!!'
+                    print '    value:', newvar[ilatlon[0],ilatlon[1]], 'distance:', \
+                      newvarindiff[newvarin[ilatlon[0],ilatlon[1]]]
+
+                    print '   iv:', iv, 'lon:', lonvs[iv],'lat:',latvs[iv]
+                    print '   mindiffLl:',mindiffLl,'ilatlon:',ilatlon
+                    quit(-1)
+#            print '  Lluis; ' + fname + ' iv:', iv, 'lon:', lonvs[iv],'lat:',latvs[iv]
+#            print '  Lluis; ' + fname + ' mindiffLl:',mindiffLl,'ilatlon:',ilatlon
+#            print '  Lluis; ' + fname + ' A found lon:',projlon[ilatlon[0], ilatlon[1]], 'lat:', projlat[ilatlon[0], ilatlon[1]]
+#            quit()
+
+                percendone(iv,Ninpts,0.5,'done:')
+                if mindiffLl > mindiff:
+                    print errormsg
+                    print '  ' + fname + ': for point #', iv,'lon,lat in ' +         \
+                      'incomplet map:', lonvs[iv], ',', latvs[iv], 'there is not ' + \
+                      'a set of lon,lat in the completed map closer than: ', mindiff,\
+                      '!!'
+                    print '    minimum difference:', mindiffLl
+                    quit()
+
+                if ilatlon[0] >= 0 and ilatlon[1] >= 0:
+                    newvar[ilatlon[0],ilatlon[1]] = ovar[iv]
+                    newvarin[ilatlon[0],ilatlon[1]] = iv
+                    newvarinpt[iv] = 1
+                    newvarindiff[iv] = mindiffLl
+#                    print 'Lluis iv:', newvarin[ilatlon[0],ilatlon[1]],              \
+#                      'localized:', newvarinpt[iv], 'values:',                       \
+#                      newvar[ilatlon[0],ilatlon[1]], 'invalues:', ovar[iv],          \
+#                      'mindist:', newvarindiff[iv], 'point:',ilatlon
+                else:
+                    print errormsg
+                    print '  ' + fname + ': point iv:', iv, 'at', lonvs[iv], ',' ,   \
+                      latvs[iv],' not relocated !!'
+                    print '    mindiffl:', mindiffLl, 'ilon:', ilatlon[1],           \
+                      'ilatlon:', ilatlon[1]
+                    quit(-1)
+
+                if np.mod(iv,100) == 0: 
+                    newnc.sync()
+#                    print '  ' + fname + 'values localized:', newvar[:].compressed()
+#                if iv > 10:
+#                    newnc.sync()
+#                    newnc.close()
+#                    quit(-1)
+
+    elif kind == 'lonlat':
         for iv in range(Ntotvals):
             difflat = np.abs(projlat - latvs[iv])
@@ -18299 +18932 @@
 
         newnc.sync()
-                        
 
     elif kind == 'lonlat_dxdyFix':
@@ -18358 +18990 @@
     return 
 
+#Partialmap_Entiremap('longitude,latitude,std,5000.,Goode,Goode_5km.nc', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/#carteveg5km.nc', 'vegetation_map')
 #Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat_dxdyFix', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
 #Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')