Python get elevation

def get_elevation(lat, lng):
    """
    Get the elevation at point (lat,lng) using the currently opened interface
    :param lat: 
    :param lng: 
    :return:
    """
    try:
        elevation = interface.lookup(lat, lng)
    except:
        return {
            'latitude': lat,
            'longitude': lng,
            'error': 'No such coordinate (%s, %s)' % (lat, lng)
        }

    return {
        'latitude': lat,
        'longitude': lng,
        'elevation': elevation
    } 

def get_min_max_elevation(gdir):
    """Reads the DEM and computes the minimal and maximal glacier surface
     elevation in meters asl, from the given (RGI) glacier outline.

    Parameters
    ----------
    gdir : :py:class:`oggm.GlacierDirectory`

    Returns
    -------
    [float, float]
        minimal and maximal glacier surface elevation [m asl.]

    """
    # open DEM file and mask the glacier surface area
    fpath = gdir.get_filepath('gridded_data')
    with ncDataset(fpath) as nc:
        mask = nc.variables['glacier_mask'][:]
        topo = nc.variables['topo'][:]
    # get relevant elevation information
    min_elev = np.min(topo[np.where(mask == 1)])
    max_elev = np.max(topo[np.where(mask == 1)])

    return min_elev, max_elev 

def get_elevation(self, sweep, copy=False):
        """
        Return an array of elevation angles for a given sweep.

        Parameters
        ----------
        sweep : int
            Sweep number to retrieve data for, 0 based.
        copy : bool, optional
            True to return a copy of the elevations. False, the default,
            returns a view of the elevations (when possible), changing this
            data will change the data in the underlying Radar object.

        Returns
        -------
        azimuths : array
            Array containing the elevation angles for a given sweep.

        """
        s = self.get_slice(sweep)
        elevation = self.elevation['data'][s]
        if copy:
            return elevation.copy()
        else:
            return elevation 

def getElevation(rX, rY, rZ, coords, interp="linear"):
    """
    This function interpolates elevation from the regular grid to the triamgular mesh
    using **SciPy** *interpn* funciton.

    Args:
        rX: numpy arrays containing the X coordinates from the regular grid.
        rY: numpy arrays containing the Y coordinates from the regular grid.
        rZ: numpy arrays containing the Z coordinates from the regular grid.
        coords: numpy float-type array containing X, Y coordinates for the TIN nodes.
        interp: interpolation method as in *SciPy interpn function* (default: 'linear')

    Returns:
        - elev - numpy array containing the updated elevations for the local domain.
    """

    # Set new elevation to 0
    elev = numpy.zeros(len(coords[:, 0]))

    # Get the TIN points elevation values using the regular grid dataset
    elev = interpn((rX, rY), rZ, (coords[:, :2]), method=interp)

    return elev 

def get_reference_elevation(input, recGrid, elevation):
    """
    The following function define the elevation from the TIN to a regular grid...

    Args:
        input: class containing XML input file parameters.
        recGrid: class describing the regular grid characteristics.
        elevation: TIN elevation mesh.

    Returns:
        - ref_elev - interpolated elevation on the regular grid
    """
    if input.searef:
        x_ref, y_ref = input.searef
        pts = recGrid.tinMesh["vertices"]
        ref_elev = griddata(
            points=pts, values=elevation, xi=[x_ref, y_ref], method="nearest"
        )
    else:
        ref_elev = 0.0

    return ref_elev 

def get_elevation_value(label):
    """
    Returns the angle in degrees represented by a elevation label int.

    Parameters
    ----------
    label: int
      Elevation label.
    """

    name = 'elevation'
    _check_is_integral(name, label)
    _check_range(name, label, -1, 8)

    if label == -1:
        return None
    else:
        return label * 5 + 30 

def get_elevation_geoid(self, lat, lon, block_size=50000):
        """
        Get the elevation value relative to the geoid.

        Parameters
        ----------
        lat : numpy.ndarray|list|tuple|int|float
        lon : numpy.ndarray|list|tuple|int|float
        block_size : int|None
            If `None`, then the entire calculation will proceed as a single block.
            Otherwise, block processing using blocks of the given size will be used.
            The minimum value used for this is 50,000, and any smaller value will be
            replaced with 50,000. Default is 50,000.

        Returns
        -------
        numpy.ndarray
            the elevation relative to the geoid
        """

        raise NotImplementedError 

