How to use the pyresample._spatial_mp.Proj function in pyresample

To help you get started, we’ve selected a few pyresample 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 pytroll / pyresample / pyresample / geometry.py View on Github external
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),
github pytroll / pyresample / pyresample / geo_filter.py View on Github external
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
github pytroll / pyresample / pyresample / grid.py View on Github external
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
github pytroll / pyresample / pyresample / geometry.py View on Github external
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
github pytroll / pyresample / pyresample / geometry.py View on Github external
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)
github pytroll / pyresample / pyresample / geometry.py View on Github external
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)
github pytroll / pyresample / pyresample / geometry.py View on Github external
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)
github pytroll / pyresample / pyresample / _spatial_mp.py View on Github external
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)