How to use the salem.wgs84 function in salem

To help you get started, we’ve selected a few salem examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github fmaussion / salem / salem / gis.py View on Github external
    def transform(self, x, y, z=None, crs=wgs84, nearest=False, maskout=False):
        """Converts any coordinates into the local grid.

        Parameters
        ----------
        x : ndarray
            the grid coordinates of the point(s) you want to convert
        y : ndarray
            the grid coordinates of the point(s) you want to convert
        z : None
            ignored (but necessary since some shapes have a z dimension)
        crs : crs
            reference system of x, y. Could be a pyproj.Proj instance or a
            Grid instance. In the latter case (x, y) are actually (i, j).
            (Default: lonlat in wgs84).
        nearest : bool
            set to True if you wish to return the closest i, j coordinates
github fmaussion / salem / salem / graphics.py View on Github external
-------
        the topography if needed (bonus)
        """

        if topo is None:
            self._shading_base()
            return

        kwargs.setdefault('interp', 'spline')

        if isinstance(topo, six.string_types):
            _, ext = os.path.splitext(topo)
            if ext.lower() == '.tif':
                g = GeoTiff(topo)
                # Spare memory
                ex = self.grid.extent_in_crs(crs=wgs84)  # l, r, b, t
                g.set_subset(corners=((ex[0], ex[2]), (ex[1], ex[3])),
                             crs=wgs84, margin=10)
                z = g.get_vardata()
                z[z < -999] = 0
                z = self.grid.map_gridded_data(z, g.grid, **kwargs)
            else:
                raise ValueError('File extension not recognised: {}'
                                 .format(ext))
        else:
            z = self._check_data(topo, crs=crs, **kwargs)

        # Gradient in m m-1
        ddx = self.grid.dx
        ddy = self.grid.dy
        if gis.proj_is_latlong(self.grid.proj):
            # we make a coarse approx of the avg dx on a sphere
github OGGM / oggm / oggm / core / gis.py View on Github external
if gdir.is_tidewater and cfg.PARAMS['clip_tidewater_border']:
        border = 10

    # Corners, incl. a buffer of N pix
    ulx = np.min(xx) - border * dx
    lrx = np.max(xx) + border * dx
    uly = np.max(yy) + border * dx
    lry = np.min(yy) - border * dx
    # n pixels
    nx = np.int((lrx - ulx) / dx)
    ny = np.int((uly - lry) / dx)

    # Back to lon, lat for DEM download/preparation
    tmp_grid = salem.Grid(proj=proj_out, nxny=(nx, ny), x0y0=(ulx, uly),
                          dxdy=(dx, -dx), pixel_ref='corner')
    minlon, maxlon, minlat, maxlat = tmp_grid.extent_in_crs(crs=salem.wgs84)

    # Open DEM
    source = entity.DEM_SOURCE if hasattr(entity, 'DEM_SOURCE') else None
    if not is_dem_source_available(source,
                                   (minlon, maxlon),
                                   (minlat, maxlat)):
        raise InvalidDEMError('Source: {} not available for glacier {}'
                              .format(source, gdir.rgi_id))
    dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
                                         rgi_region=gdir.rgi_region,
                                         rgi_subregion=gdir.rgi_subregion,
                                         dx_meter=dx,
                                         source=source)
    log.debug('(%s) DEM source: %s', gdir.rgi_id, dem_source)
    log.debug('(%s) N DEM Files: %s', gdir.rgi_id, len(dem_list))
github fmaussion / salem / salem / sio.py View on Github external
# OK, get it
    lon = ds.variables[x][:]
    lat = ds.variables[y][:]

    # double check for dubious variables
    if not utils.str_in_list([x], utils.valid_names['lon_var']) or \
            not utils.str_in_list([y], utils.valid_names['lat_var']):
        # name not usual. see if at least the range follows some conv
        if (np.max(np.abs(lon)) > 360.1) or (np.max(np.abs(lat)) > 90.1):
            return None

    # Make the grid
    dx = lon[1]-lon[0]
    dy = lat[1]-lat[0]
    args = dict(nxny=(lon.shape[0], lat.shape[0]), proj=wgs84, dxdy=(dx, dy),
                x0y0=(lon[0], lat[0]))
    return gis.Grid(**args)
github fmaussion / salem / salem / sio.py View on Github external
def is_rotated_proj_working():

    import pyproj
    srs = ('+ellps=WGS84 +proj=ob_tran +o_proj=latlon '
           '+to_meter=0.0174532925199433 +o_lon_p=0.0 +o_lat_p=80.5 '
           '+lon_0=357.5 +no_defs')

    p1 = pyproj.Proj(srs)
    p2 = wgs84

    return np.isclose(transform_proj(p1, p2, -20, -9),
                      [-22.243473889042903, -0.06328365194179102],
                      atol=1e-5).all()
github fmaussion / salem / salem / grids.py View on Github external
def googlestatic_mercator_grid(center_ll=None, nx=640, ny=640, zoom=12):
    """Mercator map centered on a specified point as seen by google API"""

    # Make a local proj
    lon, lat = center_ll
    proj_params = dict(proj='merc', datum='WGS84')
    projloc = pyproj.Proj(proj_params)

    # Meter per pixel
    mpix = (2 * np.pi * google_earth_radius) / google_pix / (2**zoom)
    xx = nx * mpix
    yy = ny * mpix

    e, n = pyproj.transform(wgs84, projloc, lon, lat)
    corner = (-xx / 2. + e, yy / 2. + n)
    dxdy = (xx / nx, - yy / ny)

    return Grid(proj=projloc, corner=corner, nxny=(nx, ny), dxdy=dxdy,
                pixel_ref='corner')