Reprojects the raster data from Geopedia's CRS (POP_WEB) to EOPatch's CRS.
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
630893
def interp_value(self, lat, lon, indexed=False):
        
        (px, py) = self.grid_coordinates.projection_to_raster_coords(lat, lon)
        if indexed:
            return self.raster_data[round(py), round(px)]
        else:
#             from scipy.interpolate import interp2d
#             f_interp = interp2d(self.grid_coordinates.x_axis, self.grid_coordinates.y_axis, self.raster_data, bounds_error=True)
#             return f_interp(lon, lat)[0]
            from scipy.ndimage import map_coordinates
            ret = map_coordinates(self.raster_data, [[py], [px]], order=1)  # linear interp
            return ret[0]
198924
def make_geo(attrs_dict, globe):
    """"""
    attr_mapping = [('satellite_height', 'perspective_point_height'),
                    ('sweep_axis', 'sweep_angle_axis')]
    kwargs = CFProjection.build_projection_kwargs(attrs_dict, attr_mapping)

    # CartoPy can't handle central latitude for Geostationary (nor should it)
    # Just remove it if it's 0.
    if not kwargs.get('central_latitude'):
        kwargs.pop('central_latitude', None)

    # If sweep_angle_axis is not present, we should look for fixed_angle_axis and adjust
    if 'sweep_axis' not in kwargs:
        kwargs['sweep_axis'] = 'x' if attrs_dict['fixed_angle_axis'] == 'y' else 'y'

    return ccrs.Geostationary(globe=globe, **kwargs)
1152488
def to_(self, data_pts):
        """"""
        data_pts = np.asarray(data_pts, dtype=np.float)
        has_z = (data_pts.shape[-1] > 2)

        if self.use_center:
            data_pts = data_pts - self.viewer.data_off

        # subtract data indexes at center reference pixel
        ref_pt = [self.viewer._org_x, self.viewer._org_y]
        if has_z:
            ref_pt.append(self.viewer._org_z)

        off_pts = np.subtract(data_pts, ref_pt)
        return off_pts
207949
def to_pandas(raster, name='value'):
    
    df = pd.DataFrame(raster.raster)
    df = df.stack()
    df = df.reset_index()
    df.columns = ['row', 'col', name]
    df['x'] = df.col.apply(lambda col: raster.geot[0]+(col)*raster.geot[1])
    df['y'] = df.row.apply(lambda row: raster.geot[3]+(row)*raster.geot[-1])
    return df
890642
def proj4_to_wkt(projection):
    """"""

    srs = osgeo.osr.SpatialReference()
    srs.ImportFromProj4(projection.srs)

    return srs.ExportToWkt()
900841
def raster_to_asc(raster_f, asc_f):
        
        raster_r = RasterUtilClass.read_raster(raster_f)
        RasterUtilClass.write_asc_file(asc_f, raster_r.data, raster_r.nCols, raster_r.nRows,
                                       raster_r.geotrans, raster_r.noDataValue)
318841
def from_record(cls, record, crs, schema=None):
        """"""
        properties = cls._to_properties(record, schema)
        vector = GeoVector(shape(record['geometry']), crs)
        if record.get('raster'):
            assets = {k: dict(type=RASTER_TYPE, product='visual', **v) for k, v in record.get('raster').items()}
        else:
            assets = record.get('assets', {})
        return cls(vector, properties, assets)
384502
def to_sky(self, wcs, origin=_DEFAULT_WCS_ORIGIN, mode=_DEFAULT_WCS_MODE):
        
        return SkyCoord.from_pixel(
            xp=self.x, yp=self.y, wcs=wcs,
            origin=origin, mode=mode,
        )
21647
def _get_tile_when_web_mercator_crs(self, x_tile, y_tile, zoom,
                                        bands=None, masked=None,
                                        resampling=Resampling.cubic):
        
        roi = GeoVector.from_xyz(x_tile, y_tile, zoom)
        coordinates = roi.get_bounds(WEB_MERCATOR_CRS)
        window = self._window(coordinates, to_round=False)
        bands = bands or list(range(1, self.num_bands + 1))
        # we know the affine the result should produce becuase we know where
        # it is located by the xyz, therefore we calculate it here
        ratio = MERCATOR_RESOLUTION_MAPPING[zoom] / self.resolution()

        # the affine should be calculated before rounding the window values
        affine = self.window_transform(window)
        affine = affine * Affine.scale(ratio, ratio)

        window = Window(round(window.col_off),
                        round(window.row_off),
                        round(window.width),
                        round(window.height))
        return self.get_window(window, bands=bands, xsize=256, ysize=256, masked=masked, affine=affine)
258837
def offset_to_sky(skydir, offset_lon, offset_lat,
                  coordsys='CEL', projection='AIT'):
    """"""

    offset_lon = np.array(offset_lon, ndmin=1)
    offset_lat = np.array(offset_lat, ndmin=1)

    w = create_wcs(skydir, coordsys, projection)
    pixcrd = np.vstack((offset_lon, offset_lat)).T

    return w.wcs_pix2world(pixcrd, 0)
