Read GeoTiff raster data and perform log transformation.
496616
def main():
    
    input_tif = "../tests/data/Jamaica_dem.tif"
    output_tif = "../tests/data/tmp_results/log_dem.tif"
    rst = RasterUtilClass.read_raster(input_tif)
    # raster data (with noDataValue as numpy.nan) as numpy array
    rst_valid = rst.validValues
    output_data = np.log(rst_valid)
    # write output raster
    RasterUtilClass.write_gtiff_file(output_tif, rst.nRows, rst.nCols, output_data, rst.geotrans,
                                     rst.srs, rst.noDataValue, rst.dataType)
738178
def files_log_graph(self, stream):
        '''  '''
        for row in list(stream):

            # dataframes['files_log'][['md5','mime_type','missing_bytes','rx_hosts','source','tx_hosts']]
            
            # If the mime-type is interesting add the uri and the host->uri->host relationships
            if row['mime_type'] not in self.exclude_mime_types:

                # Check for weird conditions
                if (row['total_bytes'] == '-'):
                    continue
                if ('-' in row['md5']):
                    continue

                # Check for missing bytes and small file
                if row['missing_bytes']:
                    labels = ['missing', 'file']
                elif row['total_bytes'] < 50*1024:
                    labels = ['small','file']
                else:
                    labels = ['file']

                # Make the file node name kewl
                name = '%6s  %s  %.0f-KB' % (row['md5'][:6], row['mime_type'], row['total_bytes']/1024.0)
                if row['missing_bytes']:
                    name += '*'
                name = name.replace('application/','')
                
                # Add the file node
                self.add_node(row['md5'], name, labels)

                # Add the tx_host
                self.add_node(row['tx_hosts'], row['tx_hosts'], ['host'])

                # Add the file->tx_host relationship
                self.add_rel(row['tx_hosts'], row['md5'], 'file')
903721
def run_all(logdir, verbose=False, num_summaries=400):
  """"""
  del verbose

  tf.compat.v1.set_random_seed(0)

  k = tf.compat.v1.placeholder(tf.float32)

  # Make a normal distribution, with a shifting mean
  mean_moving_normal = tf.random.normal(shape=[1000], mean=(5*k), stddev=1)
  # Record that distribution into a histogram summary
  histogram_summary.op("normal/moving_mean",
                       mean_moving_normal,
                       description="A normal distribution whose mean changes "
                                   "over time.")

  # Make a normal distribution with shrinking variance
  shrinking_normal = tf.random.normal(shape=[1000], mean=0, stddev=1-(k))
  # Record that distribution too
  histogram_summary.op("normal/shrinking_variance", shrinking_normal,
                       description="A normal distribution whose variance "
                                   "shrinks over time.")

  # Let's combine both of those distributions into one dataset
  normal_combined = tf.concat([mean_moving_normal, shrinking_normal], 0)
  # We add another histogram summary to record the combined distribution
  histogram_summary.op("normal/bimodal", normal_combined,
                       description="A combination of two normal distributions, "
                                   "one with a moving mean and one with  "
                                   "shrinking variance. The result is a "
                                   "distribution that starts as unimodal and "
                                   "becomes more and more bimodal over time.")

  # Add a gamma distribution
  gamma = tf.random.gamma(shape=[1000], alpha=k)
  histogram_summary.op("gamma", gamma,
                       description="A gamma distribution whose shape "
                                   "parameter, α, changes over time.")

  # And a poisson distribution
  poisson = tf.compat.v1.random_poisson(shape=[1000], lam=k)
  histogram_summary.op("poisson", poisson,
                       description="A Poisson distribution, which only "
                                   "takes on integer values.")

  # And a uniform distribution
  uniform = tf.random.uniform(shape=[1000], maxval=k*10)
  histogram_summary.op("uniform", uniform,
                       description="A simple uniform distribution.")

  # Finally, combine everything together!
  all_distributions = [mean_moving_normal, shrinking_normal,
                       gamma, poisson, uniform]
  all_combined = tf.concat(all_distributions, 0)
  histogram_summary.op("all_combined", all_combined,
                       description="An amalgamation of five distributions: a "
                                   "uniform distribution, a gamma "
                                   "distribution, a Poisson distribution, and "
                                   "two normal distributions.")

  summaries = tf.compat.v1.summary.merge_all()

  # Setup a session and summary writer
  sess = tf.compat.v1.Session()
  writer = tf.summary.FileWriter(logdir)

  # Setup a loop and write the summaries to disk
  N = num_summaries
  for step in xrange(N):
    k_val = step/float(N)
    summ = sess.run(summaries, feed_dict={k: k_val})
    writer.add_summary(summ, global_step=step)
