Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# 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)
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
# 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
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: