Python Geod.npts方法代码示例

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import npts [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 npts [as 别名]
 def build_great_circle(self, start_latlng, end_latlng, distance):
     num_points = int(distance / self.ARC_LENGTH_SPACING)
     geod = Geod(ellps='WGS84')
     start_lat, start_lon = start_latlng
     end_lat, end_lon = end_latlng
     lonlats = geod.npts(start_lon, start_lat, end_lon, end_lat, num_points)
     latlngs = util.swap_pairs(lonlats)
     return latlngs

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import npts [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 npts [as 别名]
    def midpoint_utm(self,p1,p2):
        """
        calculates midpoint between 2 points on earths surface

        input: p1,p2   location in the form (lon,lat)
        return: m      midpoint in the form (lon,lat)
        """

        g = Geod(ellps='WGS84')
        l = g.npts(p1[0],p1[0],p1[1],p1[1],1)
        
        m = l[0]
        # print 'm: ', m
        
        return m

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import npts [as 别名]
def geodesic(crs, start, end, steps):
    r"""Construct a geodesic path between two points.

    This function acts as a wrapper for the geodesic construction available in `pyproj`.

    Parameters
    ----------
    crs: `cartopy.crs`
        Cartopy Coordinate Reference System to use for the output
    start: (2, ) array_like
        A latitude-longitude pair designating the start point of the geodesic (units are
        degrees north and degrees east).
    end: (2, ) array_like
        A latitude-longitude pair designating the end point of the geodesic (units are degrees
        north and degrees east).
    steps: int, optional
        The number of points along the geodesic between the start and the end point
        (including the end points).

    Returns
    -------
    `numpy.ndarray`
        The list of x, y points in the given CRS of length `steps` along the geodesic.

    See Also
    --------
    cross_section

    """
    import cartopy.crs as ccrs
    from pyproj import Geod

    # Geod.npts only gives points *in between* the start and end, and we want to include
    # the endpoints.
    g = Geod(crs.proj4_init)
    geodesic = np.concatenate([
        np.array(start[::-1])[None],
        np.array(g.npts(start[1], start[0], end[1], end[0], steps - 2)),
        np.array(end[::-1])[None]
    ]).transpose()
    points = crs.transform_points(ccrs.Geodetic(), *geodesic)[:, :2]

    return points

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import npts [as 别名]
def ctr_generator(lons, lats, outfname, d = 10.):
    if lons.size != lats.size:
        raise ValueError('Size of longitude and latitude list must be the same')
    g           = Geod(ellps='WGS84')
    N           = lons.size
    lonlats     = []
    for i in range(N):
        lon1    = lons[i]
        lat1    = lats[i]
        if i < N-1:
            lon2    = lons[i+1]
            lat2    = lats[i+1]
        else:
            lon2    = lons[0]
            lat2    = lats[0]
        az, baz, dist   = g.inv(lon1, lat1, lon2, lat2)
        dist            = dist/1000.
        if d < dist:
            d               = dist/float(int(dist/d))
            Nd              = int(dist/d)
            lonlats         += [(lon1, lat1)]
            lonlats         += g.npts(lon1, lat1, lon2, lat2, npts=Nd-1)
        else:
            lonlats         += [(lon1, lat1)]
    with open(outfname, 'w') as fid:
        npts        = len(lonlats)
        fid.writelines('0. 0. \n')
        fid.writelines('%g \n' %npts)
        for lonlat in lonlats:
            if lonlat[0] < 0.:
                outlon  = lonlat[0]+360.
            else:
                outlon  = lonlat[0]
            outlat      = lonlat[1]
            fid.writelines('%g  %g\n' %(outlon, outlat))
        fid.writelines('%g \n' %npts)
        for i in range(npts):
            if i < npts-1:
                fid.writelines('%g  %g\n' %(i+1, i+2))
            else:
                fid.writelines('%g  %g\n' %(i+1, 1))

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import npts [as 别名]
class TestRadians(unittest.TestCase):
    """Tests issue #84"""
    def setUp(self):
        self.g = Geod(ellps='clrk66')
        self.boston_d = (-71. - (7. / 60.), 42. + (15. / 60.))
        self.boston_r = (math.radians(self.boston_d[0]), math.radians(self.boston_d[1]))
        self.portland_d = (-123. - (41. / 60.), 45. + (31. / 60.))
        self.portland_r = (math.radians(self.portland_d[0]), math.radians(self.portland_d[1]))

    def test_inv_radians(self):

        # Get bearings and distance from Boston to Portland in degrees
        az12_d, az21_d, dist_d = self.g.inv(
            self.boston_d[0],
            self.boston_d[1],
            self.portland_d[0],
            self.portland_d[1],
            radians=False)

        # Get bearings and distance from Boston to Portland in radians
        az12_r, az21_r, dist_r = self.g.inv(
            self.boston_r[0],
            self.boston_r[1],
            self.portland_r[0],
            self.portland_r[1],
            radians=True)

        # Check they are equal
        self.assertAlmostEqual(az12_d, math.degrees(az12_r))
        self.assertAlmostEqual(az21_d, math.degrees(az21_r))
        self.assertAlmostEqual(dist_d, dist_r)

    def test_fwd_radians(self):
        # Get bearing and distance to Portland
        az12_d, az21_d, dist = self.g.inv(
            self.boston_d[0],
            self.boston_d[1],
            self.portland_d[0],
            self.portland_d[1],
            radians=False)

        # Calculate Portland's lon/lat from bearing and distance in degrees
        endlon_d, endlat_d, backaz_d = self.g.fwd(
            self.boston_d[0],
            self.boston_d[1],
            az12_d,
            dist,
            radians=False)

        # Calculate Portland's lon/lat from bearing and distance in radians
        endlon_r, endlat_r, backaz_r = self.g.fwd(
            self.boston_r[0],
            self.boston_r[1],
            math.radians(az12_d),
            dist,
            radians=True)

        # Check they are equal
        self.assertAlmostEqual(endlon_d, math.degrees(endlon_r))
        self.assertAlmostEqual(endlat_d, math.degrees(endlat_r))
        self.assertAlmostEqual(backaz_d, math.degrees(backaz_r))

        # Check to make sure we're back in Portland
        self.assertAlmostEqual(endlon_d, self.portland_d[0])
        self.assertAlmostEqual(endlat_d, self.portland_d[1])

    def test_npts_radians(self):
        # Calculate 10 points between Boston and Portland in degrees
        points_d = self.g.npts(
            lon1=self.boston_d[0],
            lat1=self.boston_d[1],
            lon2=self.portland_d[0],
            lat2=self.portland_d[1],
            npts=10,
            radians=False)

        # Calculate 10 points between Boston and Portland in radians
        points_r = self.g.npts(
            lon1=self.boston_r[0],
            lat1=self.boston_r[1],
            lon2=self.portland_r[0],
            lat2=self.portland_r[1],
            npts=10,
            radians=True)

        # Check they are equal
        for index, dpoint in enumerate(points_d):
            self.assertAlmostEqual(dpoint[0], math.degrees(points_r[index][0]))
            self.assertAlmostEqual(dpoint[1], math.degrees(points_r[index][1]))

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import npts [as 别名]

#.........这里部分代码省略.........
            count       the frequency of each pair

        Once the data is set, call draw() to render the image.
        '''
        self.data = np.array((lon1, lat1, lon2, lat2)).T
        self.data_size = self.data.shape[0]

        # calculate great-circle distances between each coordinate
        # pair. This will be used for determining which order to
        # draw the lines in.
        _, _, dist = self.geo.inv(np.array(lon1),
                                  np.array(lat1),
                                  np.array(lon2),
                                  np.array(lat2))

        # if a frequency count is given, take it into account when
        # creating weights. Weights determine the coloring and
        # drawing order of each line.
        if count is not None:
            self.weight = np.array(count) / dist
        else:
            self.weight = 1 / dist

        self.order = np.argsort(self.weight)

    def draw(self):
        '''
        Render the image. Assumes that set_data has already been called.

        Returns a Python Image Library (PIL) Image object.
        '''

        img = Image.new('RGB', (self.width, self.height), self.bgcol)
        canvas = Draw(img)

        # create the projection. Here we use an equidistant cylindrical projection,
        # but others may work with tweaking of the parameters.
        proj = Proj(proj=self.proj,
                a=self.width/(2*pi), # set the radius of the earth such that our
                                     # projections work
                x_0=self.width/2,    # center horizontally on the image
                y_0=self.height/2)   # center verticallly on the image

        # two branches below will use the same sequence of commands to
        # draw a great-circle on the map, so the common elements are wrapped
        # up into a locally defined function. Given a matrix of points and
        # a pen, draw the path through the points.
        def draw_(pts, pen):
            lons, lats = pts.T
            x, y = proj(lons, lats)
            y = self.height - y
            path = reduce(operator.add, zip(x, y))
            canvas.line(path, pen)

        # loop over every coordinate pair
        for i, (lon1, lat1, lon2, lat2) in enumerate(self.data[self.order]):
            # calculate the fraction of the paths already drawn, and use
            # it to create a pen of the appropriate color
            frac = i / float(self.data_size)
            pen = Pen(self.cols(frac), self.line_width)

            # find the intermediate coordinates along a line between the two 
            # coordinates
            pts = self.geo.npts(lon1, lat1, lon2, lat2, self.gc_resolution)
            pts = np.array(pts)

            # if the longitudinal distance between the two points (travelling
            # through the prime meridian) is more than 180 degrees, it's faster
            # to *not* travel through the prime meridian, so we have to special-
            # case the drawing of the lines.
            if abs(lon1 - lon2) >= HALF_ROTATION:
                # find the index of the path where the line wraps around the image
                (cut_point,), = np.where(np.abs(np.diff(pts[:,0])) > HALF_ROTATION)
                
                # draw the two resultant lines separately
                pts1 = pts[:cut_point+1,:]
                pts2 = pts[cut_point+1:,:]

                # plot one point after the break on each sides so that the
                # paths go to the edge of the screen
                x1, y1 = pts[cut_point+2, :]
                x2, y2 = pts[cut_point+1, :]

                if x1 > 0:
                    pts1 = np.vstack((pts1, [-HALF_ROTATION, y1]))
                    pts2 = np.vstack(([HALF_ROTATION, y2], pts2))
                else:
                    pts1 = np.vstack((pts1, [HALF_ROTATION, y1]))
                    pts2 = np.vstack(([-HALF_ROTATION, y2], pts2))

                draw_(pts1, pen)
                draw_(pts2, pen)
            else:
                # the path does not wrap the image, so we can simply draw
                # it as-is
                draw_(pts, pen)
            
        canvas.flush()

        return img

# 需要导入模块: from pyproj import Geod [as 别名]
# 或者: from pyproj.Geod import npts [as 别名]
az12,az21,dist = g.inv(lon1pt,lat1pt,lon2pt,lat2pt)
print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
print 'forward transform'
print 'from proj.4 geod:'
print commands.getoutput('echo "42d15\'N 71d07\'W -66d31\'50.141 4164192.708" | geod +ellps=clrk66 -f "%.3f"')
endlon,endlat,backaz = g.fwd(lon1pt,lat1pt,az12,dist)
print 'from pyproj.Geod.fwd:'
print "%6.3f  %6.3f %13.3f" % (endlat,endlon,backaz)
print 'intermediate points:'
print 'from geod with +lat_1,+lon_1,+lat_2,+lon_2,+n_S:'
points = '+lon_1=%s +lat_1=%s +lon_2=%s +lat_2=%s' % (lon1pt,lat1pt,lon2pt,lat2pt,)
print points
print commands.getoutput('geod +ellps=clrk66 -f "%.3f" +n_S=5 '+points)
print 'from pyproj.Geod.npts:'
npts = 4
lonlats = g.npts(lon1pt,lat1pt,lon2pt,lat2pt,npts)
lonprev = lon1pt
latprev = lat1pt
print dist/(npts+1)
print '%6.3f  %7.3f' % (lat1pt, lon1pt)
for lon, lat in lonlats:
    az12,az21,dist = g.inv(lonprev,latprev,lon,lat)
    print '%6.3f  %7.3f %11.3f' % (lat, lon, dist)
    latprev = lat; lonprev = lon
az12,az21,dist = g.inv(lonprev,latprev,lon2pt,lat2pt)
print '%6.3f  %7.3f %11.3f' % (lat2pt, lon2pt, dist)

# specify the lat/lons of some cities.
boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)