Python NetCDFFile.units方法代码示例

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import units [as 别名]
def _convert_dem_from_ascii2netcdf(name_in, name_out = None,
                                   verbose = False):
    """Read Digital Elevation model from the following ASCII format (.asc)

    Internal function. See public function convert_dem_from_ascii2netcdf
    for details.
    """

    import os
    from anuga.file.netcdf import NetCDFFile

    root = name_in[:-4]

    # Read Meta data
    if verbose: log.critical('Reading METADATA from %s' % (root + '.prj'))

    metadatafile = open(root + '.prj')
    metalines = metadatafile.readlines()
    metadatafile.close()

    L = metalines[0].strip().split()
    assert L[0].strip().lower() == 'projection'
    projection = L[1].strip()                   #TEXT

    L = metalines[1].strip().split()
    assert L[0].strip().lower() == 'zone'
    zone = int(L[1].strip())

    L = metalines[2].strip().split()
    assert L[0].strip().lower() == 'datum'
    datum = L[1].strip()                        #TEXT

    L = metalines[3].strip().split()
    assert L[0].strip().lower() == 'zunits'     #IGNORE
    zunits = L[1].strip()                       #TEXT

    L = metalines[4].strip().split()
    assert L[0].strip().lower() == 'units'
    units = L[1].strip()                        #TEXT

    L = metalines[5].strip().split()
    assert L[0].strip().lower() == 'spheroid'   #IGNORE
    spheroid = L[1].strip()                     #TEXT

    L = metalines[6].strip().split()
    assert L[0].strip().lower() == 'xshift'
    false_easting = float(L[1].strip())

    L = metalines[7].strip().split()
    assert L[0].strip().lower() == 'yshift'
    false_northing = float(L[1].strip())

    if name_in[-4:] != '.asc':
        raise IOError('Input file %s should be of type .asc.' % name_in)

    #Read DEM data
    datafile = open(name_in)

    if verbose: log.critical('Reading DEM from %s' % (name_in))

    lines = datafile.readlines()
    datafile.close()

    if verbose: log.critical('Got %d lines' % len(lines))

    ncols = int(lines[0].split()[1].strip())
    nrows = int(lines[1].split()[1].strip())

    # Do cellsize (line 4) before line 2 and 3
    cellsize = float(lines[4].split()[1].strip())

    # Checks suggested by Joaquim Luis
    # Our internal representation of xllcorner
    # and yllcorner is non-standard.
    xref = lines[2].split()
    if xref[0].strip() == 'xllcorner':
        xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
    elif xref[0].strip() == 'xllcenter':
        xllcorner = float(xref[1].strip())
    else:
        msg = 'Unknown keyword: %s' % xref[0].strip()
        raise Exception, msg

    yref = lines[3].split()
    if yref[0].strip() == 'yllcorner':
        yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
    elif yref[0].strip() == 'yllcenter':
        yllcorner = float(yref[1].strip())
    else:
        msg = 'Unknown keyword: %s' % yref[0].strip()
        raise Exception, msg

    NODATA_value = int(float(lines[5].split()[1].strip()))

    assert len(lines) == nrows + 6

    if name_out == None:
        netcdfname = name_in[:-4]+'.dem'
    else:
        netcdfname = name_out + '.dem'
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import units [as 别名]
def dem2dem(name_in, stencil, cellsize_new, name_out=None,
                 verbose=False):
    """Read Digitial Elevation model from the following NetCDF format (.dem)

    Example:

    ncols         3121
    nrows         1800
    xllcorner     722000
    yllcorner     5893000
    cellsize      25
    NODATA_value  -9999
    138.3698 137.4194 136.5062 135.5558 ..........

    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
    """

    import os
    from anuga.file.netcdf import NetCDFFile

    if name_in[-4:] != '.dem':
        raise IOError('Input file %s should be of type .dem.' % name_in)

    if name_out != None and basename_out[-4:] != '.dem':
        raise IOError('Input file %s should be of type .dem.' % name_out)

    #Open existing netcdf file to read
    infile = NetCDFFile(name_in, netcdf_mode_r)

    if verbose: log.critical('Reading DEM from %s' % inname)

    # Read metadata (convert from numpy.int32 to int where appropriate)
    ncols = int(infile.ncols)
    nrows = int(infile.nrows)
    xllcorner = infile.xllcorner
    yllcorner = infile.yllcorner
    cellsize = int(infile.cellsize)
    NODATA_value = int(infile.NODATA_value)
    zone = int(infile.zone)
    false_easting = infile.false_easting
    false_northing = infile.false_northing
    projection = infile.projection
    datum = infile.datum
    units = infile.units

    dem_elevation = infile.variables['elevation']

    #Get output file name
    if name_out == None:
        outname = name_in[:-4] + '_' + repr(cellsize_new) + '.dem'
    else:
        outname = name_out

    if verbose: log.critical('Write decimated NetCDF file to %s' % outname)

    #Determine some dimensions for decimated grid
    (nrows_stencil, ncols_stencil) = stencil.shape
    x_offset = ncols_stencil / 2
    y_offset = nrows_stencil / 2
    cellsize_ratio = int(cellsize_new / cellsize)
    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio

    #print type(ncols_new), ncols_new
    
    #Open netcdf file for output
    outfile = NetCDFFile(outname, netcdf_mode_w)

    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'NetCDF DEM format for compact and portable ' \
                          'storage of spatial point data'

    #Georeferencing
    outfile.zone = zone
    outfile.projection = projection
    outfile.datum = datum
    outfile.units = units

    outfile.cellsize = cellsize_new
    outfile.NODATA_value = NODATA_value
    outfile.false_easting = false_easting
    outfile.false_northing = false_northing

    outfile.xllcorner = xllcorner + (x_offset * cellsize)
    outfile.yllcorner = yllcorner + (y_offset * cellsize)
    outfile.ncols = ncols_new
    outfile.nrows = nrows_new

    # dimension definition
    #print nrows_new, ncols_new, nrows_new*ncols_new
    #print type(nrows_new), type(ncols_new), type(nrows_new*ncols_new)
    outfile.createDimension('number_of_points', nrows_new*ncols_new)

    # variable definition
    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))

    # Get handle to the variable
    elevation = outfile.variables['elevation']