440271
def process(self, raster):
        
        dim = len(raster.shape)
        if dim == 3:
            for dim in range(raster.shape[2]):
                raster[:, :, dim] = self.morph_operation(raster[:, :, dim], self.struct_elem)
        elif dim == 4:
            for time, dim in it.product(range(raster.shape[0]), range(raster.shape[3])):
                raster[time, :, :, dim] = self.morph_operation(raster[time, :, :, dim], self.struct_elem)
        else:
            raise ValueError('Invalid number of dimensions: {}'.format(dim))

        return raster
147524
def read_geo(fid, key):
    """"""
    dsid = GEO_NAMES[key.name]
    add_epoch = False
    if "time" in key.name:
        days = fid["/L1C/" + dsid["day"]].value
        msecs = fid["/L1C/" + dsid["msec"]].value
        data = _form_datetimes(days, msecs)
        add_epoch = True
        dtype = np.float64
    else:
        data = fid["/L1C/" + dsid].value
        dtype = np.float32
    data = xr.DataArray(da.from_array(data, chunks=CHUNK_SIZE),
                        name=key.name, dims=['y', 'x']).astype(dtype)

    if add_epoch:
        data.attrs['sensing_time_epoch'] = EPOCH

    return data
1329668
def fit_all(xy,uv,mode='rscale',center=None,verbose=True):
    
    if mode not in ['general', 'shift', 'rscale']:
        mode = 'rscale'
    if not isinstance(xy,np.ndarray):
        # cast input list as numpy ndarray for fitting
        xy = np.array(xy)
    if not isinstance(uv,np.ndarray):
        # cast input list as numpy ndarray for fitting
        uv = np.array(uv)

    if mode == 'shift':
        logstr = 'Performing "shift" fit'
        if verbose:
            print(logstr)
        else:
            log.info(logstr)
        result = fit_shifts(xy, uv)

    elif mode == 'general':
        logstr = 'Performing "general" fit'
        if verbose:
            print(logstr)
        else:
            log.info(logstr)
        result = fit_general(xy, uv)

    else:
        logstr = 'Performing "rscale" fit'
        if verbose:
            print(logstr)
        else:
            log.info(logstr)
        result = geomap_rscale(xy, uv, center=center)

    return result
964624
def map_data(self, map_name):
    """"""
    with gfile.Open(os.path.join(self.data_dir, "Maps", map_name), "rb") as f:
      return f.read()
1264671
def from_rasters(cls, rasters, relative_to_vrt=True, destination_file=None, nodata=None, mask_band=None):
        """"""
        if isinstance(rasters, list):
            doc = raster_list_vrt(rasters, relative_to_vrt, nodata, mask_band).tostring()
        else:
            doc = raster_collection_vrt(rasters, relative_to_vrt, nodata, mask_band).tostring()
        filename = cls._save_to_destination_file(doc, destination_file)
        return GeoRaster2.open(filename)
647
def _raster_opener(cls, filename, *args, **kwargs):
        """"""
        with rasterio.Env(**cls.get_gdal_env(filename)):
            try:
                return rasterio.open(filename, *args, **kwargs)
            except (rasterio.errors.RasterioIOError, rasterio._err.CPLE_BaseError) as e:
                raise GeoRaster2IOError(e)
391361
def extract(self, point_x, point_y, radius=0):
        '''
        
        '''
        row, col = map_pixel(point_x, point_y, self.x_cell_size, self.y_cell_size,
                             self.xmin, self.ymax)
        col2 = np.abs(radius/self.x_cell_size).astype(int)
        row2 = np.abs(radius/self.y_cell_size).astype(int)
        return GeoRaster(self.raster[max(row-row2, 0):min(row+row2+1, self.shape[0]), \
                        max(col-col2, 0):min(col+col2+1, self.shape[1])], self.geot,
                        nodata_value=self.nodata_value,\
                        projection=self.projection, datatype=self.datatype)
