Python pyproj.Geod() Examples

def distance(lat1: T, lon1: T, lat2: T, lon2: T, *args, **kwargs) -> T:
    geod = Geod(ellps="WGS84")
    angle1, angle2, dist1 = geod.inv(lon1, lat1, lon2, lat2, *args, **kwargs)
    return dist1 

def distMatrix(inCoords, distanceMetric=False):
    """
    Compute distance matrix between points
    coords : nparray shape[nPoints,2], with first column X, and Y. Proj 4326(WGS84)
    Return matrix of distance matrix between points.
    """
    if distanceMetric:
        from pyproj import Geod
        geod = Geod(ellps='WGS84')

        distArray = np.zeros((len(inCoords), len(inCoords)))
        for n, p in enumerate(np.nditer(inCoords.T.copy(), flags=[
                              'external_loop'], order='F')):
            for i in range(len(inCoords)):
                x1, y1 = p
                x2, y2 = inCoords[i]
                angle1, angle2, dist = geod.inv(x1, y1, x2, y2)

                distArray[n, i] = dist

    else:
        from scipy.spatial import distance

        distArray = distance.cdist(inCoords, inCoords, 'euclidean')

    return distArray 

def bearing(lat1: T, lon1: T, lat2: T, lon2: T, *args, **kwargs) -> T:
    geod = Geod(ellps="WGS84")
    angle1, angle2, dist1 = geod.inv(lon1, lat1, lon2, lat2, *args, **kwargs)
    return angle1 

def destination(
    lat: T, lon: T, bearing: T, distance: T, *args, **kwargs
) -> Tuple[T, T, T]:
    geod = Geod(ellps="WGS84")
    lon_, lat_, back_ = geod.fwd(lon, lat, bearing, distance, *args, **kwargs)
    return lat_, lon_, back_ 

def greatcircle(
    lat1: T, lon1: T, lat2: T, lon2: T, *args, **kwargs
) -> List[Tuple[T, T]]:
    geod = Geod(ellps="WGS84")
    return [
        (lat, lon)
        for (lon, lat) in geod.npts(lon1, lat1, lon2, lat2, *args, **kwargs)
    ] 

def _compute_uniform_shape(self):
        """Compute the height and width of a domain to have uniform resolution across dimensions."""
        g = Geod(ellps='WGS84')

        def notnull(arr):
            try:
                return arr.where(arr.notnull(), drop=True)
            except AttributeError:
                return arr[np.isfinite(arr)]
        leftlons = self.lons[:, 0]
        rightlons = self.lons[:, -1]
        middlelons = self.lons[:, int(self.lons.shape[1] / 2)]
        leftlats = self.lats[:, 0]
        rightlats = self.lats[:, -1]
        middlelats = self.lats[:, int(self.lats.shape[1] / 2)]
        try:
            import dask.array as da
        except ImportError:
            pass
        else:
            leftlons, rightlons, middlelons, leftlats, rightlats, middlelats = da.compute(leftlons, rightlons,
                                                                                          middlelons, leftlats,
                                                                                          rightlats, middlelats)
        leftlons = notnull(leftlons)
        rightlons = notnull(rightlons)
        middlelons = notnull(middlelons)
        leftlats = notnull(leftlats)
        rightlats = notnull(rightlats)
        middlelats = notnull(middlelats)

        az1, az2, width1 = g.inv(leftlons[0], leftlats[0], rightlons[0], rightlats[0])
        az1, az2, width2 = g.inv(leftlons[-1], leftlats[-1], rightlons[-1], rightlats[-1])
        az1, az2, height = g.inv(middlelons[0], middlelats[0], middlelons[-1], middlelats[-1])
        width = min(width1, width2)
        vresolution = height * 1.0 / self.lons.shape[0]
        hresolution = width * 1.0 / self.lons.shape[1]
        resolution = min(vresolution, hresolution)
        width = int(width * 1.1 / resolution)
        height = int(height * 1.1 / resolution)
        return height, width 