#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import units [as 别名]
def _dem2pts(name_in, name_out=None, verbose=False,
            easting_min=None, easting_max=None,
            northing_min=None, northing_max=None):
    """Read Digitial Elevation model from the following NetCDF format (.dem)

    Internal function. See public function dem2pts for details.
    """

    # FIXME: Can this be written feasibly using write_pts?

    import os
    from anuga.file.netcdf import NetCDFFile

    root = name_in[:-4]

    if name_in[-4:] == '.asc':
        intermediate = root + '.dem'
        if verbose:
            log.critical('Preconvert %s from asc to %s' % \
                                    (name_in, intermediate))
        asc2dem(name_in)
        name_in = intermediate
    elif name_in[-4:] != '.dem':
        raise IOError('Input file %s should be of type .asc or .dem.' % name_in)

    if name_out != None and basename_out[-4:] != '.pts':
        raise IOError('Input file %s should be of type .pts.' % name_out)

    # Get NetCDF
    infile = NetCDFFile(name_in, netcdf_mode_r) 

    if verbose: log.critical('Reading DEM from %s' % (name_in))

    ncols = int(infile.ncols)
    nrows = int(infile.nrows)
    xllcorner = float(infile.xllcorner)  # Easting of lower left corner
    yllcorner = float(infile.yllcorner)  # Northing of lower left corner
    cellsize = float(infile.cellsize)
    NODATA_value = float(infile.NODATA_value)

    dem_elevation = infile.variables['elevation']

    zone = int(infile.zone)
    false_easting = float(infile.false_easting)
    false_northing = float(infile.false_northing)

    #print ncols, nrows, xllcorner,yllcorner, cellsize, NODATA_value, zone


    # Text strings
    projection = infile.projection
    datum = infile.datum
    units = infile.units

    #print projection, datum, units

    # Get output file
    if name_out == None:
        ptsname = root + '.pts'
    else:
        ptsname = name_out

    if verbose: log.critical('Store to NetCDF file %s' % ptsname)

    # NetCDF file definition
    outfile = NetCDFFile(ptsname, netcdf_mode_w)

    # Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'NetCDF pts format for compact and portable ' \
                          'storage of spatial point data'

    # Assign default values
    if easting_min is None: easting_min = xllcorner
    if easting_max is None: easting_max = xllcorner + ncols*cellsize
    if northing_min is None: northing_min = yllcorner
    if northing_max is None: northing_max = yllcorner + nrows*cellsize


    #print easting_min, easting_max, northing_min, northing_max

    # Compute offsets to update georeferencing
    easting_offset = xllcorner - easting_min
    northing_offset = yllcorner - northing_min

    # Georeferencing
    outfile.zone = zone
    outfile.xllcorner = easting_min # Easting of lower left corner
    outfile.yllcorner = northing_min # Northing of lower left corner
    outfile.false_easting = false_easting
    outfile.false_northing = false_northing

    outfile.projection = projection
    outfile.datum = datum
    outfile.units = units

    # Grid info (FIXME: probably not going to be used, but heck)
    outfile.ncols = ncols
    outfile.nrows = nrows