194047
def _reproject(self, eopatch, src_raster):
        
        height, width = src_raster.shape

        dst_raster = np.ones((height, width), dtype=self.raster_dtype)

        src_bbox = transform_bbox(eopatch.bbox, CRS.POP_WEB)
        src_transform = rasterio.transform.from_bounds(*src_bbox, width=width, height=height)

        dst_bbox = eopatch.bbox
        dst_transform = rasterio.transform.from_bounds(*dst_bbox, width=width, height=height)

        rasterio.warp.reproject(src_raster, dst_raster,
                                src_transform=src_transform, src_crs={'init': CRS.ogc_string(CRS.POP_WEB)},
                                src_nodata=0,
                                dst_transform=dst_transform, dst_crs={'init': CRS.ogc_string(eopatch.bbox.crs)},
                                dst_nodata=self.no_data_val)

        return dst_raster
686070
def raster_shpclip(r_fn, shp_fn, extent='raster', bbox=False, pad=None, invert=False, verbose=False):
    
    from pygeotools.lib import iolib
    from pygeotools.lib import warplib

    r_ds = iolib.fn_getds(r_fn)
    r_srs = get_ds_srs(r_ds)
    r_extent = ds_extent(r_ds)
    r_extent_geom = bbox2geom(r_extent)

    #NOTE: want to add spatial filter here to avoid reprojeting global RGI polygons, for example
    
    shp_ds = ogr.Open(shp_fn)
    lyr = shp_ds.GetLayer()
    shp_srs = lyr.GetSpatialRef()
    if not r_srs.IsSame(shp_srs):
        shp_ds = lyr_proj(lyr, r_srs)
        lyr = shp_ds.GetLayer()
    #This returns xmin, ymin, xmax, ymax
    shp_extent = lyr_extent(lyr)
    shp_extent_geom = bbox2geom(shp_extent)

    #Define the output - can set to either raster or shp
    #Could accept as cl arg
    out_srs = r_srs

    if extent == 'raster':
        out_extent = r_extent 
    elif extent == 'shp':
        out_extent = shp_extent
    elif extent == 'intersection':
        out_extent = geom_intersection([r_extent_geom, shp_extent_geom])
    elif extent == 'union':
        out_extent = geom_union([r_extent_geom, shp_extent_geom])
    else:
        print("Unexpected extent specification, reverting to input raster extent")
        out_extent = 'raster'

    #Add padding around shp_extent
    #Should implement buffer here
    if pad is not None:
        out_extent = pad_extent(out_extent, width=pad)

    print("Raster to clip: %s\nShapefile used to clip: %s" % (r_fn, shp_fn))
    if verbose:
        print(shp_extent) 
        print(r_extent)
        print(out_extent)

    r_ds = warplib.memwarp(r_ds, extent=out_extent, t_srs=out_srs, r='cubic')
    r = iolib.ds_getma(r_ds)

    #If bbox, return without clipping, otherwise, clip to polygons
    if not bbox:
        #Create binary mask from shp
        mask = shp2array(shp_fn, r_ds)
        if invert:
            mask = ~(mask)
        #Now apply the mask
        r = np.ma.array(r, mask=mask)
    #Return both the array and the dataset, needed for writing out
    #Should probably just write r to r_ds and return r_ds 
    return r, r_ds
731428
def read_data( self, step, filename = None ):
        """"""
        filename = get_default( filename, self.filename )

        out = {}

        fd = open( self.filename, 'r' )
        while 1:
            line = skip_read_line(fd, no_eof=True).split()
            if line[0] == 'POINT_DATA':
                break

        n_nod = int(line[1])

        while 1:
            line = skip_read_line(fd)
            if not line:
                break

            line = line.split()

            if line[0] == 'SCALARS':
                name, dtype, nc = line[1:]
                assert_(int(nc) == 1)
                fd.readline() # skip lookup table line
                
                data = nm.zeros((n_nod,), dtype=nm.float64)
                ii = 0
                while ii < n_nod:
                    data[ii] = float(fd.readline())
                    ii += 1

                out[name] = Struct( name = name,
                                    mode = 'vertex',
                                    data = data,
                                    dofs = None )

            elif line[0] == 'VECTORS':
                name, dtype = line[1:]
                data = []
                ii = 0
                while ii < n_nod:
                    data.append([float(val) for val in fd.readline().split()])
                    ii += 1

                out[name] = Struct( name = name,
                                    mode = 'vertex',
                                    data = nm.array(data, dtype=nm.float64),
                                    dofs = None )

            elif line[0] == 'CELL_DATA':
                break

            line = fd.readline()

        fd.close()

        return out
