Absolute elevation difference range filter using values from a source raster file and a reference raster file
326452
def dz_fltr(dem_fn, refdem_fn, perc=None, rangelim=(0,30), smooth=False):
    
    try:
        open(refdem_fn)
    except IOError:
        sys.exit('Unable to open reference DEM: %s' % refdem_fn)

    from pygeotools.lib import warplib
    dem_ds, refdem_ds = warplib.memwarp_multi_fn([dem_fn, refdem_fn], res='first', extent='first', t_srs='first')
    dem = iolib.ds_getma(dem_ds)
    refdem = iolib.ds_getma(refdem_ds)
    out = dz_fltr_ma(dem, refdem, perc, rangelim, smooth)
    return out
1259465
def rd2xy(input,ra=None,dec=None,coordfile=None,colnames=None,
            precision=6,output=None,verbose=True):
    
    single_coord = False
    if coordfile is not None:
        if colnames in blank_list:
            colnames = ['c1','c2']
        elif isinstance(colnames,type('a')):
            colnames = colnames.split(',')
        # convert input file coordinates to lists of decimal degrees values
        xlist,ylist = tweakutils.readcols(coordfile,cols=colnames)
    else:
        if isinstance(ra,np.ndarray):
            ralist = ra.tolist()
            declist = dec.tolist()
        elif not isinstance(ra, list):
            ralist = [ra]
            declist = [dec]
        else:
            ralist = ra
            declist = dec
        xlist  = [0]*len(ralist)
        ylist = [0]*len(ralist)
        if len(xlist) == 1:
            single_coord = True
        for i,(r,d) in enumerate(zip(ralist,declist)):
            # convert input value into decimal degrees value
            xval,yval = tweakutils.parse_skypos(r,d)
            xlist[i] = xval
            ylist[i] = yval

    # start by reading in WCS+distortion info for input image
    inwcs = wcsutil.HSTWCS(input)
    if inwcs.wcs.is_unity():
        print("####\nNo valid WCS found in {}.\n  Results may be invalid.\n####\n".format(input))

    # Now, convert pixel coordinates into sky coordinates
    try:
        outx,outy = inwcs.all_world2pix(xlist,ylist,1)
    except RuntimeError:
        outx,outy = inwcs.wcs_world2pix(xlist,ylist,1)

    # add formatting based on precision here...
    xstr = []
    ystr = []
    fmt = "%."+repr(precision)+"f"
    for x,y in zip(outx,outy):
        xstr.append(fmt%x)
        ystr.append(fmt%y)

    if verbose or (not verbose and util.is_blank(output)):
        print ('# Coordinate transformations for ',input)
        print('# X      Y         RA             Dec\n')
        for x,y,r,d in zip(xstr,ystr,xlist,ylist):
            print("%s  %s    %s  %s"%(x,y,r,d))

    # Create output file, if specified
    if output:
        f = open(output,mode='w')
        f.write("# Coordinates converted from %s\n"%input)
        for x,y in zip(xstr,ystr):
            f.write('%s    %s\n'%(x,y))
        f.close()
        print('Wrote out results to: ',output)

    if single_coord:
        outx = outx[0]
        outy = outy[0]

    return outx, outy
715214
def pixel_difference_filter(pixel, x, y)
      chans = []
      [:r, :b, :g].each do |chan| 
        chans << channel_difference(chan, pixel, x, y)
      end
      @images.last[x,y] = ChunkyPNG::Color.rgb(chans[0], chans[1], chans[2])
    end