def get_elevation_hae(self, lat, lon, block_size=50000):
        """
        Get the elevation value relative to the WGS-84 ellipsoid.

        Parameters
        ----------
        lat : numpy.ndarray|list|tuple|int|float
        lon : numpy.ndarray|list|tuple|int|float
        block_size : int|None
            If `None`, then the entire calculation will proceed as a single block.
            Otherwise, block processing using blocks of the given size will be used.
            The minimum value used for this is 50,000, and any smaller value will be
            replaced with 50,000. Default is 50,000.

        Returns
        -------
        numpy.ndarray
            the elevation relative to the WGS-84 ellipsoid.
        """

        raise NotImplementedError 

def get_elevation(lat, lon):
    hgt_file = get_file_name(lat, lon)
    if hgt_file:
        return read_elevation_from_file(hgt_file, lat, lon)
    # Treat it as data void as in SRTM documentation
    # if file is absent
    return -32768 

def getTokenElevationType(handle=-1):

    token = HANDLE(0)
    etype = DWORD(0)
    outsize = DWORD(0)
    if not advapi32.OpenProcessToken(handle, TOKEN_QUERY, addressof(token)):
        raise Exception('Invalid Process Handle: %d' % handle)

    advapi32.GetTokenInformation(token, TokenElevationType, addressof(etype), 4, addressof(outsize))

    return etype.value 

def getSlantRangeElevation(self, groundRange, z):
        """Convert ground range (great circle distance) and height above
        the earth's surface to slant range (along the beam) and elevation angle.
        Follows Doviak and Zrnic 1993, eq. 2.28"""

        lat = self.ctrLat * pi / 180.0

        #figure out earth's radius at radar's lat ... non-spherical earth model
        e2 = self.eccen           # First eccentricity squared - WGS-84 value = 0.00669437999013
        a = self.Requator         # Equatorial radius  - WGS-84 value = 6378137.0
        Rearth = a/sqrt(1-e2*(sin(lat))**2) # radius of curvature

        Rprime = self.effectiveRadiusMultiplier * Rearth

        h = array(z - self.ctrAlt, dtype='float64')
        s = array(groundRange, dtype='float64')

        # Use law of cosines (Side-Angle-Side triangle theorem) with
        # R', R'+h as sides and s/R' as the angle to get slant range
        r  = sqrt(Rprime**2.0 + (Rprime+h)**2.0 - 2*(Rprime+h)*Rprime*cos(s/Rprime))
        # Inverse of eq. 2.28c in Doviak and Zrnic 1993
        # Will return NaN for r=0, and only positive angles
        el = atleast_1d(arccos((Rprime+h) * sin(s/Rprime) / r))
        # Below gives all negative angles
        # el = arcsin((Rprime+h) * sin(s/Rprime) / r) - (pi/2.0)

        # If elevation angle is negative, the triangle will be acute
        acute = atleast_1d( (Rprime+h)*(Rprime+h) < (Rprime*Rprime + r*r) )
        el[acute] *= -1
        el *= 180.0 / pi

        return r, el 

def get_elevation(cls, origin, destination):
        # avoid error when origin and destination are exactly on top of each other
        if destination.x == origin.x and destination.y == origin.y and destination.z == origin.z:
            return 0

        # calculate distance to destination
        dist_xy = PositionVector.get_distance_xy(origin, destination)

        # calculate elevation
        elevation = -math.atan2(destination.z-origin.z, dist_xy)
        return elevation

    # get_lon_scale - return lon scaling factor to account for curvature of earth 

def get_free_surface_elevation(self, result, free_surface, keep_details=False):
        """Compute the elevation of the free surface on a mesh for a previously solved problem.

        Parameters
        ----------
        result : LinearPotentialFlowResult
            the return of the solver
        free_surface : FreeSurface
            a meshed free surface
        keep_details : bool, optional
            if True, keep the free surface elevation in the LinearPotentialFlowResult (default:False)

        Returns
        -------
        array of shape (free_surface.nb_faces,)
            the free surface elevation on each faces of the meshed free surface

        Raises
        ------
        Exception: if the :code:`Result` object given as input does not contain the source distribution.
        """
        fs_elevation = 1j*result.omega/result.g * self.get_potential_on_mesh(result, free_surface.mesh)
        if keep_details:
            result.fs_elevation[free_surface] = fs_elevation
        return fs_elevation


# LEGACY INTERFACE 

def get_elevation(self, e, n):
        if not self.interp:
            return 0.0
        tmp = self.interp(e, n)
        if np.isnan(tmp):
            return 0.0
        else:
            return float(tmp) 