def get_bounding_box(self):
        """Get the bounding box of this file."""
        from pyproj import Geod
        geod = Geod(ellps='WGS84')
        dataset_group = DATASET_KEYS[self.datasets[0]]
        idx = 0
        lons_ring = None
        lats_ring = None
        while True:
            path = 'Data_Products/{dataset_group}/{dataset_group}_Gran_{idx}/attr/'
            prefix = path.format(dataset_group=dataset_group, idx=idx)
            try:
                lats = self.file_content[prefix + 'G-Ring_Latitude']
                lons = self.file_content[prefix + 'G-Ring_Longitude']
                if lons_ring is None:
                    lons_ring = lons
                    lats_ring = lats
                else:
                    prev_lon = lons_ring[0]
                    prev_lat = lats_ring[0]
                    dists = list(geod.inv(lon, lat, prev_lon, prev_lat)[2] for lon, lat in zip(lons, lats))
                    first_idx = np.argmin(dists)
                    if first_idx == 2 and len(lons) == 8:
                        lons_ring = np.hstack((lons[:3], lons_ring[:-2], lons[4:]))
                        lats_ring = np.hstack((lats[:3], lats_ring[:-2], lats[4:]))
                    else:
                        raise NotImplementedError("Don't know how to handle G-Rings of length %d" % len(lons))

            except KeyError:
                break
            idx += 1

        return lons_ring, lats_ring 

def toLonLatAlt(self, r, az, el):
        """Convert slant range r, azimuth az, and elevation el to ECEF system"""
        geoSys = GeographicSystem()
        geodetic = proj4.Geod(ellps=self.ellps)

        try:
            n = max((az.size, r.size))
        except AttributeError:
            n = max((len(az), len(r)))

        dist, z = self.getGroundRangeHeight(r,el)
        lon, lat, backAz = geodetic.fwd([self.ctrLon]*n, [self.ctrLat]*n, az, dist)
        return lon, lat, z 

def fromECEF(self, x, y, z):
        """Convert ECEF system to slant range r, azimuth az, and elevation el"""
        # x = np.atleast1d(x)
        geoSys = GeographicSystem()
        geodetic = proj4.Geod(ellps=self.ellps)

        try:
            n = x.size
        except AttributeError:
            n = len(x)

        lon, lat, z = geoSys.fromECEF(x, y, z)
        radarToGateAz, gateToRadarAz, dist = geodetic.inv([self.ctrLon]*n, [self.ctrLat]*n, lon, lat)
        az = array(radarToGateAz)   #radarToGateAz may be a list.
        # change negative azimuths to positive
        az[az < 0.0] += 360.0

        #have height, ground range, azimuth. need to get elev angle and slant range from ground range and height
        r, el = self.getSlantRangeElevation(dist, z)

        return r, az, el 

def expandRect(minLat, minLong, maxLat, maxLong, distance):
    geod = pyproj.Geod(ellps="WGS84")
    midLat  = (minLat + maxLat) / 2.0
    midLong = (minLong + maxLong) / 2.0

    try:
        availDistance = geod.inv(midLong, maxLat, midLong,
                                 +90)[2]
        if availDistance >= distance:
            x,y,angle = geod.fwd(midLong, maxLat, 0, distance)
            maxLat = y
        else:
            maxLat = +90
    except:
        maxLat = +90 # Can't expand north.

    try:
        availDistance = geod.inv(maxLong, midLat, +180,
                                 midLat)[2]
        if availDistance >= distance:
            x,y,angle = geod.fwd(maxLong, midLat, 90,
                                 distance)
            maxLong = x
        else:
            maxLong = +180
    except:
        maxLong = +180 # Can't expand east.

    try:
        availDistance = geod.inv(midLong, minLat, midLong,
                                 -90)[2]
        if availDistance >= distance:
            x,y,angle = geod.fwd(midLong, minLat, 180,
                                 distance)
            minLat = y
        else:
            minLat = -90
    except:
        minLat = -90 # Can't expand south.

    try:
        availDistance = geod.inv(maxLong, midLat, -180,
                                 midLat)[2]
        if availDistance >= distance:
            x,y,angle = geod.fwd(minLong, midLat, 270,
                                 distance)
            minLong = x
        else:
            minLong = -180
    except:
        minLong = -180 # Can't expand west.

    return (minLat, minLong, maxLat, maxLong)

############################################################################# 