237378
def mask_raster(in_raster, mask, out_raster):
        
        if is_string(in_raster) and is_string(out_raster):
            in_raster = [str(in_raster)]
            out_raster = [str(out_raster)]
        if len(in_raster) != len(out_raster):
            raise RuntimeError('input raster and output raster must have the same size.')

        maskr = RasterUtilClass.read_raster(mask)
        rows = maskr.nRows
        cols = maskr.nCols
        maskdata = maskr.data
        temp = maskdata == maskr.noDataValue
        for inr, outr in zip(in_raster, out_raster):
            origin = RasterUtilClass.read_raster(inr)
            if origin.nRows == rows and origin.nCols == cols:
                masked = numpy.where(temp, origin.noDataValue, origin.data)
            else:
                masked = numpy.ones((rows, cols)) * origin.noDataValue
                # TODO, the following loop should be optimized by numpy or numba
                for i in range(rows):
                    for j in range(cols):
                        if maskdata[i][j] == maskr.noDataValue:
                            continue
                        # get the center point coordinate of current cell
                        tempx, tempy = maskr.get_central_coors(i, j)
                        tempv = origin.get_value_by_xy(tempx, tempy)
                        if tempv is None:
                            continue
                        masked[i][j] = tempv
            RasterUtilClass.write_gtiff_file(outr, maskr.nRows, maskr.nCols, masked,
                                             maskr.geotrans, maskr.srs,
                                             origin.noDataValue, origin.dataType)
87450
function adjustRange(w, bisect) {
  var r0 = w.i0,
      r1 = w.i1 - 1,
      c = w.compare,
      d = w.data,
      n = d.length - 1;

  if (r0 > 0 && !c(d[r0], d[r0-1])) w.i0 = bisect.left(d, d[r0]);
  if (r1 < n && !c(d[r1], d[r1+1])) w.i1 = bisect.right(d, d[r1]);
}
1026964
def raster_statistics(raster_file):
        
        ds = gdal_Open(raster_file)
        band = ds.GetRasterBand(1)
        minv, maxv, meanv, std = band.ComputeStatistics(False)
        return minv, maxv, meanv, std
428057
def get_ras(repi, theta, mag, coeff):
    
    rx = 100.
    ras = 200.
    #
    # calculate the difference between epicentral distances
    dff = fnc(ras, repi, theta, mag, coeff)
    while abs(dff) > 1e-3:
        # update the value of distance computed
        if dff > 0.:
            ras = ras - rx
        else:
            ras = ras + rx
        dff = fnc(ras, repi, theta, mag, coeff)
        rx = rx / 2.
        if rx < 1e-3:
            break
    return ras
467770
def set_elevation(self, elevation_grid_path, mask_shapefile):
        
        # ADD ELEVATION FILE
        ele_file = ElevationGridFile(project_file=self.project_manager,
                                     session=self.db_session)
        ele_file.generateFromRaster(elevation_grid_path,
                                    mask_shapefile,
                                    load_raster_to_db=self.load_rasters_to_db)
591956
def max_diff(dset1,dset2):
    ''''''
    for dset in [dset1,dset2]:
        if not os.path.exists(dset):
            nl.notify('Error: Could not find file: %s' % dset,level=nl.level.error)
            return float('inf')
    try:
        dset1_d = nib.load(dset1)
        dset2_d = nib.load(dset2)
        dset1_data = dset1_d.get_data()
        dset2_data = dset2_d.get_data()
    except IOError:
        nl.notify('Error: Could not read files %s and %s' % (dset1,dset2),level=nl.level.error)
        return float('inf')
    try:
        old_err = np.seterr(divide='ignore',invalid='ignore')
        max_val = 100*np.max(np.ma.masked_invalid(np.double(dset1_data - dset2_data) / ((dset1_data+dset2_data)/2)))
        np.seterr(**old_err)
        return max_val
    except ValueError:
        return float('inf')
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)
450673
def compare(what, imt, calc_ids, files, samplesites=100, rtol=.1, atol=1E-4):
    
    sids, imtls, poes, arrays = getdata(what, calc_ids, samplesites)
    try:
        levels = imtls[imt]
    except KeyError:
        sys.exit(
            '%s not found. The available IMTs are %s' % (imt, list(imtls)))
    imt2idx = {imt: i for i, imt in enumerate(imtls)}
    head = ['site_id'] if files else ['site_id', 'calc_id']
    if what == 'hcurves':
        array_imt = arrays[:, :, imtls(imt)]
        header = head + ['%.5f' % lvl for lvl in levels]
    else:  # hmaps
        array_imt = arrays[:, :, imt2idx[imt]]
        header = head + [str(poe) for poe in poes]
    rows = collections.defaultdict(list)
    diff_idxs = get_diff_idxs(array_imt, rtol, atol)
    if len(diff_idxs) == 0:
        print('There are no differences within the tolerance of %d%%' %
              (rtol * 100))
        return
    arr = array_imt.transpose(1, 0, 2)  # shape (N, C, L)
    for sid, array in sorted(zip(sids[diff_idxs], arr[diff_idxs])):
        for calc_id, cols in zip(calc_ids, array):
            if files:
                rows[calc_id].append([sid] + list(cols))
            else:
                rows['all'].append([sid, calc_id] + list(cols))
    if files:
        fdict = {calc_id: open('%s.txt' % calc_id, 'w')
                 for calc_id in calc_ids}
        for calc_id, f in fdict.items():
            f.write(views.rst_table(rows[calc_id], header))
            print('Generated %s' % f.name)
    else:
        print(views.rst_table(rows['all'], header))