def get_elevation_data(self, num_lines=80, elevation_pts=300, viewpoint="south"):
        """Fetch elevation data and return a numpy array.

        Parameters
        ----------
        num_lines : int
            Number of horizontal lines to draw
        elevation_pts : int
            Number of points on each line to request. There's some limit to
            this that srtm enforces, but feel free to go for it!
        viewpoint : str in ["south", "west", "north", "east"] (default "south")
            The compass direction from which the map will be visualised.

        Returns
        -------
        np.ndarray
        """
        if viewpoint in ["east", "west"]:
            num_lines, elevation_pts = elevation_pts, num_lines
        values = self._srtm_data.get_image(
            (elevation_pts, num_lines), self.lats, self.longs, 5280, mode="array"
        )

        switch = {"south": 0, "west": 3, "north": 2, "east": 1}
        rotations = switch[viewpoint]
        values = np.rot90(m=values, k=rotations)
        return values 

def get_tree_height_elevation(self, loc='top'):
        ''' Sets tree height and elevation in tree dataframe

        Parameters
        ----------
        loc :    str, optional
                 initial or corrected tree top location: `top` or `top_cor`
        '''
        lons, lats = self._tree_lonlat(loc)
        self.trees[f'{loc}_height'] = self._get_z(
            lons, lats, self.chm, self.resolution)
        self.trees[f'{loc}_elevation'] = self._get_z(
            lons, lats, self.dtm, self.resolution) 

def get_elevation(gmaps, latitude, longitude):
    elevation = gmaps.elevation((latitude, longitude))
    elevation_metres = None
    if(len(elevation)>0):
        elevation_metres = elevation[0]["elevation"]
    return elevation_metres

# Get the timezone including any DST offset for the time the GPS position was recorded. 

def get_elevation_degrees(cls, scalar_label):
        """
        Returns the elevation, in degrees, corresponding to an integer
        elevation label.
        """
        scalar_label = int(scalar_label)
        assert scalar_label >= 0
        assert scalar_label < 9
        return 30 + 5 * scalar_label 

def get_elevation_angles(
        self, scans: Optional[int] = None
    ) -> Union[np.ndarray, float]:
        if scans is None:
            return self.el
        else:
            return self.el[scans] 

def get_elevation_from_latlng(self, latitude, longitude):
        """
        Get the elevation in meters of a geo point given lat/long

        :param latitude: Latitude
        :type latitude: float

        :param longitude: Longitude
        :type longitude: float
        """

        response = requests.get('https://maps.googleapis.com/maps/api/elevation/json',
                                params = {
                                    'locations': '{},{}'.format(latitude, longitude),
                                    'key': self.api_key,
                                }).json()

        elevation = None

        if response.get('results'):
            elevation = response['results'][0]['elevation']

        return { 'elevation': elevation }


# vim:sw=4:ts=4:et: 

def getBlockAtElevation(self, elevation):
        """
        Returns the block at a specified axial dimension elevation (given in cm)

        If height matches the exact top of the block, the block is considered at that
        height.

        Used as a way to determine what block the control rod will be modifying with a
        mergeBlocks.

        Parameters
        ----------
        elevation : float
            The elevation of interest to grab a block (cm)

        Returns
        -------
        targetBlock : block
            The block that exists at the specified height in the reactor
        """
        bottomOfBlock = 0.0
        for b in self:
            topOfBlock = bottomOfBlock + b.getHeight()
            if (
                topOfBlock > elevation
                or abs(topOfBlock - elevation) / elevation < 1e-10
            ) and bottomOfBlock < elevation:
                return b
            bottomOfBlock = topOfBlock
        return None 

def getElevationBoundariesByBlockType(self, blockType=None):
        """
        Gets of list of elevations, ordered from bottom to top of all boundaries of the block of specified type

        Useful for determining location of the top of the upper grid plate or active
        fuel, etc by using [0] to get the lowest boundary and [-1] to get highest

        Notes
        -----
        The list will have duplicates when blocks of the same type share a boundary.
        this is intentional. It makes it easy to grab pairs off the list and know that
        the first item in a pair is the bottom boundary and the second is the top.

        Parameters
        ----------
        blockType : str
            Block type to find. empty accepts all

        Returns
        -------
        elevation : list of floats
            Every float in the list is an elevation of a block boundary for the block
            type specified (has duplicates)

        """
        elevation, elevationsWithBlockBoundaries = 0.0, []

        # loop from bottom to top, stopping at the first instance of blockType
        for b in self:
            if b.hasFlags(blockType):
                elevationsWithBlockBoundaries.append(elevation)  # bottom Boundary
                elevationsWithBlockBoundaries.append(
                    elevation + b.getHeight()
                )  # top Boundary
            elevation += b.getHeight()

        return elevationsWithBlockBoundaries 