#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import units [as 别名]
    def _generic_dem2pts(
        self,
        name_in,
        name_out=None,
        quantity_name=None,
        verbose=False,
        easting_min=None,
        easting_max=None,
        northing_min=None,
        northing_max=None,
    ):
        """Read raster from the following NetCDF format (.dem)

        Internal function. See public function generic_dem2pts for details.
        """

        # FIXME: Can this be written feasibly using write_pts?

        import os
        from anuga.file.netcdf import NetCDFFile

        root = name_in[:-4]

        if name_in[-4:] == ".asc":
            intermediate = root + ".dem"
            if verbose:
                log.critical("Preconvert %s from asc to %s" % (name_in, intermediate))
            asc2dem(name_in)
            name_in = intermediate
        elif name_in[-4:] != ".dem":
            raise IOError("Input file %s should be of type .asc or .dem." % name_in)

        if name_out != None and basename_out[-4:] != ".pts":
            raise IOError("Input file %s should be of type .pts." % name_out)

        # Get NetCDF
        infile = NetCDFFile(name_in, netcdf_mode_r)

        if verbose:
            log.critical("Reading raster from %s" % (name_in))

        ncols = int(infile.ncols)
        nrows = int(infile.nrows)
        xllcorner = float(infile.xllcorner)  # Easting of lower left corner
        yllcorner = float(infile.yllcorner)  # Northing of lower left corner
        cellsize = float(infile.cellsize)
        NODATA_value = float(infile.NODATA_value)

        dem_elevation = infile.variables[quantity_name]

        zone = int(infile.zone)
        false_easting = float(infile.false_easting)
        false_northing = float(infile.false_northing)

        # print ncols, nrows, xllcorner,yllcorner, cellsize, NODATA_value, zone

        # Text strings
        projection = infile.projection
        datum = infile.datum
        units = infile.units

        # print projection, datum, units

        # Get output file
        if name_out == None:
            ptsname = root + ".pts"
        else:
            ptsname = name_out

        if verbose:
            log.critical("Store to NetCDF file %s" % ptsname)

        # NetCDF file definition
        outfile = NetCDFFile(ptsname, netcdf_mode_w)

        # Create new file
        outfile.institution = "Geoscience Australia"
        outfile.description = "NetCDF pts format for compact and portable " "storage of spatial point data"

        # Assign default values
        if easting_min is None:
            easting_min = xllcorner
        if easting_max is None:
            easting_max = xllcorner + ncols * cellsize
        if northing_min is None:
            northing_min = yllcorner
        if northing_max is None:
            northing_max = yllcorner + nrows * cellsize

        # print easting_min, easting_max, northing_min, northing_max

        # Compute offsets to update georeferencing
        easting_offset = xllcorner - easting_min
        northing_offset = yllcorner - northing_min

        # Georeferencing
        outfile.zone = zone
        outfile.xllcorner = easting_min  # Easting of lower left corner
        outfile.yllcorner = northing_min  # Northing of lower left corner
        outfile.false_easting = false_easting
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import units [as 别名]
    def test_decimate_dem(self):
        """Test decimation of dem file
        """

        import os
        from anuga.file.netcdf import NetCDFFile

        # Write test dem file
        root = "decdemtest"

        filename = root + ".dem"
        fid = NetCDFFile(filename, netcdf_mode_w)

        fid.institution = "Geoscience Australia"
        fid.description = "NetCDF DEM format for compact and portable " + "storage of spatial point data"

        nrows = 15
        ncols = 18

        fid.ncols = ncols
        fid.nrows = nrows
        fid.xllcorner = 2000.5
        fid.yllcorner = 3000.5
        fid.cellsize = 25
        fid.NODATA_value = -9999

        fid.zone = 56
        fid.false_easting = 0.0
        fid.false_northing = 0.0
        fid.projection = "UTM"
        fid.datum = "WGS84"
        fid.units = "METERS"

        fid.createDimension("number_of_points", nrows * ncols)

        fid.createVariable("elevation", netcdf_float, ("number_of_points",))

        elevation = fid.variables["elevation"]

        elevation[:] = num.arange(nrows * ncols)

        fid.close()

        # generate the elevation values expected in the decimated file
        ref_elevation = [
            (0 + 1 + 2 + 18 + 19 + 20 + 36 + 37 + 38) / 9.0,
            (4 + 5 + 6 + 22 + 23 + 24 + 40 + 41 + 42) / 9.0,
            (8 + 9 + 10 + 26 + 27 + 28 + 44 + 45 + 46) / 9.0,
            (12 + 13 + 14 + 30 + 31 + 32 + 48 + 49 + 50) / 9.0,
            (72 + 73 + 74 + 90 + 91 + 92 + 108 + 109 + 110) / 9.0,
            (76 + 77 + 78 + 94 + 95 + 96 + 112 + 113 + 114) / 9.0,
            (80 + 81 + 82 + 98 + 99 + 100 + 116 + 117 + 118) / 9.0,
            (84 + 85 + 86 + 102 + 103 + 104 + 120 + 121 + 122) / 9.0,
            (144 + 145 + 146 + 162 + 163 + 164 + 180 + 181 + 182) / 9.0,
            (148 + 149 + 150 + 166 + 167 + 168 + 184 + 185 + 186) / 9.0,
            (152 + 153 + 154 + 170 + 171 + 172 + 188 + 189 + 190) / 9.0,
            (156 + 157 + 158 + 174 + 175 + 176 + 192 + 193 + 194) / 9.0,
            (216 + 217 + 218 + 234 + 235 + 236 + 252 + 253 + 254) / 9.0,
            (220 + 221 + 222 + 238 + 239 + 240 + 256 + 257 + 258) / 9.0,
            (224 + 225 + 226 + 242 + 243 + 244 + 260 + 261 + 262) / 9.0,
            (228 + 229 + 230 + 246 + 247 + 248 + 264 + 265 + 266) / 9.0,
        ]

        # generate a stencil for computing the decimated values
        stencil = num.ones((3, 3), num.float) / 9.0

        dem2dem(filename, stencil=stencil, cellsize_new=100)

        # Open decimated NetCDF file
        fid = NetCDFFile(root + "_100.dem", netcdf_mode_r)

        # Get decimated elevation
        elevation = fid.variables["elevation"]

        # Check values
        assert num.allclose(elevation, ref_elevation)

        # Cleanup
        fid.close()

        os.remove(root + ".dem")
        os.remove(root + "_100.dem")