How to use the salem.gis.transform_proj 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 / grids.py View on Github external
# Here we have to accept xarray's datasets too
    try:
        nx = len(nc.dimensions['west_east'])
        ny = len(nc.dimensions['south_north'])
    except AttributeError:
        # maybe an xarray dataset
        nx = nc.dims['west_east']
        ny = nc.dims['south_north']
    if hasattr(nc, 'PROJ_ENVI_STRING'):
        # HAR
        x0 = nc.GRID_X00
        y0 = nc.GRID_Y00
    else:
        # Normal WRF file
        e, n = gis.transform_proj(wgs84, proj, cen_lon, cen_lat)
        x0 = -(nx-1) / 2. * dx + e  # DL corner
        y0 = -(ny-1) / 2. * dy + n  # DL corner
    grid = gis.Grid(nxny=(nx, ny), ll_corner=(x0, y0), dxdy=(dx, dy),
                    proj=proj)

    if tmp_check_wrf:
        #  Temporary asserts
        if 'XLONG' in nc.variables:
            # Normal WRF
            mylon, mylat = grid.ll_coordinates
            reflon = nc.variables['XLONG']
            reflat = nc.variables['XLAT']
            if len(reflon.shape) == 3:
                reflon = reflon[0, :, :]
                reflat = reflat[0, :, :]
            assert np.allclose(reflon, mylon, atol=1e-4)
github fmaussion / salem / salem / wrftools.py View on Github external
pwrf = pwrf.format(**pargs)
    elif map_proj == 'MERCATOR':
        pwrf = '+proj=merc +lat_ts={lat_1} +lon_0={lon_0} ' \
               '+x_0=0 +y_0=0 +a=6370000 +b=6370000'
        pwrf = pwrf.format(**pargs)
    elif map_proj == 'POLAR':
        pwrf = '+proj=stere +lat_ts={lat_1} +lat_0=90.0 +lon_0={lon_0} ' \
               '+x_0=0 +y_0=0 +a=6370000 +b=6370000'
        pwrf = pwrf.format(**pargs)
    else:
        raise NotImplementedError('WRF proj not implemented yet: '
                                  '{}'.format(map_proj))
    pwrf = gis.check_crs(pwrf)

    # get easting and northings from dom center (probably unnecessary here)
    e, n = gis.transform_proj(wgs84, pwrf, pargs['ref_lon'], pargs['lat_0'])

    # LL corner
    nx, ny = e_we[0]-1, e_sn[0]-1
    x0 = -(nx-1) / 2. * dx + e  # -2 because of staggered grid
    y0 = -(ny-1) / 2. * dy + n

    # parent grid
    grid = gis.Grid(nxny=(nx, ny), x0y0=(x0, y0), dxdy=(dx, dy), proj=pwrf)

    # child grids
    out = [grid]
    for ips, jps, pid, ratio, we, sn in zip(i_parent_start, j_parent_start,
                                            parent_id, parent_ratio,
                                            e_we, e_sn):
        if ips == 1:
            continue
github fmaussion / salem / salem / sio.py View on Github external
# Here we have to accept xarray and netCDF4 datasets
    try:
        nx = len(ds.dimensions['west_east'])
        ny = len(ds.dimensions['south_north'])
    except AttributeError:
        # maybe an xarray dataset
        nx = ds.dims['west_east']
        ny = ds.dims['south_north']
    if hasattr(ds, 'PROJ_ENVI_STRING'):
        # HAR
        x0 = ds['west_east'][0]
        y0 = ds['south_north'][0]
    else:
        # Normal WRF file
        e, n = gis.transform_proj(wgs84, proj, cen_lon, cen_lat)
        x0 = -(nx-1) / 2. * dx + e  # DL corner
        y0 = -(ny-1) / 2. * dy + n  # DL corner
    grid = gis.Grid(nxny=(nx, ny), x0y0=(x0, y0), dxdy=(dx, dy), proj=proj)

    if tmp_check_wrf:
        #  Temporary asserts
        if 'XLONG' in ds.variables:
            # Normal WRF
            reflon = ds.variables['XLONG']
            reflat = ds.variables['XLAT']
        elif 'XLONG_M' in ds.variables:
            # geo_em
            reflon = ds.variables['XLONG_M']
            reflat = ds.variables['XLAT_M']
        elif 'lon' in ds.variables:
            # HAR
github fmaussion / salem / salem / datasets.py View on Github external
Notes
        -----
        To obtain the exact domain specified in `x` and `y` you may have to
        play with the `size_x` and `size_y` kwargs.
        """

        global API_KEY

        if 'zoom' in kwargs or 'center_ll' in kwargs:
            raise ValueError('incompatible kwargs.')

        # Transform to lonlat
        crs = gis.check_crs(crs)
        if isinstance(crs, pyproj.Proj):
            lon, lat = gis.transform_proj(crs, wgs84, x, y)
        elif isinstance(crs, Grid):
            lon, lat = crs.ij_to_crs(x, y, crs=wgs84)
        else:
            raise NotImplementedError()

        # surely not the smartest way to do but should be enough for now
        mc = (np.mean(lon), np.mean(lat))
        zoom = 20
        while zoom >= 0:
            grid = gis.googlestatic_mercator_grid(center_ll=mc, nx=size_x,
                                                  ny=size_y, zoom=zoom,
                                                  scale=scale)
            dx, dy = grid.transform(lon, lat, maskout=True)
            if np.any(dx.mask):
                zoom -= 1
            else: