Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
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
reflon = ds.variables['lon']
reflat = ds.variables['lat']
else:
operations to do with the same grid, set ``lut`` to an existing
table obtained from a previous call to :py:meth:`Grid.grid_lookup`
return_lut : bool, optional
set to True if you want to return the lookup table for later use.
in this case, returns a tuple
Returns
-------
An aggregated ndarray of the data, in ``self`` coordinates.
If ``return_lut==True``, also return the lookup table
"""
# Input checks
if grid is None:
grid = check_crs(data) # xarray
if not isinstance(grid, Grid):
raise ValueError('grid should be a Grid instance')
if hasattr(data, 'values'):
data = data.values # xarray
# dimensional check
in_shape = data.shape
ndims = len(in_shape)
if (ndims < 2) or (ndims > 4):
raise ValueError('data dimension not accepted')
if (in_shape[-1] != grid.nx) or (in_shape[-2] != grid.ny):
raise ValueError('data dimension not compatible')
if lut is None:
lut = self.grid_lookup(grid)
# Prepare the output
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
ips -= 1
jps -= 1
we -= 1
sn -= 1
nx = we / ratio
ny = sn / ratio
if nx != (we / ratio):
raise RuntimeError('e_we and ratios are incompatible: '
else:
ny = nx * yy / xx
nx = np.rint(nx)
ny = np.rint(ny)
e, n = transform_proj(wgs84, projloc, lon, lat)
if origin== 'upper-left':
corner = (-xx / 2. + e, yy / 2. + n)
dxdy = (xx / nx, - yy / ny)
else:
corner = (-xx / 2. + e, -yy / 2. + n)
dxdy = (xx / nx, yy / ny)
return Grid(proj=projloc, x0y0=corner, nxny=(nx, ny), dxdy=dxdy,
pixel_ref='corner')
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)
we -= 1
sn -= 1
nx = we / ratio
ny = sn / ratio
if nx != (we / ratio):
raise RuntimeError('e_we and ratios are incompatible: '
'(e_we - 1) / ratio must be integer!')
if ny != (sn / ratio):
raise RuntimeError('e_sn and ratios are incompatible: '
'(e_sn - 1) / ratio must be integer!')
prevgrid = out[pid - 1]
xx, yy = prevgrid.corner_grid.x_coord, prevgrid.corner_grid.y_coord
dx = prevgrid.dx / ratio
dy = prevgrid.dy / ratio
grid = gis.Grid(nxny=(we, sn),
x0y0=(xx[ips], yy[jps]),
dxdy=(dx, dy),
pixel_ref='corner',
proj=pwrf)
out.append(grid.center_grid)
maps = None
if do_maps:
from salem import Map
import shapely.geometry as shpg
if map_kwargs is None:
map_kwargs = {}
maps = []
for i, g in enumerate(out):
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)
assert np.allclose(reflat, mylat, atol=1e-4)
if 'lon' in nc.variables:
# HAR