def spark_matchup_driver(tile_ids, bounding_wkt, primary_ds_name, matchup_ds_names, parameter, depth_min, depth_max,
                         time_tolerance, radius_tolerance, platforms, match_once, sc=None):
    from functools import partial

    with DRIVER_LOCK:
        # Broadcast parameters
        primary_b = sc.broadcast(primary_ds_name)
        matchup_b = sc.broadcast(matchup_ds_names)
        depth_min_b = sc.broadcast(float(depth_min) if depth_min is not None else None)
        depth_max_b = sc.broadcast(float(depth_max) if depth_max is not None else None)
        tt_b = sc.broadcast(time_tolerance)
        rt_b = sc.broadcast(float(radius_tolerance))
        platforms_b = sc.broadcast(platforms)
        bounding_wkt_b = sc.broadcast(bounding_wkt)
        parameter_b = sc.broadcast(parameter)

        # Parallelize list of tile ids
        rdd = sc.parallelize(tile_ids, determine_parllelism(len(tile_ids)))

    # Map Partitions ( list(tile_id) )
    rdd_filtered = rdd.mapPartitions(
        partial(match_satellite_to_insitu, primary_b=primary_b, matchup_b=matchup_b, parameter_b=parameter_b, tt_b=tt_b,
                rt_b=rt_b, platforms_b=platforms_b, bounding_wkt_b=bounding_wkt_b, depth_min_b=depth_min_b,
                depth_max_b=depth_max_b), preservesPartitioning=True) \
        .filter(lambda p_m_tuple: abs(
        iso_time_to_epoch(p_m_tuple[0].time) - iso_time_to_epoch(p_m_tuple[1].time)) <= time_tolerance)

    if match_once:
        # Only the 'nearest' point for each primary should be returned. Add an extra map/reduce which calculates
        # the distance and finds the minimum

        # Method used for calculating the distance between 2 DomsPoints
        from pyproj import Geod

        def dist(primary, matchup):
            wgs84_geod = Geod(ellps='WGS84')
            lat1, lon1 = (primary.latitude, primary.longitude)
            lat2, lon2 = (matchup.latitude, matchup.longitude)
            az12, az21, distance = wgs84_geod.inv(lon1, lat1, lon2, lat2)
            return distance

        rdd_filtered = rdd_filtered \
            .map(lambda (primary, matchup): tuple([primary, tuple([matchup, dist(primary, matchup)])])) \
            .reduceByKey(lambda match_1, match_2: match_1 if match_1[1] < match_2[1] else match_2) \
            .mapValues(lambda x: [x[0]])
    else:
        rdd_filtered = rdd_filtered \
            .combineByKey(lambda value: [value],  # Create 1 element list
                          lambda value_list, value: value_list + [value],  # Add 1 element to list
                          lambda value_list_a, value_list_b: value_list_a + value_list_b)  # Add two lists together

    result_as_map = rdd_filtered.collectAsMap()

    return result_as_map 

def proj4_radius_parameters(proj4_dict):
    """Calculate 'a' and 'b' radius parameters.

    Arguments:
        proj4_dict (str or dict): PROJ.4 parameters

    Returns:
        a (float), b (float): equatorial and polar radius
    """
    if CRS is not None:
        import math
        crs = CRS(proj4_dict)
        a = crs.ellipsoid.semi_major_metre
        b = crs.ellipsoid.semi_minor_metre
        if not math.isnan(b):
            return a, b
        # older versions of pyproj didn't always have a valid minor radius
        proj4_dict = crs.to_dict()

    if isinstance(proj4_dict, str):
        new_info = proj4_str_to_dict(proj4_dict)
    else:
        new_info = proj4_dict.copy()

    # load information from PROJ.4 about the ellipsis if possible

    from pyproj import Geod

    if 'ellps' in new_info:
        geod = Geod(**new_info)
        new_info['a'] = geod.a
        new_info['b'] = geod.b
    elif 'a' not in new_info or 'b' not in new_info:

        if 'rf' in new_info and 'f' not in new_info:
            new_info['f'] = 1. / float(new_info['rf'])

        if 'a' in new_info and 'f' in new_info:
            new_info['b'] = float(new_info['a']) * (1 - float(new_info['f']))
        elif 'b' in new_info and 'f' in new_info:
            new_info['a'] = float(new_info['b']) / (1 - float(new_info['f']))
        elif 'R' in new_info:
            new_info['a'] = new_info['R']
            new_info['b'] = new_info['R']
        else:
            geod = Geod(**{'ellps': 'WGS84'})
            new_info['a'] = geod.a
            new_info['b'] = geod.b

    return float(new_info['a']), float(new_info['b']) 

