Changeset 1484 in lmdz_wrf for trunk

Timestamp:

Apr 1, 2017, 10:31:26 PM (8 years ago)

Author:

lfita

Message:

Adding:

`lonlat_polygon': Function to define a lon/lat region giving the coordinates of the vertexs of a given polygon

Location:

trunk/tools

Files:

• documentation/ncmanage/lonlat_polygon.html (added)
• documentation/ncmanage/ncmanage.html (modified) (1 diff)
• nc_var.py (modified) (4 diffs)
• nc_var_tools.py (modified) (2 diffs)

trunk/tools/documentation/ncmanage/ncmanage.html

@@ -10 +10 @@
       <A CLASS="lc" HREF="curve_section.html" TARGET="value">curve_section</A><BR>
       <A CLASS="lc" HREF="field_stats.html" TARGET="value">field_stats</A><BR>
+      <A CLASS="lc" HREF="lonlat_polygon.html" TARGET="value">lonlat_polygon</A><BR>
     </DIV>
     <DIV CLASS="valmenu">

trunk/tools/nc_var.py

@@ -8 +8 @@
 #   Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0)
 #
+
+## e.g. # nc_var.py -o lonlat_polygon -f wrfout_d01_1995-01-01_00:00:00 -S star.dat -v XLONG,XLAT
 
 ## e.g. ccrc468-17 # ./nc_var.py -v time -f 123/CCRC_NARCliM_Sydney_All_1990-1999_pr10max.nc -o out -S 1:-1
@@ -102 +104 @@
 # ivattrs: Give all the attributes of a variable and its type
 # LMDZ_toCF: Function to pass a LMDZ original file to CF-conventions
+# lonlat_polygon: Function to define a lon/lat region giving the coordinates of the vertexs of a given polygon 
 # maskvar: Function to mask a variable using a mask. It is assumed that mask[...,dimM,dimJ,dimK] and 
 # model_characteristics: Function to provide major characterisitcs of a given model output
@@ -171 +174 @@
   'getvalues_lonlat', 'getvars_tofile', 'grattr',                                    \
   'grmattr', 'idims', 'igattrs', 'increaseDimvar', 'isgattrs',                       \
-  'isvattrs', 'ivars', 'ivattrs', 'LMDZ_toCF', 'maskvar', 'model_characteristics',   \
+  'isvattrs', 'ivars', 'ivattrs', 'LMDZ_toCF', 'lonlat_polygon', 'maskvar',          \
+  'model_characteristics',                                                           \
   'mthDYNAMICO_toCF', 'ncreplace', 'ncstepdiff', 'netcdf_concatenation',             \
   'netcdf_fold_concatenation',                                                       \
@@ -372 +376 @@
 elif oper == 'LMDZ_toCF':
     ncvar.LMDZ_toCF(opts.values, opts.ncfile)
+elif oper == 'lonlat_polygon':
+    ncvar.lonlat_polygon(opts.values, opts.ncfile, opts.varname)
 elif oper == 'maskvar':
     ncvar.maskvar(opts.values, opts.ncfile, opts.varname)

trunk/tools/nc_var_tools.py

@@ -89 +89 @@
 # load_ncvariable_lastdims: Function to load the last [Nlastdims] dimensions of a variable [varn] from a netCDF object at the other dimensions at [prevdims]
 # lonlatProj: Function to compute longitudes and latitudes for a given projection following subroutine calc_resolution_in from ORCHIDEE/src_global/module_InterpWeight
+# lonlat_polygon: Function to define a lon/lat region giving the coordinates of the vertexs of a given polygon 
 # maskvar: Function to mask a variable using a mask. It is assumed that mask[...,dimM,dimJ,dimK] and var[...,dimM,dimJ,dimK] share the last dimensions
 # ModelChar: Class object to determine major characterisitcs of a given model output
@@ -19135 +19136 @@
     return varunits
 
+def lonlat_polygon(values,ncfile,varn):
+    """ Function to define a lon/lat region giving the coordinates of the vertexs of a given polygon 
+      lonlat_polygon(values,ncfile,varn)
+      values= [asciipolygfile]
+        [asciipoligfile]: ASCII file with the clocwise consecutive poligon vertexs as: ('#' for coments)
+          [lon1] [lat1]
+          (...)
+          [lonN] [latN]
+      ncfile= name of the file with the variable to fill
+      varn= [lonvname],[latvname] name of the longitude and latitude variable in the file
+    """
+    import matplotlib.path as mpltPath
+    fname = 'lonlat_polygon'
+
+    ofile = fname + '.nc'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print lonlat_polygon.__doc__
+        quit()
+
+    asciipolygfile = values
+
+    lonvname = varn.split(',')[0]
+    latvname = varn.split(',')[1]
+
+    if not os.path.isfile(asciipolygfile):
+        print errormsg
+        print '  ' + fname + ": ASCII vertex polygon file '" + asciipolygfile +      \
+          "' does not exist !!"
+        quit(-1)
+
+    # getting lon,lat
+    onc = NetCDFFile(ncfile, 'r')
+    if not onc.variables.has_key(lonvname):
+        print errormsg
+        print '  ' + fname + ": netCDF file '" + ncfile + "' does not have variable"+\
+          " longitude '" + lonvname + "' !!"
+        quit(-1)
+    if not onc.variables.has_key(latvname):
+        print errormsg
+        print '  ' + fname + ": netCDF file '" + ncfile + "' does not have variable"+\
+          " latitude '" + latvname + "' !!"
+        quit(-1)
+
+    lon, lat = gen.lonlat2D(onc.variables[lonvname], onc.variables[latvname])
+
+    onc.close()
+
+    # getting vertices
+    overtf = open(asciipolygfile, 'r')
+    lLvertices = []
+    for line in overtf:
+        if line[0:1] != '#' and len(line) > 3:
+            linevals = line.replace('\n','').split(' ')
+            lLvertices.append([np.float(linevals[0]), np.float(linevals[1])])
+    overtf.close()
+
+    # Getting ij in lon,lat space of the vertices
+    Nvertices = len(lLvertices)
+    ijvertices = []
+    for iv in range(Nvertices):
+        vertex = lLvertices[iv]
+        ijvertices.append(gen.lonlatFind(lon,lat,vertex[0],vertex[1])[0])
+        print iv, vertex, ijvertices[iv]
+
+    plane = np.zeros((lon.shape), dtype=int)
+    # Getting the points which are corssed by the faces of the polygon
+    for iiv in range(Nvertices):
+        lL = lLvertices[iiv]
+        l = lL[0]
+        L = lL[1]
+        ij = ijvertices[iiv]
+        iv = ij[1]
+        jv = ij[0]
+
+        # Next point
+        if iiv == Nvertices - 1:
+            nextiiv = 0
+        else:
+            nextiiv = iiv+1
+        
+        lLdest = lLvertices[nextiiv]
+        ldest = lLdest[0]
+        Ldest = lLdest[1]
+        ijdest = ijvertices[nextiiv]
+        ivdest = ijdest[1]
+        jvdest = ijdest[0]
+
+        # connecting line y = a + by
+        dlon = ldest - l
+        dlat = Ldest - L
+        if dlon != 0.:
+            b = 1.*dlat/dlon
+            sdx = dlon/np.abs(dlon)
+        else:
+            b = 0.
+            sdx = 1
+        if dlat != 0.:
+            sdy = dlat/np.abs(dlat)
+        else:
+            sdy = 1
+
+        # Get more grid points (just in case)
+        Nlon = np.abs(int((ivdest - iv)*1.5))
+        Nlat = np.abs(int((jvdest - jv)*1.5))
+
+        if dlon != 0. and Nlon == 0: Nlon = Nlat
+        if Nlat > Nlon: Nlon = Nlat
+        
+        #print iiv, ',', lL, lLdest, ':', dlon, dlat, b, 'Npts:', Nlon, Nlat
+        # Moving across the line folloing lon,lat
+        if dlon != 0.:
+            for iil in range(Nlon + 1):
+                ilon = dlon*iil/Nlon
+                ilat = L + ilon*b
+                [iy,ix], dist = gen.lonlatFind(lon,lat,l+ilon,ilat)
+                #print '  ',iil,l+ilon,ilat,ix,iy,dist
+                plane[iy,ix] = 1
+        else:
+            ilon = l
+            for iiL in range(Nlat + 1):
+                ilat = L + dlat*iiL/Nlat
+                [iy,ix], dist = gen.lonlatFind(lon,lat,ilon,ilat)
+                #print '  ',iiL,ilon,ilat,ix,iy,dist
+                plane[iy,ix] = 1
+
+        #print '    last point of the face:', l+ilon, ilat, 'destiny:', ldest, Ldest
+
+    # Filling the polygon
+    # FROM: http://stackoverflow.com/questions/36399381/whats-the-fastest-way-of-checking-if-a-point-is-inside-a-polygon-in-python
+    lonvals = lon.flatten()
+    latvals = lat.flatten()
+    lonlatvals = zip(lonvals, latvals)
+    
+    path = mpltPath.Path(lLvertices)
+    inside = path.contains_points(lonlatvals)
+    insidemat = inside.reshape(lon.shape)
+
+    insideones = np.where(insidemat, 1, 0)
+
+    # output file
+    newnc = NetCDFFile(ofile, 'w')
+
+    # dimensions
+    newdim = newnc.createDimension('lon',lon.shape[1])
+    newdim = newnc.createDimension('lat',lon.shape[0])
+    newdim = newnc.createDimension('vertex',Nvertices)
+    
+    # variable-dimensions
+    newvar = newnc.createVariable('lon', 'f8', ('lat', 'lon'))
+    newvar[:] = lon
+    basicvardef(newvar,'longitude','Longitude','degrees_East')
+    newvar.setncattr('axis', 'X')
+    newvar.setncattr('_CoordinateAxisType', 'Lon')
+    newvar = newnc.createVariable('lat', 'f8', ('lat', 'lon'))
+    newvar[:] = lat
+    basicvardef(newvar,'latitude','Latitude','degrees_North')
+    newvar.setncattr('axis', 'Y')
+    newvar.setncattr('_CoordinateAxisType', 'Lat')
+    lnewvar = newnc.createVariable('vertexlon', 'f8', ('vertex'))
+    basicvardef(lnewvar,'vertexlon','longitude coordinate of the vertex of polygon','degrees_East')
+    Lnewvar = newnc.createVariable('vertexlat', 'f8', ('vertex'))
+    basicvardef(Lnewvar,'vertexlat','latitude coordinate of the vertex of polygon','degrees_North')
+    inewvar = newnc.createVariable('ivertex', 'i', ('vertex'))
+    basicvardef(inewvar,'ivertex','coordinate on axis-x of the vertex of polygon','1')
+    jnewvar = newnc.createVariable('jvertex', 'i', ('vertex'))
+    basicvardef(jnewvar,'jvertex','coordinate on axis-y of the vertex of polygon','1')
+    for i in range(Nvertices):
+        lL = lLvertices[i]
+        ij = ijvertices[i]
+        lnewvar[i] = lL[0]
+        Lnewvar[i] = lL[1]
+        inewvar[i] = ij[0]
+        jnewvar[i] = ij[1]
+    newattr = set_attribute(lnewvar, 'src_polygon_file', asciipolygfile)
+    newattr = set_attribute(Lnewvar, 'src_polygon_file', asciipolygfile)
+    newattr = set_attribute(inewvar, 'src_projection_file', ncfile)
+    newattr = set_attribute(jnewvar, 'src_projection_file', ncfile)
+
+    # variable
+    newvar = newnc.createVariable('polygon_countour', 'i', ('lat', 'lon'))
+    newvar[:] = plane
+    basicvardef(newvar,'polygon_countour','countour of the Polygon','1')
+    newvar.setncattr('coordinates', 'lon lat')
+    newvar = newnc.createVariable('polygon_sfc', 'i', ('lat', 'lon'))
+    newvar[:] = insideones
+    basicvardef(newvar,'polygon_surface','surface of the Polygon','1')
+    newvar.setncattr('coordinates', 'lon lat')
+
+    newnc.sync()
+
+    # Global attributes
+    add_global_PyNCplot(newnc, main, fname, '1.0')
+
+    newnc.close()
+
+    print fname + ": succesfull written of file with polygon surface '" + ofile + "' !!"
+
+    return
+
+#lonlat_polygon('/home/lluis/PY/square.dat','/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00', 'XLONG,XLAT')
+#lonlat_polygon('/home/lluis/PY/star.dat','/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00', 'XLONG,XLAT')
+
+
 #quit()