def get_elevation(self, interpolation='nearest_neighbor'):
        assert interpolation in ('nearest_neighbor', 'bivariate')

        if self._elevation.get(interpolation, None) is None:
            self._log.debug('Getting elevation...')
            # region = llLon, urLon, llLat, urLon
            coords = self.frame.coordinates
            lons = coords[:, 0]
            lats = coords[:, 1]

            region = (lons.min(), lons.max(), lats.min(), lats.max())
            if not srtmgl3.covers(region):
                raise AssertionError(
                    'Region is outside of SRTMGL3 topo dataset')

            tile = srtmgl3.get(region)
            if not tile:
                raise AssertionError('Cannot get SRTMGL3 topo dataset')

            if interpolation == 'nearest_neighbor':
                iy = num.rint((lats - tile.ymin) / tile.dy).astype(num.intp)
                ix = num.rint((lons - tile.xmin) / tile.dx).astype(num.intp)

                elevation = tile.data[(iy, ix)]

            elif interpolation == 'bivariate':
                interp = interpolate.RectBivariateSpline(
                    tile.y(), tile.x(), tile.data)
                elevation = interp(lats, lons, grid=False)

            elevation = elevation.reshape(self.rows, self.cols)
            self._elevation[interpolation] = elevation

        return self._elevation[interpolation] 

def get_elevation(self, lat, lon, block_size=50000):
        """
        Interpolate the elevation values for lat/lon. This is relative to the EGM96
        geoid by DTED specification.

        Parameters
        ----------
        lat : numpy.ndarray
        lon : numpy.ndarray
        block_size : None|int
            If `None`, then the entire calculation will proceed as a single block.
            Otherwise, block processing using blocks of the given size will be used.
            The minimum value used for this is 50,000, and any smaller value will be
            replaced with 50,000. Default is 50,000.

        Returns
        -------
        numpy.ndarray
            Elevation values of the same shape as lat/lon.
        """

        o_shape, lat, lon = _argument_validation(lat, lon)

        out = numpy.full(lat.shape, numpy.nan, dtype=numpy.float64)
        if block_size is None:
            boolc = self.in_bounds(lat, lon)
            if numpy.any(boolc):
                out[boolc] = self._get_elevation(lat[boolc], lon[boolc])
        else:
            block_size = min(50000, int(block_size))
            start_block = 0
            while start_block < lat.size:
                end_block = min(lat.size, start_block + block_size)
                lat1 = lat[start_block:end_block]
                lon1 = lon[start_block:end_block]
                boolc = self.in_bounds(lat1, lon1)
                out1 = numpy.full(lat1.shape, numpy.nan, dtype=numpy.float64)
                out1[boolc] = self._get_elevation(lat1[boolc], lon[boolc])
                out[start_block:end_block] = out1
                start_block = end_block

        if o_shape == ():
            return float(out[0])
        else:
            return numpy.reshape(out, o_shape) 

def get_raster_elevation(dataset, resample=None, **kwargs):
    """Return surface elevation corresponding to raster dataset
       The resampling algorithm is chosen based on scale ratio

    Parameters
    ----------
    dataset : gdal.Dataset
        raster image with georeferencing (GeoTransform at least)
    resample : GDALResampleAlg
        If None the best algorithm is chosen based on scales.
        GRA_NearestNeighbour = 0, GRA_Bilinear = 1, GRA_Cubic = 2,
        GRA_CubicSpline = 3, GRA_Lanczos = 4, GRA_Average = 5, GRA_Mode = 6,
        GRA_Max = 8, GRA_Min = 9, GRA_Med = 10, GRA_Q1 = 11, GRA_Q3 = 12
    kwargs : keyword arguments
        passed to wradlib.io.dem.get_strm()

    Returns
    -------
    elevation : :class:`numpy:numpy.ndarray`
        Array of shape (rows, cols, 2) containing elevation
    """
    extent = get_raster_extent(dataset)
    src_ds = wradlib.io.dem.get_srtm(extent, **kwargs)

    driver = gdal.GetDriverByName("MEM")
    dst_ds = driver.CreateCopy("ds", dataset)

    if resample is None:
        src_gt = src_ds.GetGeoTransform()
        dst_gt = dst_ds.GetGeoTransform()
        src_scale = min(abs(src_gt[1]), abs(src_gt[5]))
        dst_scale = min(abs(dst_gt[1]), abs(dst_gt[5]))
        ratio = dst_scale / src_scale

        resample = gdal.GRA_Bilinear
        if ratio > 2:
            resample = gdal.GRA_Average
        if ratio < 0.5:
            resample = gdal.GRA_NearestNeighbour

    gdal.ReprojectImage(
        src_ds, dst_ds, src_ds.GetProjection(), dst_ds.GetProjection(), resample
    )
    elevation = read_gdal_values(dst_ds)

    return elevation