371803
private BufferedImage resample(BufferedImage pSource, BufferedImage pDest, InterpolationFilter pFilter) {
        final int dstWidth = pDest.getWidth();
        final int dstHeight = pDest.getHeight();

        final int srcWidth = pSource.getWidth();
        final int srcHeight = pSource.getHeight();

        /* create intermediate column to hold horizontal dst column zoom */
        final ColorModel cm = pSource.getColorModel();
//        final WritableRaster work = cm.createCompatibleWritableRaster(1, srcHeight);
        final WritableRaster work = ImageUtil.createCompatibleWritableRaster(pSource, cm, 1, srcHeight);

        double xscale = (double) dstWidth / (double) srcWidth;
        double yscale = (double) dstHeight / (double) srcHeight;

        ContributorList[] contribY = new ContributorList[dstHeight];
        for (int i = 0; i < contribY.length; i++) {
            contribY[i] = new ContributorList();
        }

        // TODO: What to do when fwidth > srcHeight or dstHeight
        double fwidth = pFilter.support();
        if (yscale < 1.0) {
            double width = fwidth / yscale;
            double fscale = 1.0 / yscale;

            if (width <= .5) {
                // Reduce to point sampling.
                width = .5 + 1.0e-6;
                fscale = 1.0;
            }

            for (int i = 0; i < dstHeight; i++) {
                //contribY[i].n = 0;
                contribY[i].p = new Contributor[(int) (width * 2.0 + 1 + 0.5)];

                double center = (double) i / yscale;
                int left = (int) Math.ceil(center - width);
                int right = (int) Math.floor(center + width);

                double density = 0.0;

                for (int j = left; j <= right; j++) {
                    double weight = center - (double) j;
                    weight = pFilter.filter(weight / fscale) / fscale;
                    int n;
                    if (j < 0) {
                        n = -j;
                    }
                    else if (j >= srcHeight) {
                        n = (srcHeight - j) + srcHeight - 1;
                    }
                    else {
                        n = j;
                    }

                    /**/
                    if (n >= srcHeight) {
                        n = n % srcHeight;
                    }
                    else if (n < 0) {
                        n = srcHeight - 1;
                    }
                    /**/

                    int k = contribY[i].n++;
                    contribY[i].p[k] = new Contributor();
                    contribY[i].p[k].pixel = n;
                    contribY[i].p[k].weight = weight;

                    density += weight;
                }

                if ((density != 0.0) && (density != 1.0)) {
                    //Normalize.
                    density = 1.0 / density;
                    for (int k = 0; k < contribY[i].n; k++) {
                        contribY[i].p[k].weight *= density;
                    }
                }
            }
        }
        else {
            for (int i = 0; i < dstHeight; ++i) {
                //contribY[i].n = 0;
                contribY[i].p = new Contributor[(int) (fwidth * 2 + 1 + 0.5)];

                double center = (double) i / yscale;
                double left = Math.ceil(center - fwidth);
                double right = Math.floor(center + fwidth);
                for (int j = (int) left; j <= right; ++j) {
                    double weight = center - (double) j;
                    weight = pFilter.filter(weight);
                    int n;
                    if (j < 0) {
                        n = -j;
                    }
                    else if (j >= srcHeight) {
                        n = (srcHeight - j) + srcHeight - 1;
                    }
                    else {
                        n = j;
                    }

                    /**/
                    if (n >= srcHeight) {
                        n = n % srcHeight;
                    }
                    else if (n < 0) {
                        n = srcHeight - 1;
                    }
                    /**/

                    int k = contribY[i].n++;
                    contribY[i].p[k] = new Contributor();
                    contribY[i].p[k].pixel = n;
                    contribY[i].p[k].weight = weight;
                }
            }
        }

        final Raster raster = pSource.getRaster();
        final WritableRaster out = pDest.getRaster();

        // TODO: This is not optimal for non-byte-packed rasters...
        // (What? Maybe I implemented the fix, but forgot to remove the TODO?)
        final int numChannels = raster.getNumBands();
        final int[] channelMax = new int[numChannels];
        for (int k = 0; k < numChannels; k++) {
            channelMax[k] = (1 << pSource.getColorModel().getComponentSize(k)) - 1;
        }

        for (int xx = 0; xx < dstWidth; xx++) {
            ContributorList contribX = calcXContrib(xscale, fwidth, srcWidth, pFilter, xx);
            /* Apply horiz filter to make dst column in tmp. */
            for (int k = 0; k < srcHeight; k++) {
                for (int channel = 0; channel < numChannels; channel++) {

                    double weight = 0.0;
                    boolean bPelDelta = false;
                    // TODO: This line throws index out of bounds, if the image
                    // is smaller than filter.support()
                    double pel = raster.getSample(contribX.p[0].pixel, k, channel);
                    for (int j = 0; j < contribX.n; j++) {
                        double pel2 = j == 0 ? pel : raster.getSample(contribX.p[j].pixel, k, channel);
                        if (pel2 != pel) {
                            bPelDelta = true;
                        }
                        weight += pel2 * contribX.p[j].weight;
                    }
                    weight = bPelDelta ? round(weight) : pel;

                    if (weight < 0) {
                        weight = 0;
                    }
                    else if (weight > channelMax[channel]) {
                        weight = channelMax[channel];
                    }

                    work.setSample(0, k, channel, weight);

                }
            }/* next row in temp column */

            /* The temp column has been built. Now stretch it vertically into dst column. */
            for (int i = 0; i < dstHeight; i++) {
                for (int channel = 0; channel < numChannels; channel++) {

                    double weight = 0.0;
                    boolean bPelDelta = false;
                    double pel = work.getSample(0, contribY[i].p[0].pixel, channel);

                    for (int j = 0; j < contribY[i].n; j++) {
                        // TODO: This line throws index out of bounds, if the image
                        // is smaller than filter.support()
                        double pel2 = j == 0 ? pel : work.getSample(0, contribY[i].p[j].pixel, channel);
                        if (pel2 != pel) {
                            bPelDelta = true;
                        }
                        weight += pel2 * contribY[i].p[j].weight;
                    }
                    weight = bPelDelta ? round(weight) : pel;
                    if (weight < 0) {
                        weight = 0;
                    }
                    else if (weight > channelMax[channel]) {
                        weight = channelMax[channel];
                    }

                    out.setSample(xx, i, channel, weight);
                }
            }/* next dst row */
        }/* next dst column */
        return pDest;
    }
