How to use the pyproj.Geod function in pyproj

def test_geometry_area_perimeter__linestring():
    geod = Geod(ellps="WGS84")
    assert_almost_equal(
        geod.geometry_area_perimeter(LineString([Point(1, 2), Point(3, 4)])),
        (0.0, 627176.7944251911),
        decimal=2,
    )

def test_geometry_area_perimeter__linestring__radians():
    geod = Geod(ellps="WGS84")
    assert_almost_equal(
        geod.geometry_area_perimeter(
            LineString(
                [
                    Point(math.radians(1), math.radians(2)),
                    Point(math.radians(3), math.radians(4)),
                ]
            ),
            radians=True,
        ),
        (0.0, 627176.7944251911),
        decimal=2,
    )

subfault=deepcopy(subfault_list)
        else:
            subfault=deepcopy(subfault_list[ksub])
            
    
        #Assing useful object properties
        lon=subfault.longitude
        lat=subfault.latitude
        depth=subfault.depth
        strike=subfault.strike
        dip=subfault.dip
        L=subfault.length
        W=subfault.width
        
        #Projection object
        P=Geod(ellps='WGS84')
        
        #Lon lat of pseudo centers
        lon_P1,lat_P1,baz=P.fwd(lon,lat,strike,L/4.0)
        lon_P2,lat_P2,baz=P.fwd(lon,lat,strike,-L/4.0)
        
        # Pseudo widths
        wH=W*cos(deg2rad(dip))
        wV=W*sin(deg2rad(dip))
        
        #Get final lon lat depth of 4 new faults
        lon1,lat1,baz=P.fwd(lon_P1,lat_P1,strike-90,wH/4.0)
        lon2,lat2,baz=P.fwd(lon_P1,lat_P1,strike-90,-wH/4.0)
        z1=depth-wV/4.
        z2=depth+wV/4.
        
        lon3,lat3,baz=P.fwd(lon_P2,lat_P2,strike-90,wH/4.0)

# Convert UTM WGS84 coordinates to UTM NAD27.
wgs84 = pyproj.Proj('+proj=utm +zone=18 +datum=WGS84')
nad27 = pyproj.Proj('+proj=utm +zone=18 +datum=NAD27')
x, y = pyproj.transform(wgs84, nad27, 580744.32, 4504695.26)
print(x, y)


#######################  8.3.2 Great-circle calculations  #####################

# Set lat/lon coordinates for Los Angeles and Berlin.
la_lat, la_lon = 34.0500, -118.2500
berlin_lat, berlin_lon = 52.5167, 13.3833

# Create a WGS84 Geod.
geod = pyproj.Geod(ellps='WGS84')

# Get the bearings and distance between LA and Berlin
forward, back, dist = geod.inv(la_lon, la_lat, berlin_lon, berlin_lat)
print('forward: {}\nback: {}\ndist: {}'.format(forward, back, dist))

# Get your final coordinates if you start in Berlin and go dist distance in
# the back direction. These coordinates should match LA.
x, y, bearing = geod.fwd(berlin_lon, berlin_lat, back, dist)
print('{}, {}\n{}'.format(x, y, bearing))

# Get a list of equally spaced coordinates along the great circel from LA
# to Berlin.
coords = geod.npts(la_lon, la_lat, berlin_lon, berlin_lat, 100)

# Only print the first 3.
for i in range(3):

from shapely.geometry import asPolygon, Polygon
from info_dataless import locs_from_dataless
from descartes.patch import PolygonPatch
from matplotlib.colors import LogNorm
from scipy.spatial import ConvexHull
from scipy.cluster.vq import kmeans
from shapely.affinity import scale
from matplotlib.path import Path
from shapely import geometry

#------------------------------------------------------------------------------
# VARIABLES
#------------------------------------------------------------------------------

# Reference elipsoid to calculate distance.
wgs84 = pyproj.Geod(ellps='WGS84')

#------------------------------------------------------------------------------
# CLASSES
#------------------------------------------------------------------------------

class InShape:
    """
    Class defined in order to define a shapefile boundary AND quickly check
    if a given set of coordinates is contained within it. This class uses
    the shapely module.
    """    
    
    def __init__(self, input_shape, coords=0.):
        #initialise boundary shapefile location string input
        self.boundary = input_shape
        #initialise coords shape input        

from shapely import geometry

#-----------------------------------------------------------------------------
# GENERATE SECOND SET OF VARIABLES AND STATES
#-----------------------------------------------------------------------------

verbose = False
#Enter path to boundary shape file.
shape_path = "/home/boland/Dropbox/University/UniMelb\
/AGOS/PROGRAMS/ANT/Versions/26.04.2015/shapefiles/aus.shp"
# Enter number of stations required.
N = 130
# Enter km spacing between path density points.
km_points = 100.0
# Reference elipsoid to calculate distance.
wgs84 = pyproj.Geod(ellps='WGS84')
# Enter number of bins for 2D Histogram density calculation. 
nbins = 200
# Enter estimated average shear wave velocity. 3kms-1 is the default!
velocity = 3.0
# Define your ambient noise period range OR individual period in seconds.
global period_range
period_range = [1,40]

#-----------------------------------------------------------------------------
#SHAPEFILE FUNCTIONS
#-----------------------------------------------------------------------------
def shape_(input_shape):
    with fiona.open(input_shape) as fiona_collection:
        # In this case, we'll assume the shapefile only has one record/layer 
        shapefile_record = fiona_collection.next()
        # Use Shapely to create the polygon

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 shapely import geometry

#-----------------------------------------------------------------------------
# GENERATE SECOND SET OF VARIABLES AND STATES
#-----------------------------------------------------------------------------

verbose = False
#Enter path to boundary shape file.
shape_path = "/home/boland/Dropbox/University/UniMelb\
/AGOS/PROGRAMS/ANT/Versions/26.04.2015/shapefiles/aus.shp"
# Enter number of stations required.
N = 130
# Enter km spacing between path density points.
km_points = 100.0
# Reference elipsoid to calculate distance.
wgs84 = pyproj.Geod(ellps='WGS84')
# Enter number of bins for 2D Histogram density calculation. 
nbins = 200
# Enter estimated average shear wave velocity. 3kms-1 is the default!
velocity = 3.0
# Define your ambient noise period range OR individual period in seconds.
global period_range
period_range = [1,40]

#-----------------------------------------------------------------------------
#SHAPEFILE FUNCTIONS
#-----------------------------------------------------------------------------
def shape_(input_shape):
    with fiona.open(input_shape) as fiona_collection:
        # In this case, we'll assume the shapefile only has one record/layer 
        shapefile_record = fiona_collection.next()
        # Use Shapely to create the polygon

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}