Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
thelons = self.lons[:, int(cols / 2)]
thelons = thelons.where(thelons.notnull(), drop=True)
thelats = self.lats[:, int(cols / 2)]
thelats = thelats.where(thelats.notnull(), drop=True)
lon1, lon2 = np.asanyarray(thelons[[0, -1]])
lines = len(thelats)
lat1, lat, lat2 = np.asanyarray(thelats[[0, int(lines / 2), -1]])
proj_dict2points = {'proj': 'omerc', 'lat_0': lat, 'ellps': ellipsoid,
'lat_1': lat1, 'lon_1': lon1,
'lat_2': lat2, 'lon_2': lon2,
'no_rot': True
}
# We need to compute alpha-based omerc for geotiff support
lonc, lat0 = Proj(**proj_dict2points)(0, 0, inverse=True)
az1, az2, _ = Geod(**proj_dict2points).inv(lonc, lat0, lon2, lat2)
azimuth = az1
az1, az2, _ = Geod(**proj_dict2points).inv(lonc, lat0, lon1, lat1)
if abs(az1 - azimuth) > 1:
if abs(az2 - azimuth) > 1:
logger.warning("Can't find appropriate azimuth.")
else:
azimuth += az2
azimuth /= 2
else:
azimuth += az1
azimuth /= 2
if abs(azimuth) > 90:
azimuth = 180 + azimuth
prj_params = {'proj': 'omerc', 'alpha': float(azimuth), 'lat_0': float(lat0), 'lonc': float(lonc),
lons (numpy array): Longitude degrees array
lats (numpy array): Latitude degrees array
Returns:
Boolean numpy array of same shape as lons and lats
"""
lons = geometry_def.lons[:]
lats = geometry_def.lats[:]
# Get projection coords
if self.nprocs > 1:
proj = _spatial_mp.Proj_MP(**self.area_def.proj_dict)
else:
proj = _spatial_mp.Proj(**self.area_def.proj_dict)
x_coord, y_coord = proj(lons, lats, nprocs=self.nprocs)
# Find array indices of coordinates
target_x = ((x_coord / self.area_def.pixel_size_x) +
self.area_def.pixel_offset_x).astype(np.int32)
target_y = (self.area_def.pixel_offset_y -
(y_coord / self.area_def.pixel_size_y)).astype(np.int32)
# Create mask for pixels outside array (invalid pixels)
target_x_valid = (target_x >= 0) & (target_x < self.area_def.width)
target_y_valid = (target_y >= 0) & (target_y < self.area_def.height)
# Set index of invalid pixels to 0
target_x[np.invert(target_x_valid)] = 0
target_y[np.invert(target_y_valid)] = 0
source_area_def : object
Source definition as AreaDefinition object
nprocs : int, optional
Number of processor cores to be used
Returns
-------
(row_indices, col_indices) : tuple of numpy arrays
Arrays for resampling area by array indexing
"""
# Proj.4 definition of source area projection
if nprocs > 1:
source_proj = _spatial_mp.Proj_MP(**source_area_def.proj_dict)
else:
source_proj = _spatial_mp.Proj(**source_area_def.proj_dict)
# get cartesian projection values from longitude and latitude
source_x, source_y = source_proj(lons, lats, nprocs=nprocs)
# Find corresponding pixels (element by element conversion of ndarrays)
source_pixel_x = (source_area_def.pixel_offset_x +
source_x / source_area_def.pixel_size_x).astype(np.int32)
source_pixel_y = (source_area_def.pixel_offset_y -
source_y / source_area_def.pixel_size_y).astype(np.int32)
return source_pixel_y, source_pixel_x
raise ValueError("Can't specify 'nprocs' and 'chunks' at the same time")
# Proj.4 definition of target area projection
proj_def = self.crs_wkt if hasattr(self, 'crs_wkt') else self.proj_dict
if hasattr(target_x, 'chunks'):
# we are using dask arrays, map blocks to th
from dask.array import map_blocks
res = map_blocks(invproj, target_x, target_y,
chunks=(target_x.chunks[0], target_x.chunks[1], 2),
new_axis=[2], proj_dict=proj_def).astype(dtype)
return res[:, :, 0], res[:, :, 1]
if nprocs > 1:
target_proj = Proj_MP(proj_def)
else:
target_proj = Proj(proj_def)
# Get corresponding longitude and latitude values
lons, lats = target_proj(target_x, target_y, inverse=True, nprocs=nprocs)
lons = np.asanyarray(lons, dtype=dtype)
lats = np.asanyarray(lats, dtype=dtype)
if cache and data_slice is None:
# Cache the result if requested
self.lons = lons
self.lats = lats
return lons, lats
def colrow2lonlat(self, cols, rows):
"""Return lons and lats for the given image columns and rows.
Both scalars and arrays are supported. To be used with scarse
data points instead of slices (see get_lonlats).
"""
p = Proj(self.proj_str)
x = self.projection_x_coords
y = self.projection_y_coords
return p(x[cols], y[rows], inverse=True)
xmax, ymax = get_geostationary_angle_extent(geos_area)
# generate points around the north hemisphere in satellite projection
# make it a bit smaller so that we stay inside the valid area
x = np.cos(np.linspace(-np.pi, 0, int(nb_points / 2.0))) * (xmax - 0.0001)
y = -np.sin(np.linspace(-np.pi, 0, int(nb_points / 2.0))) * (ymax - 0.0001)
ll_x, ll_y, ur_x, ur_y = geos_area.area_extent
x *= geos_area.proj_dict['h']
y *= geos_area.proj_dict['h']
x = np.clip(np.concatenate([x, x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x))
y = np.clip(np.concatenate([y, -y]), min(ll_y, ur_y), max(ll_y, ur_y))
return Proj(**geos_area.proj_dict)(x, y, inverse=True)
if proj_info is not None:
# this is now our complete projection information
proj_dict.update(proj_info)
projection = proj_dict
if self.optimize_projection:
return lonslats.compute_optimal_bb_area(proj_dict)
if resolution is None:
resolution = self.resolution
if shape is None:
shape = self.shape
height, width = shape
shape = None if None in shape else shape
area_extent = self.area_extent
if not area_extent or not width or not height:
proj4 = Proj(proj_dict)
try:
lons, lats = lonslats
except (TypeError, ValueError):
lons, lats = lonslats.get_lonlats()
xarr, yarr = proj4(np.asarray(lons), np.asarray(lats))
xarr[xarr > 9e29] = np.nan
yarr[yarr > 9e29] = np.nan
corners = [np.nanmin(xarr), np.nanmin(yarr),
np.nanmax(xarr), np.nanmax(yarr)]
area_extent, width, height = self.compute_domain(corners, resolution, shape)
return AreaDefinition(self.area_id, self.description, '',
projection, width, height,
area_extent, self.rotation)
def __call__(self, data1, data2, inverse=False, radians=False,
errcheck=False, nprocs=1):
if self.is_latlong():
return data1, data2
return super(Proj, self).__call__(data1, data2, inverse=inverse,
radians=radians, errcheck=errcheck)