def _compute_omerc_parameters(self, ellipsoid):
        """Compute the oblique mercator projection bouding box parameters."""
        lines, cols = self.lons.shape
        lon1, lon2 = np.asanyarray(self.lons[[0, -1], int(cols / 2)])
        lat1, lat, lat2 = np.asanyarray(
            self.lats[[0, int(lines / 2), -1], int(cols / 2)])
        if any(np.isnan((lon1, lon2, lat1, lat, lat2))):
            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),
                      'gamma': 0,
                      'ellps': ellipsoid}

        return prj_params 

def _find_zcta_closest_isd_stations(zcta_metadata, isd_station_metadata, limit=None):
    if limit is None:
        limit = 10
    import pyproj

    geod = pyproj.Geod(ellps="WGS84")

    isd_usaf_ids, isd_lats, isd_lngs = zip(
        *[
            (
                isd_station["usaf_id"],
                float(isd_station["latitude"]),
                float(isd_station["longitude"]),
            )
            for isd_station in isd_station_metadata.values()
        ]
    )

    isd_lats = np.array(isd_lats)
    isd_lngs = np.array(isd_lngs)

    for zcta in zcta_metadata.values():
        zcta_lats = np.tile(zcta["latitude"], isd_lats.shape)
        zcta_lngs = np.tile(zcta["longitude"], isd_lngs.shape)

        dists = geod.inv(zcta_lngs, zcta_lats, isd_lngs, isd_lats)[2]
        sorted_dists = np.argsort(dists)[:limit]

        closest_isd_stations = []
        for i, idx in enumerate(sorted_dists):
            usaf_id = isd_usaf_ids[idx]
            isd_station = isd_station_metadata[usaf_id]
            closest_isd_stations.append(
                {
                    "usaf_id": usaf_id,
                    "distance_meters": int(round(dists[idx])),
                    "rank": i + 1,
                    "iecc_climate_zone_match": (
                        zcta.get("iecc_climate_zone")
                        == isd_station.get("iecc_climate_zone")
                    ),
                    "iecc_moisture_regime_match": (
                        zcta.get("iecc_moisture_regime")
                        == isd_station.get("iecc_moisture_regime")
                    ),
                    "ba_climate_zone_match": (
                        zcta.get("ba_climate_zone")
                        == isd_station.get("ba_climate_zone")
                    ),
                    "ca_climate_zone_match": (
                        zcta.get("ca_climate_zone")
                        == isd_station.get("ca_climate_zone")
                    ),
                }
            )
        zcta["closest_isd_stations"] = closest_isd_stations 

def distance(
    da: Union[xarray.DataArray, xarray.Dataset],
    *,
    lon: Union[float, Sequence[float], xarray.DataArray],
    lat: Union[float, Sequence[float], xarray.DataArray],
):
    """Return distance to a point in meters.

    Parameters
    ----------
    da : Union[xarray.DataArray, xarray.Dataset]
      Input data.
    lon : Union[float, Sequence[float], xarray.DataArray]
      Longitude coordinate.
    lat : Union[float, Sequence[float], xarray.DataArray]
      Latitude coordinate.

    Returns
    -------
    xarray.DataArray
      Distance in meters to point.

    Note
    ----
    To get the indices from closest point, use:
    >>> da = xr.open_dataset(path_to_pr_file).pr
    >>> d = distance(da, lon=-75, lat=45)
    >>> k = d.argmin()
    >>> i, j, _ = np.unravel_index(k, d.shape)
    """
    ptdim = lat.dims[0]

    g = Geod(ellps="WGS84")  # WGS84 ellipsoid - decent globally

    def func(lons, lats, lon, lat):
        return g.inv(lons, lats, lon, lat)[2]

    out = xarray.apply_ufunc(
        func,
        *xarray.broadcast(da.lon.load(), da.lat.load(), lon, lat),
        input_core_dims=[[ptdim]] * 4,
        output_core_dims=[[ptdim]],
    )
    out.attrs["units"] = "m"
    return out