1157340
def shp2array(shp_fn, r_ds=None, res=None, extent=None, t_srs=None):
    
    shp_ds = ogr.Open(shp_fn)
    lyr = shp_ds.GetLayer()
    #This returns xmin, ymin, xmax, ymax
    shp_extent = lyr_extent(lyr)
    shp_srs = lyr.GetSpatialRef()
    # dst_dt = gdal.GDT_Byte
    ndv = 0
    if r_ds is not None:
        r_extent = ds_extent(r_ds)
        res = get_res(r_ds, square=True)[0] 
        if extent is None:
            extent = r_extent
        r_srs = get_ds_srs(r_ds)
        r_geom = ds_geom(r_ds)
        # dst_ns = r_ds.RasterXSize
        # dst_nl = r_ds.RasterYSize
        #Convert raster extent to shp_srs
        cT = osr.CoordinateTransformation(r_srs, shp_srs)
        r_geom_reproj = geom_dup(r_geom)
        r_geom_reproj.Transform(cT)
        r_geom_reproj.AssignSpatialReference(t_srs)
        lyr.SetSpatialFilter(r_geom_reproj)
        #lyr.SetSpatialFilter(ogr.CreateGeometryFromWkt(wkt))
    else:
        #TODO: clean this up
        if res is None:
            sys.exit("Must specify input res")
        if extent is None:
            print("Using input shp extent")
            extent = shp_extent
    if t_srs is None:
        t_srs = r_srs
    if not shp_srs.IsSame(t_srs):
        print("Input shp srs: %s" % shp_srs.ExportToProj4())
        print("Specified output srs: %s" % t_srs.ExportToProj4())
        out_ds = lyr_proj(lyr, t_srs)
        outlyr = out_ds.GetLayer()
    else:
        outlyr = lyr
    #outlyr.SetSpatialFilter(r_geom)
    m_ds = mem_ds(res, extent, srs=t_srs, dtype=gdal.GDT_Byte)
    b = m_ds.GetRasterBand(1)
    b.SetNoDataValue(ndv)
    gdal.RasterizeLayer(m_ds, [1], outlyr, burn_values=[1])
    a = b.ReadAsArray()
    a = ~(a.astype('Bool'))
    return a
866484
def cli(source_f, raster_f, output, verbose):
    
    with fiona.open(source_f, 'r') as source:
        source_driver = source.driver
        source_crs = source.crs
        sink_schema = source.schema.copy()

        source_geom = source.schema['geometry']
        if source_geom == 'Point':
            sink_schema['geometry'] = '3D Point'
        elif source_geom == 'LineString':
            sink_schema['geometry'] = '3D LineString'
        elif source_geom == '3D Point' or source_geom == '3D LineString':
            pass
        else:
            click.BadParameter("Source geometry type {} not implemented".format(source_geom))

        with rasterio.open(raster_f) as raster:
            if source_crs != raster.crs:
                click.BadParameter("Features and raster have different CRS.")
            if raster.count > 1:
                warnings.warn("Found {0} bands in {1}, expected a single band raster".format(raster.bands, raster_f))
            supported = ['int16', 'int32', 'float32', 'float64']
            if raster.dtypes[0] not in supported:
                warnings.warn("Found {0} type in {1}, expected one of {2}".format(raster.dtypes[0]), raster_f, supported)
            with fiona.open(
                output, 'w',
                driver=source_driver,
                crs=source_crs,
                schema=sink_schema) as sink:

                for feature in source:
                    try:
                        feature_z = drapery.drape(raster, feature)
                        sink.write({
                            'geometry': mapping(feature_z),
                            'properties': feature['properties'],
                        })
                    except Exception:
                        logging.exception("Error processing feature %s:", feature['id'])
863096
def _reads_in_peaks(bam_file, peaks_file, sample):
    """"""
    if not peaks_file:
        return {}
    rip = number_of_mapped_reads(sample, bam_file, bed_file = peaks_file)
    return {"metrics": {"RiP": rip}}
763679
def _read(self, directory, filename, session, path, name, extension, spatial, spatialReferenceID, replaceParamFile):
        
        # Assign file extension attribute to file object
        self.fileExtension = extension

        # Open file and parse into a data structure
        with open(path, 'r') as f:
            for line in f:
                sline = line.strip().split()

                if len(sline) == 1:
                    self.numLocations = sline[0]
                else:
                    # Create GSSHAPY OutputLocation object
                    location = OutputLocation(linkOrCellI=sline[0],
                                              nodeOrCellJ=sline[1])

                    # Associate OutputLocation with OutputLocationFile
                    location.outputLocationFile = self