517870
def median_fltr(dem, fsize=7, origmask=False):
    
    print("Applying median filter with size %s" % fsize)
    from scipy.ndimage.filters import median_filter
    dem_filt_med = median_filter(dem.filled(np.nan), fsize)
    #Now mask all nans
    out = np.ma.fix_invalid(dem_filt_med, copy=False, fill_value=dem.fill_value)
    if origmask:
        out = np.ma.array(out, mask=dem.mask, fill_value=dem.fill_value)
    out.set_fill_value(dem.fill_value)
    return out
417706
def lowPassFilter(self, threshold):
        '''
        
        '''
        if not threshold:
            return
        rows, cols = self.img.shape
        tx = int(cols * threshold * 0.25)
        ty = int(rows * threshold * 0.25)
        # upper side
        self.fshift[rows - tx:rows, :] = 0
        # lower side
        self.fshift[0:tx, :] = 0
        # left side
        self.fshift[:, 0:ty] = 0
        # right side
        self.fshift[:, cols - ty:cols] = 0
910003
def _residual(x, in_filter, out_filter, stride,
              activate_before_residual=False):
  """"""
  if activate_before_residual:
    with tf.variable_scope('shared_activation'):
      x = _batch_norm('init_bn', x)
      x = _relu(x, 0.1)
      orig_x = x
  else:
    with tf.variable_scope('residual_only_activation'):
      orig_x = x
      x = _batch_norm('init_bn', x)
      x = _relu(x, 0.1)

  with tf.variable_scope('sub1'):
    x = _conv('conv1', x, 3, in_filter, out_filter, stride)

  with tf.variable_scope('sub2'):
    x = _batch_norm('bn2', x)
    x = _relu(x, 0.1)
    x = _conv('conv2', x, 3, out_filter, out_filter, [1, 1, 1, 1])

  with tf.variable_scope('sub_add'):
    if in_filter != out_filter:
      orig_x = tf.nn.avg_pool(orig_x, stride, stride, 'VALID')
      orig_x = tf.pad(
          orig_x, [[0, 0], [0, 0],
                   [0, 0], [(out_filter - in_filter) // 2,
                            (out_filter - in_filter) // 2]])
    x += orig_x

  tf.logging.debug('image after unit %s', x.get_shape())
  return x
1331348
def from_file(filename, **kwargs):
    
    ndv, xsize, ysize, geot, projection, datatype = get_geo_info(filename, **kwargs)
    data = gdalnumeric.LoadFile(filename, **kwargs)
    data = np.ma.masked_array(data, mask=data == ndv, fill_value=ndv)
    return GeoRaster(data, geot, nodata_value=ndv, projection=projection, datatype=datatype)
1070527
def getRastersAsPngs(self, session, tableName, rasterIds, postGisRampString, rasterField='raster', rasterIdField='id',  cellSize=None, resampleMethod='NearestNeighbour'):
        
        # Validate
        VALID_RESAMPLE_METHODS = ('NearestNeighbour', 'Bilinear', 'Cubic', 'CubicSpline', 'Lanczos')

        # Validate
        if resampleMethod not in VALID_RESAMPLE_METHODS:
            print('RASTER CONVERSION ERROR: {0} is not a valid resampleMethod.'
                  ' Please use either {1}'.format(resampleMethod,
                                                  ', '.join(VALID_RESAMPLE_METHODS)))

        if cellSize is not None:
            if not self.isNumber(cellSize):
                raise ValueError('RASTER CONVERSION ERROR: cellSize must be a number or None.')

        # Convert raster ids into formatted string
        rasterIdsString = '({0})'.format(', '.join(rasterIds))

        if cellSize is not None:
            statement = '''
                        SELECT ST_AsPNG(ST_Transform(ST_ColorMap(ST_Rescale({0}, {5}, '{6}'), 1, '{4}'), 4326, 'Bilinear')) As png
                        FROM {1}
                        WHERE {2} IN {3};
                        '''.format(rasterField, tableName, rasterIdField, rasterIdsString, postGisRampString, cellSize,
                                   resampleMethod)
        else:
            statement = '''
                        SELECT ST_AsPNG(ST_Transform(ST_ColorMap({0}, 1, '{4}'), 4326, 'Bilinear')) As png
                        FROM {1}
                        WHERE {2} IN {3};
                        '''.format(rasterField, tableName, rasterIdField, rasterIdsString, postGisRampString)
        result = session.execute(statement)
        return result