Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
The extent of the image ``(x0, y0, x1, y1)`` referenced in
the ``img_proj`` coordinate system.
"""
url = 'https://www.nasa.gov/sites/default/files/pia17037.jpg'
img_handle = BytesIO(urlopen(url).read())
raw_image = Image.open(img_handle)
# The image is extremely high-resolution, which takes a long time to
# plot. Sub-sampling reduces the time taken to plot while not
# significantly altering the integrity of the result.
smaller_image = raw_image.resize([raw_image.size[0] // 10,
raw_image.size[1] // 10])
img = np.asarray(smaller_image)
# We define the semimajor and semiminor axes, but must also tell the
# globe not to use the WGS84 ellipse, which is its default behaviour.
img_globe = ccrs.Globe(semimajor_axis=285000., semiminor_axis=229000.,
ellipse=None)
img_proj = ccrs.PlateCarree(globe=img_globe)
img_extent = (-895353.906273091, 895353.906273091,
447676.9531365455, -447676.9531365455)
return img, img_globe, img_proj, img_extent
def _globe_from_proj4(proj4_terms):
"""Create a `Globe` object from PROJ.4 parameters."""
globe_terms = filter(lambda term: term[0] in _GLOBE_PARAMS,
proj4_terms.items())
globe = ccrs.Globe(**{_GLOBE_PARAMS[name]: value for name, value in
globe_terms})
return globe
'site_x_coordinate': site_x_coordinate,
'node_limit': node_limit,
'site_y_coordinate': site_y_coordinate
}
fargs = (site_list, orography, land_sea_mask)
kwargs = {k: v for (k, v) in args.items() if v is not None}
# Deal with coordinate systems for sites other than PlateCarree.
if 'site_coordinate_system' in kwargs.keys():
scrs = kwargs['site_coordinate_system']
if scrs not in PROJECTION_LIST:
raise ValueError('invalid projection {}'.format(scrs))
site_crs = getattr(ccrs, scrs)
scrs_opts = json.loads(kwargs.pop('site_coordinate_options', '{}'))
if 'globe' in scrs_opts:
crs_globe = ccrs.Globe(**scrs_opts['globe'])
del scrs_opts['globe']
else:
crs_globe = ccrs.Globe()
kwargs['site_coordinate_system'] = site_crs(
globe=crs_globe, **scrs_opts)
# Call plugin to generate neighbour cubes
if all_methods:
methods = [
{**kwargs, 'land_constraint': False, 'minimum_dz': False},
{**kwargs, 'land_constraint': True, 'minimum_dz': False},
{**kwargs, 'land_constraint': False, 'minimum_dz': True},
{**kwargs, 'land_constraint': True, 'minimum_dz': True}
]
all_methods = iris.cube.CubeList([])
for method in methods:
def cartopy_globe(self):
"""Initialize a `cartopy.crs.Globe` from the metadata."""
if 'earth_radius' in self._attrs:
kwargs = {'ellipse': 'sphere', 'semimajor_axis': self._attrs['earth_radius'],
'semiminor_axis': self._attrs['earth_radius']}
else:
attr_mapping = [('semimajor_axis', 'semi_major_axis'),
('semiminor_axis', 'semi_minor_axis'),
('inverse_flattening', 'inverse_flattening')]
kwargs = self._map_arg_names(self._attrs, attr_mapping)
# WGS84 with semi_major==semi_minor is NOT the same as spherical Earth
# Also need to handle the case where we're not given any spheroid
kwargs['ellipse'] = None if kwargs else 'sphere'
return ccrs.Globe(**kwargs)
if v == 'lcc':
cl = ccrs.LambertConformal
if v == 'merc':
cl = ccrs.Mercator
if v == 'utm':
cl = ccrs.UTM
if k in km_proj:
kw_proj[km_proj[k]] = v
if k in km_globe:
kw_globe[km_globe[k]] = v
if k in km_std:
kw_std[km_std[k]] = v
globe = None
if kw_globe:
globe = ccrs.Globe(**kw_globe)
if kw_std:
kw_proj['standard_parallels'] = (kw_std['lat_1'], kw_std['lat_2'])
# mercatoooor
if cl.__name__ == 'Mercator':
kw_proj.pop('false_easting', None)
kw_proj.pop('false_northing', None)
return cl(globe=globe, **kw_proj)
###########################################
# Get a Dataset view of the data (essentially a NetCDF-like interface to the
# underlying data). Pull out the data, (x, y) coordinates, and the projection
# information.
ds = f.to_dataset()
x = ds.variables['x'][:]
y = ds.variables['y'][:]
dat = ds.variables['WV']
proj_var = ds.variables[dat.grid_mapping]
print(proj_var)
###########################################
# Create CartoPy projection information for the file
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=proj_var.earth_radius,
semiminor_axis=proj_var.earth_radius)
proj = ccrs.LambertConformal(central_longitude=proj_var.longitude_of_central_meridian,
central_latitude=proj_var.latitude_of_projection_origin,
standard_parallels=[proj_var.standard_parallel],
globe=globe)
###########################################
# Plot the image
fig = plt.figure(figsize=(10, 12))
add_metpy_logo(fig, 125, 145)
ax = fig.add_subplot(1, 1, 1, projection=proj)
wv_norm, wv_cmap = ctables.registry.get_with_range('WVCIMSS', 100, 260)
wv_cmap.set_under('k')
im = ax.imshow(dat[:], cmap=wv_cmap, norm=wv_norm, zorder=0,
extent=ds.img_extent, origin='upper')
fig.clf()
glmx = glm_grids.x.data[:]
glmy = glm_grids.y.data[:]
proj_var = glm_grids['goes_imager_projection']
x = glmx * proj_var.perspective_point_height
y = glmy * proj_var.perspective_point_height
glm_xlim = x.min(), x.max()
glm_ylim = y.min(), y.max()
country_boundaries = cfeature.NaturalEarthFeature(category='cultural',
name='admin_0_countries',
scale='50m', facecolor='none')
state_boundaries = cfeature.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lakes',
scale='50m', facecolor='none')
globe = ccrs.Globe(semimajor_axis=proj_var.semi_major_axis, semiminor_axis=proj_var.semi_minor_axis)
proj = ccrs.Geostationary(central_longitude=proj_var.longitude_of_projection_origin,
satellite_height=proj_var.perspective_point_height, globe=globe)
cbars=[]
for fi, f in enumerate(fields):
glm_norm = display_params[f]['glm_norm']
product_label = display_params[f]['product_label']
file_tag = display_params[f]['file_tag']
max_format = display_params[f]['format_string']
ax = fig.add_subplot(subplots[0], subplots[1], fi+1, projection=proj)
ax.background_patch.set_facecolor(axes_facecolor)
# ax.set_aspect('auto', adjustable=None)
ax.coastlines('10m', color=map_color)
* crs (:class:`cartopy.crs.Projection`):
The coordinate reference system.
* x, y (array):
Locations at which to calculate the differentials,
defined in 'crs' coordinate reference system.
Returns:
(abs(ds/dx), abs(ds/dy)).
Numerically approximated partial differentials,
i.e. scaling factors between changes in distance and changes in
coordinate values.
"""
# Make a true-latlon coordinate system for distance calculations.
crs_latlon = ccrs.Geodetic(globe=ccrs.Globe(ellipse="sphere"))
# Transform points to true-latlon (just to get the true latitudes).
_, true_lat = _transform_xy(crs, x, y, crs_latlon)
# Get coordinate differentials w.r.t. true-latlon.
dlon_dx, dlat_dx, dlon_dy, dlat_dy = _inter_crs_differentials(
crs, x, y, crs_latlon
)
# Calculate effective scalings of X and Y coordinates.
lat_factor = np.cos(np.deg2rad(true_lat)) ** 2
ds_dx = np.sqrt(dlat_dx * dlat_dx + dlon_dx * dlon_dx * lat_factor)
ds_dy = np.sqrt(dlat_dy * dlat_dy + dlon_dy * dlon_dy * lat_factor)
return ds_dx, ds_dy