366614
def reconstruct(self, D=None, X=None):
        """"""

        if D is None:
            D = self.getdict(crop=False)
        if X is None:
            X = self.getcoef()
        Df = sl.rfftn(D, self.xstep.cri.Nv, self.xstep.cri.axisN)
        Xf = sl.rfftn(X, self.xstep.cri.Nv, self.xstep.cri.axisN)
        DXf = sl.inner(Df, Xf, axis=self.xstep.cri.axisM)
        return sl.irfftn(DXf, self.xstep.cri.Nv, self.xstep.cri.axisN)
1246027
function () {
        var localTransform = this.localTransform.array;
        var worldTransform = this.worldTransform.array;
        if (this._parent) {
            mat4.multiplyAffine(
                worldTransform,
                this._parent.worldTransform.array,
                localTransform
            );
        }
        else {
            mat4.copy(worldTransform, localTransform);
        }
    }
800990
def _update_raster_info(self, **update_props):
        

        tree_to_update = update_props['tree_to_update']
        prop = update_props['prop']
        values = update_props.pop('values')

        # Update number of dimensions at raster_info root (applies to all dimensions below)

        xroot, xpath = None, self._data_map['_ri_num_dims']
        raster_info = [update_property(tree_to_update, xroot, xpath, prop, values.get('dimensions', u''))]

        # Derive vertical, longitude, and latitude dimensions from raster_info

        xpath_root = self._get_xroot_for(prop)
        xpath_map = self._data_structures[prop]

        v_dimension = {}
        if values.get('vertical_count'):
            v_dimension = v_dimension.fromkeys(xpath_map, u'')
            v_dimension['type'] = 'vertical'
            v_dimension['size'] = values.get('vertical_count', u'')

        x_dimension = {}
        if values.get('column_count') or values.get('x_resolution'):
            x_dimension = x_dimension.fromkeys(xpath_map, u'')
            x_dimension['type'] = 'column'
            x_dimension['size'] = values.get('column_count', u'')
            x_dimension['value'] = values.get('x_resolution', u'')

        y_dimension = {}
        if values.get('row_count') or values.get('y_resolution'):
            y_dimension = y_dimension.fromkeys(xpath_map, u'')
            y_dimension['type'] = 'row'
            y_dimension['size'] = values.get('row_count', u'')
            y_dimension['value'] = values.get('y_resolution', u'')

        # Update derived dimensions as complex list, and append affected elements for return

        update_props['prop'] = RASTER_DIMS
        update_props['values'] = [v_dimension, x_dimension, y_dimension]

        raster_info += update_complex_list(xpath_root=xpath_root, xpath_map=xpath_map, **update_props)

        return raster_info
140120
def transform(data, target_wd, target_ht, is_train, box):
    """"""
    if box is not None:
        x, y, w, h = box
        data = data[y:min(y+h, data.shape[0]), x:min(x+w, data.shape[1])]

    # Resize to target_wd * target_ht.
    data = mx.image.imresize(data, target_wd, target_ht)

    # Normalize in the same way as the pre-trained model.
    data = data.astype(np.float32) / 255.0
    data = (data - mx.nd.array([0.485, 0.456, 0.406])) / mx.nd.array([0.229, 0.224, 0.225])

    if is_train:
        if random.random() < 0.5:
            data = nd.flip(data, axis=1)
        data, _ = mx.image.random_crop(data, (224, 224))
    else:
        data, _ = mx.image.center_crop(data, (224, 224))

    # Transpose from (target_wd, target_ht, 3)
    # to (3, target_wd, target_ht).
    data = nd.transpose(data, (2, 0, 1))

    # If image is greyscale, repeat 3 times to get RGB image.
    if data.shape[0] == 1:
        data = nd.tile(data, (3, 1, 1))
    return data.reshape((1,) + data.shape)
1282010
def get_dem(bounds, out_file="dem.tif", src_crs="EPSG:3005", dst_crs="EPSG:3005", resolution=25):
    
    bbox = ",".join([str(b) for b in bounds])
    # todo: validate resolution units are equivalent to src_crs units
    # build request
    payload = {
        "service": "WCS",
        "version": "1.0.0",
        "request": "GetCoverage",
        "coverage": "pub:bc_elevation_25m_bcalb",
        "Format": "GeoTIFF",
        "bbox": bbox,
        "CRS": src_crs,
        "RESPONSE_CRS": dst_crs,
        "resx": str(resolution),
        "resy": str(resolution),
    }
    # request data from WCS
    r = requests.get(bcdata.WCS_URL, params=payload)
    # save to tiff
    if r.status_code == 200:
        with open(out_file, "wb") as file:
            file.write(r.content)
        return out_file
    else:
        raise RuntimeError(
            "WCS request failed with status code {}".format(str(r.status_code))
        )
495165
def execute(self, eopatch):
        
        bbox_map = self._get_submap(eopatch)
        height, width = self._get_shape(eopatch)
        data_transform = rasterio.transform.from_bounds(*eopatch.bbox, width=width, height=height)

        if self.feature_name in eopatch[self.feature_type]:
            raster = eopatch[self.feature_type][self.feature_name].squeeze()
        else:
            raster = np.ones((height, width), dtype=self.raster_dtype) * self.no_data_value

        if not bbox_map.empty:
            rasterio.features.rasterize([(bbox_map.cascaded_union.buffer(0), self.raster_value)], out=raster,
                                        transform=data_transform, dtype=self.raster_dtype)

        eopatch[self.feature_type][self.feature_name] = raster[..., np.newaxis]

        return eopatch