Python Geod.inv方法代码示例

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
class GeodSharedMemoryBugTestIssue64(unittest.TestCase):
    def setUp(self):
        self.g = Geod(ellps='clrk66')
        self.ga = self.g.a
        self.mercury = Geod(a=2439700) # Mercury 2000 ellipsoid
                                       # Mercury is much smaller than earth.

    def test_not_shared_memory(self):
        self.assertEqual(self.ga, self.g.a)
        # mecury must have a different major axis from earth
        self.assertNotEqual(self.g.a, self.mercury.a)
        self.assertNotEqual(self.g.b, self.mercury.b)
        self.assertNotEqual(self.g.sphere, self.mercury.sphere)
        self.assertNotEqual(self.g.f, self.mercury.f)
        self.assertNotEqual(self.g.es, self.mercury.es)

        # initstrings were not shared in issue #64
        self.assertNotEqual(self.g.initstring, self.mercury.initstring)

    def test_distances(self):
        # note calculated distance was not an issue with #64, but it still a shared memory test
        boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
        portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)

        az12,az21,dist_g = self.g.inv(boston_lon,boston_lat,portland_lon,portland_lat)

        az12,az21,dist_mercury = self.mercury.inv(boston_lon,boston_lat,portland_lon,portland_lat)
        self.assertLess(dist_mercury, dist_g)

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def lat_lon_grid_deltas(longitude, latitude, **kwargs):
    r"""Calculate the delta between grid points that are in a latitude/longitude format.

    Calculate the signed delta distance between grid points when the grid spacing is defined by
    delta lat/lon rather than delta x/y

    Parameters
    ----------
    longitude : array_like
        array of longitudes defining the grid
    latitude : array_like
        array of latitudes defining the grid
    kwargs
        Other keyword arguments to pass to :class:`~pyproj.Geod`

    Returns
    -------
    dx, dy:
        at least two dimensional arrays of signed deltas between grid points in the x and y
        direction

    Notes
    -----
    Accepts 1D, 2D, or higher arrays for latitude and longitude
    Assumes [..., Y, X] for >=2 dimensional arrays

    """
    from pyproj import Geod

    # Inputs must be the same number of dimensions
    if latitude.ndim != longitude.ndim:
        raise ValueError('Latitude and longitude must have the same number of dimensions.')

    # If we were given 1D arrays, make a mesh grid
    if latitude.ndim < 2:
        longitude, latitude = np.meshgrid(longitude, latitude)

    geod_args = {'ellps': 'sphere'}
    if kwargs:
        geod_args = kwargs

    g = Geod(**geod_args)

    forward_az, _, dy = g.inv(longitude[..., :-1, :], latitude[..., :-1, :],
                              longitude[..., 1:, :], latitude[..., 1:, :])
    dy[(forward_az < -90.) | (forward_az > 90.)] *= -1

    forward_az, _, dx = g.inv(longitude[..., :, :-1], latitude[..., :, :-1],
                              longitude[..., :, 1:], latitude[..., :, 1:])
    dx[(forward_az < 0.) | (forward_az > 180.)] *= -1

    return dx * units.meter, dy * units.meter

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def test_geod_inverse_transform():
    gg = Geod(ellps="clrk66")
    lat1pt = 42.0 + (15.0 / 60.0)
    lon1pt = -71.0 - (7.0 / 60.0)
    lat2pt = 45.0 + (31.0 / 60.0)
    lon2pt = -123.0 - (41.0 / 60.0)
    """
    distance between boston and portland, clrk66:
    -66.531 75.654  4164192.708
    distance between boston and portland, WGS84:
    -66.530 75.654  4164074.239
    testing pickling of Geod instance
    distance between boston and portland, clrk66 (from pickle):
    -66.531 75.654  4164192.708
    distance between boston and portland, WGS84 (from pickle):
    -66.530 75.654  4164074.239
    inverse transform
    from proj.4 invgeod:
    b'-66.531\t75.654\t4164192.708\n'

    """
    print("from pyproj.Geod.inv:")
    az12, az21, dist = gg.inv(lon1pt, lat1pt, lon2pt, lat2pt)
    assert_almost_equal((az12, az21, dist), (-66.531, 75.654, 4164192.708), decimal=3)

    print("forward transform")
    print("from proj.4 geod:")
    endlon, endlat, backaz = gg.fwd(lon1pt, lat1pt, az12, dist)
    assert_almost_equal((endlon, endlat, backaz), (-123.683, 45.517, 75.654), decimal=3)
    print("intermediate points:")
    print("from geod with +lat_1,+lon_1,+lat_2,+lon_2,+n_S:")
    npts = 4
    lonlats = gg.npts(lon1pt, lat1pt, lon2pt, lat2pt, npts)
    lonprev = lon1pt
    latprev = lat1pt
    print(dist / (npts + 1))
    print("%6.3f  %7.3f" % (lat1pt, lon1pt))
    result_dists = (
        (-66.53059478766238, 106.79071710136431, 832838.5416198927),
        (-73.20928289863558, 99.32289055927389, 832838.5416198935),
        (-80.67710944072617, 91.36325611787134, 832838.5416198947),
        (-88.63674388212858, 83.32809401477382, 832838.5416198922),
    )
    for (lon, lat), (res12, res21, resdist) in zip(lonlats, result_dists):
        az12, az21, dist = gg.inv(lonprev, latprev, lon, lat)
        assert_almost_equal((az12, az21, dist), (res12, res21, resdist))
        latprev = lat
        lonprev = lon
    az12, az21, dist = gg.inv(lonprev, latprev, lon2pt, lat2pt)
    assert_almost_equal(
        (lat2pt, lon2pt, dist), (45.517, -123.683, 832838.542), decimal=3
    )

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def build_bins(pts, spacing=6.25, width=50, runin=0,
               isequence0=1000, ellps='WGS84'):
    """
    Build bins along a line of points.
    """
    gd = Geod(ellps=ellps)
    div_pts = divide_line(pts, spacing=spacing,
                          runin=runin, ellps=ellps)
    bins = []
    ndiv = len(div_pts)
    ibin = isequence0
    align_to_last_bin = False
    for i in range(0, ndiv - 1):
        lon0, lat0, x0, i0 = div_pts[i]
        lon1, lat1, x1, i1 = div_pts[i + 1]
        faz, baz, dist = gd.inv(lon0, lat0, lon1, lat1)
        # bin corners
        _bin = _build_bin(lon0, lat0, lon1, lat1, width, gd)
        # bin center
        _center = _calculate_center(_bin)
        # put it all together
        bins.append([ibin, None, _center, _bin])
        # handle bends in the line
        if align_to_last_bin:
            # align bins
            bins[-1][3][0] = bins[-2][3][1]
            bins[-1][3][3] = bins[-2][3][2]
            # recalculate center and offset
            bins[-1][2] = _calculate_center(bins[-1][3])
            align_to_last_bin = False
        if i0 == i1:
            ibin -= 1
            i += 1
            _bin = bins[-1]
            del bins[-1]
            align_to_last_bin = True
        # distance on line and line azimuth
        if i == 0:
            bins[-1][1] = div_pts[0][2]
            bins[-1] += [None]
        else:
            az, _, dx = gd.inv(bins[-1][2][0], bins[-1][2][1],
                              bins[-2][2][0], bins[-2][2][1])
            bins[-1][1] = bins[-2][1] + dx
            bins[-1] += [az]
        # increment bin number
        ibin += 1
    # assume first azimuth is the same as 2nd azimuth
    bins[0][4] = [bins[1][4]]
    return bins

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def calculate_distances(d):
    g = Geod(ellps='WGS84')
    far_left_lat = d['elevation_latlon_list'][0][0]
    far_left_lon = d['elevation_latlon_list'][1][0]
    for i, v in d['twod_vertices'].iteritems():
        pt_lat = v[0]
        pt_lon = v[1]
        az1, az2, dist = g.inv(far_left_lon, far_left_lat, pt_lon, pt_lat)
        km_dist = dist / 1000.0
        d['twod_vertices'][i] = [v[0], v[1], v[2], km_dist]

    d['gps_dist'] = []
    for pt_lat, pt_lon in zip(d['gps_lat'], d['gps_lon']):
        az1, az2, dist = g.inv(far_left_lon, far_left_lat, pt_lon, pt_lat)
        km_dist = dist / 1000.0
        d['gps_dist'].append(km_dist)

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
	def __init__(self, srs, bbox, width=None, height=None, format=None, resource_id=None):
		super(WmsQuery, self).__init__()
		self.query_type = 'WMS'
		self.srs = srs
		self.bbox = bbox
		self.width = width
		self.height = height
		self.format = format
		self.resource_id = resource_id
		if width is not None and height is not None:
			# calculate resolution... this should slow things down, yay... :-(
			p = Proj(init=srs.lower())
			if not p.is_latlong():
				min_lon, min_lat = p(bbox.min_x,bbox.min_y, inverse=True)
				max_lon, max_lat = p(bbox.max_x,bbox.max_y, inverse=True)
			else:
				min_lon, min_lat = bbox.min_x, bbox.min_y
				max_lon, max_lat = bbox.max_x, bbox.max_y
			g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid. 
			_,_,diagonal = g.inv(min_lon, min_lat, max_lon, max_lat)
			# distance calculated geodesic
			dist_x = sqrt(diagonal**2 / (1 + float(height)/float(width)) )
			dist_y = dist_x * (float(height)/float(width))
			self.x_res = dist_x / float(width)
			self.y_res = dist_y / float(height)
		else:
			self.x_res = None
			self.y_res = None

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def c1ompare_points_old(a, b, dx, proj='LongLat'):
    if isinstance(dx, float):
        dx = [dx, dx]
    tolerance = [0.6 * x for x in dx]
    if (proj == 'horizontal'):
        pa = project(a, projection_type='proj_cartesian')
        pb = project(b, projection_type='proj_cartesian')
        #print tolerance, pa, pb
        if ( not (abs(pa[1] - pb[1]) < tolerance[1]) ):
            return False
        elif (abs(pa[0] - pb[0]) < tolerance[0]):
            return True
        else:
            return False
    else:
        from pyproj import Geod
        wgs84_geod = Geod(ellps='WGS84')
        az12,az21,dist = wgs84_geod.inv(a[0],a[1],a[0],a[1])
        return dist < tolerance[0] * 1e5

        if ( not (abs(a[1] - b[1]) < tolerance[1]) ):
            #AddComment('lat differ')
            return False
        elif (abs(a[0] - b[0]) < tolerance[0]):
            #AddComment('long same')
            return True
        elif ((abs(abs(a[0]) - 180) < tolerance[0]) and (abs(abs(b[0]) - 180) < tolerance[0])):
            #AddComment('long +/-180')
            return True
        else:
            #AddComment('not same %g %g' % (abs(abs(a[0]) - 180), abs(abs(b[0]) - 180) ) )
            return False

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
    def gdlComp(self, lons_lats, km_pts=20):
        """
        Compute geodesic line

            lons_lats: input coordinates.
            (start longitude, start latitude, end longitude, end latitude)

            km_pts: compute one point each 20 km (default).

        """

        try:
            lon_1, lat_1, lon_2, lat_2 = lons_lats

            pygd = Geod(ellps='WGS84')

            res = pygd.inv(lon_1, lat_1, lon_2, lat_2)
            dist = res[2]

            pts  = int(math.ceil(dist) / (km_pts * 1000))

            coords = pygd.npts(lon_1, lat_1, lon_2, lat_2, pts)

            coords_se = [(lon_1, lat_1)] + coords
            coords_se.append((lon_2, lat_2))

            self.__logger.info("Geodesic line succesfully created!")
            self.__logger.info("Total points = {:,}".format(pts))
            self.__logger.info("{:,.4f} km".format(dist / 1000.))

            return coords_se

        except Exception as e:
            self.__logger.error("Error: {0}".format(e.message))

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
    def getAzimuth(self, point):
        """
		Get azimuth (in degrees) between current point and provided point (point).
		"""
        g = Geod(ellps="sphere")
        forward_azimuth, back_azimuth, distance = g.inv(self.longitude, self.latitude, point.longitude, point.latitude)
        return forward_azimuth

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def find_closest_soundings(l1b_file, tgt_latitude, tgt_longitude, max_distance, log_output=sys.stdout):
    l1b_obj = acos_file.L1B(l1b_file)

    sounding_ids = l1b_obj.get_sounding_ids()

    # Average over band first since this is likely to be consistent for
    # different types of L1B files
    latitudes  = l1b_obj.get_sounding_info('sounding_latitude')
    longitudes = l1b_obj.get_sounding_info('sounding_longitude')

    # Average over any non sounding id sized dimensions
    while len(latitudes.shape) > 1:
        extra_dims = numpy.where(numpy.array(latitudes.shape) != sounding_ids.shape[0])
        latitudes  = numpy.average(latitudes, extra_dims[0][0])
        longitudes = numpy.average(longitudes, extra_dims[0][0])

    g = Geod(ellps='WGS84')

    # Find all distances in file
    distances = numpy.zeros(len(sounding_ids), dtype=float)
    for dist_idx, lat_lon_tuple in enumerate(zip(latitudes, longitudes)):
        curr_lat, curr_lon = lat_lon_tuple 
        az12, az21, dist = g.inv(tgt_longitude,tgt_latitude,curr_lon,curr_lat)
        
        # Convert to km
        distances[dist_idx] = dist/1000.

    closest = numpy.where(distances <= max_distance)
    if len(closest[0]) > 0:
        print >>log_output, "%s" % l1b_file
        for close_idx in closest[0]:
            print >>log_output, '%d %f' % (sounding_ids[close_idx], distances[close_idx])
        print >>log_output, ""
    else:
        print >>sys.stderr, "No soundings found in %s closer than %f km" % (l1b_file, max_distance)

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def compute_lagdistances(sta,stnum,lon,lat):
    '''
    Compute the lag distances between all stations in the given set
    Input:
        sta:            Array with strings of station names (n x 1)
        stnum:          Array with unique station numbers (n x 1)
        lon:            Array with station longitudes (n x 1)
        lat:            Array with station latitudes (n x 1)
    Output:
        lagdistance:   Upper triangular matrix with lag distances for all station pairs (n x n)
    '''
    
    from pyproj import Geod
    import numpy as np
    
    # Make projection:
    p = Geod(ellps='WGS84')
    
    # Make lag matrix:
    lagdistance_full = np.zeros((len(stnum),len(stnum)))
    
    ## Start to fill in lag matrix
    # Loop over all stations, make a matrix with the lon and lat of just that station, and compute the distance to all other stations (lon,lat):
    for stationi in range(len(stnum)):
        azimuthi,backazimuthi,distancei = p.inv(lon[stationi]*np.ones(len(stnum)),lat[stationi]*np.ones(len(stnum)), lon, lat)
        
        # Fill the matrix with these distances for this station:
        lagdistance_full[stationi,:] = distancei/1000
    
    # Turn it into an upper triangular matrix:
    lagdistance = np.triu(lagdistance_full)    
                
    # Return lag distance:
    return lagdistance

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def distance_matrix(pts, lon='lon', lat='lat', ellps='WGS84'):
    """
    Calculate distance between all points

    Parameters
    ----------
    pts: pandas.DataFrame
        Table of points with at least the collumns given by ``lon`` and ``lat``
    lon, lat: str, optional
        Column names for the longitude and latitude fields. Defaults are
        'lon' and 'lat'.
    ellps: str, optional
        Name of the ellipsoid. See :class:`pyproj.Geod` for valid names.

    Returns
    -------
    distances: numpy.ndarray
        len(pts) x len(pts) array of distances between all points in ``pts``.
    """
    gd = Geod(ellps=ellps)
    npts = len(pts)
    G = np.zeros((npts, npts))
    for i in range(npts):
        for j in range(npts):
            _, _, G[i][j] = gd.inv(pts.ix[i][lon], pts.ix[i][lat],
                                   pts.ix[j][lon], pts.ix[j][lat])
    return G

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
def grcrcl1(lon_1, lat_1, lon_2, lat_2):
    
    g = Geod(ellps='WGS84')
    
    az, az_inv, dist = g.inv(lon_1, lat_1, lon_2, lat_2)
    
    return dist

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
 def midpoint_longest(north_lat, west_lon, south_lat, east_lon):
     g = Geod(ellps='WGS84')
     af, ab, dist = g.inv(west_lon, north_lat, east_lon, south_lat)
     rlon, rlat, az = g.fwd(west_lon, north_lat, af, dist/2)
     rlon += 180 if rlon < 0 else -180
     rlon = round(rlon, 6)
     rlat = round(rlat, 6)
     return rlat, rlon

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import inv [as 别名]
    def addLengthMeters(self, stream_network):
        """
        Adds length field in meters to network (The added field name will be 'LENGTH_M'). 
        
        .. note:: This may be needed for generating the kfac file depending on the units of your raster. See: :doc:`gis_tools`.       

        Parameters:
            stream_network(str): Path to stream network file.
        
        Here is an example of how to use this:
        
        .. code:: python
        
            import os
            from RAPIDpy.gis.taudem import TauDEM
            
            td = TauDEM()
                        
            output_directory = '/path/to/output/files'
            td.addLengthMeters(os.path.join(output_directory,"stream_reach_file.shp"))

        """
        network_shapefile = ogr.Open(stream_network, 1)
        network_layer = network_shapefile.GetLayer()
        network_layer_defn = network_layer.GetLayerDefn()
        
        #make sure projection EPSG:4326
        network_layer_proj = network_layer.GetSpatialRef()
        geographic_proj = osr.SpatialReference()
        geographic_proj.ImportFromEPSG(4326)
        proj_transform = None
        if network_layer_proj != geographic_proj:
            proj_transform = osr.CoordinateTransformation(network_layer_proj, geographic_proj)

        #check for field
        create_field=True
        for i in xrange(network_layer_defn.GetFieldCount()):
            field_name = network_layer_defn.GetFieldDefn(i).GetName()
            if field_name == 'LENGTH_M':
                create_field=False
                break
            
        if create_field:
            network_layer.CreateField(ogr.FieldDefn('LENGTH_M', ogr.OFTReal))

        geo_manager = Geod(ellps="WGS84")
        for network_feature in network_layer:
            feat_geom = network_feature.GetGeometryRef()
            #make sure coordinates are geographic
            if proj_transform:
                feat_geom.Transform(proj_transform)
                
            line = shapely_loads(feat_geom.ExportToWkb())
            lon_list, lat_list = line.xy
            az1, az2, dist = geo_manager.inv(lon_list[:-1], lat_list[:-1], lon_list[1:], lat_list[1:])
            network_feature.SetField('LENGTH_M', sum(dist))
            network_layer.SetFeature(network_feature)