How to use the pyproj.transform function in pyproj

'PARAMETER["longitude_of_center",21.5],'
                     'PARAMETER["latitude_of_center",8.5],UNIT["Meter",1.0]]')
        aeqd_sr = osr.SpatialReference()
        aeqd_sr.ImportFromWkt(aeqd_proj)
        p_grid = pyproj.Proj(aeqd_sr.ExportToProj4())
        p_geo = pyproj.Proj(geo_sr.ExportToProj4())
    
        # test locations in Africa
        points = [(-31.627336, 30.306273),
                  (-14.589038, -43.880131),
                  (79.423313, -35.261658),
                  (23.456413, 10.457987)]
    
        for i, pt in enumerate(points):
            # from lat/lon to aeqd
            aeqd_x, aeqd_y = pyproj.transform(p_geo, p_grid, pt[0], pt[1])
            # from aeqd to lat/lon
            lon, lat = pyproj.transform(p_grid, p_geo, aeqd_x, aeqd_y)
            # print info
    
            print("testing location {}:".format(i))
            _info = "   ({:f},{:f}) -> ({:f},{:f}) -> ({:f},{:f})"
            print(_info.format(pt[0], pt[1], aeqd_x, aeqd_y, lon, lat))
            print("    difference: ({:f},{:f})".format(lon - pt[0], lat - pt[1]))
            nptest.assert_allclose(pt[0], lon)
            nptest.assert_allclose(pt[1], lat)

def _get_chip(self, window):
        no_shift = self.x_shift_meters == 0.0 and self.y_shift_meters == 0.0
        yes_shift = not no_shift
        if yes_shift:
            ymin, xmin, ymax, xmax = window.tuple_format()
            width = window.get_width()
            height = window.get_height()

            # Transform image coordinates into world coordinates
            transform = self.image_dataset.transform
            xmin2, ymin2 = transform * (xmin, ymin)

            # Transform from world coordinates to WGS84
            if self.crs != wgs84_proj4 and self.proj:
                lon, lat = pyproj.transform(self.proj, wgs84, xmin2, ymin2)
            else:
                lon, lat = xmin2, ymin2

            # Shift.  This is performed by computing the shifts in
            # meters to shifts in degrees.  Those shifts are then
            # applied to the WGS84 coordinate.
            #
            # Courtesy of https://gis.stackexchange.com/questions/2951/algorithm-for-offsetting-a-latitude-longitude-by-some-amount-of-meters  # noqa
            lat_radians = math.pi * lat / 180.0
            dlon = Decimal(self.x_shift_meters) / Decimal(
                meters_per_degree * math.cos(lat_radians))
            dlat = Decimal(self.y_shift_meters) / Decimal(meters_per_degree)
            lon = float(Decimal(lon) + dlon)
            lat = float(Decimal(lat) + dlat)

            # Transform from WGS84 to world coordinates

# The lat long are in epsg 4326 which is WGS84
        inProj = Proj(init='epsg:4326')
        outProj = Proj(init=EPSG_string)
        this_Lat = self.Latitude[0]
        this_Lon = self.Longitude[0]
        
        print("Lat-long: ")
        print(this_Lat)
        print(this_Lon)
        
        easting =[]
        northing = []
        
        for idx, Lon in enumerate(self.Longitude):
            Lat = self.Latitude[idx]
            ea,no = transform(inProj,outProj,Lon,Lat)
            easting.append(ea)
            northing.append(no)
            
        return easting,northing

def transform_point(in_point, in_crs, out_crs):
    """
    export a Shapely geom to GeoJSON Feature and
    transform to a new coordinate system with pyproj
    :param in_point: shapely geometry as point
    :param in_crs: pyproj crs definition
    :param out_crs:  pyproj output crs definition
    :return: GeoJSON transformed to out_crs
    """
    geojs_geom = in_point.__geo_interface__

    x1 = geojs_geom['coordinates'][0]
    y1 = geojs_geom['coordinates'][1]

    # transform the coordinate
    x, y = transform(in_crs, out_crs, x1, y1)

    # creat output new point
    new_point = dict(type='Feature', properties=dict(id=1))
    new_point['geometry'] = geojs_geom
    new_coord = (x, y)
    # add newly transformed coordinate
    new_point['geometry']['coordinates'] = new_coord

    return new_point

def proj_trans1(lon, lat):
    #p3 = pyproj.Proj("epsg:4326") 
    p1 = pyproj.Proj(proj='utm',zone=10,ellps='WGS84', preserve_units=False)
    p2 = pyproj.Proj("epsg:3857") 
    x1, y1 = p1(lon, lat)
    x2, y2 = pyproj.transform(p1, p2, x1, y1, radians=True)
    return x2, y2

def transform(self, data, src_srid, dst_srid):
        src = self.get_proj(src_srid)
        dst = self.get_proj(dst_srid)
        if src == dst:
            return data
        if data["type"] == "Point":
            x, y = data["coordinates"]
            data["coordinates"] = pyproj.transform(src, dst, x, y)
        elif data["type"] == "LineString":
            data["coordinates"] = [pyproj.transform(src, dst, x, y) for x, y in data["coordinates"]]
        return data

point2 = (point2_y, point2_x)

            # TODO improve by finding nearest edge, routing to/from node at either end
            origin_node = ox.get_nearest_node(G, point1)
            destination_node = ox.get_nearest_node(G, point2)

            try:
                if origin_node != destination_node:
                    # Find the shortest path over the network between these nodes
                    route = nx.shortest_path(G, origin_node, destination_node, weight='length')

                    # Retrieve route nodes and lookup geographical location
                    routeline = []
                    routeline.append((origin_x, origin_y))
                    for node in route:
                        routeline.append((transform(projWGS84, projUTM, G.nodes[node]['x'], G.nodes[node]['y'])))
                    routeline.append((dest_x, dest_y))
                    line = routeline
                else:
                    line = [(origin_x, origin_y), (dest_x, dest_y)]
            except nx.exception.NetworkXNoPath:
                line = [(origin_x, origin_y), (dest_x, dest_y)]

            # Map to line
            links.append({
                'type': "Feature",
                'geometry': {
                    "type": "LineString",
                    "coordinates": line
                },
                'properties': {
                    "origin": origin['properties']['id'],

if not os.access(outputDir, os.W_OK):
        raise IOError(errno.EACCES, "Not allowed to write to output directory %s" % (outputDir,))
    outputDir = os.path.abspath(outputDir)
    
    inRasterFilepath = os.path.join(outputDir, inRasterFilename)
    outRasterFilepath = os.path.join(outputDir, outRasterFilename)
    
    if not os.path.exists(outRasterFilepath):
        # Convert into coordinate system of inRaster
        inRasterSrs = getSpatialReferenceForRaster(inRasterFilepath)
        inSrs = osr.SpatialReference()
        assert(inSrs.ImportFromWkt(inRasterSrs[4]) == 0)
        p_in = Proj(init="EPSG:4326")
        p_out = Proj(inSrs.ExportToProj4())
        (ulX, ulY) = transform(p_in, p_out, bbox['minX'], bbox['maxY'])
        (lrX, lrY) = transform(p_in, p_out, bbox['maxX'], bbox['minY'])
        
        gdalCommand = "%s -q -stats -of GTiff -co 'COMPRESS=LZW' -projwin %f %f %f %f %s %s" \
                        % (gdalCmdPath, ulX, ulY, lrX, lrY, \
                           inRasterFilepath, outRasterFilepath)
        #print gdalCommand
        returnCode = os.system(gdalCommand)
        if returnCode != 0:
            raise Exception("GDAL command %s failed." % (gdalCommand,))

total = rad * 10800.0  # 3 hr rad to total rad
            else:
                total += rad * 10800.0
        now += interval

    ccursor.execute(
        """
        SELECT station, ST_x(geom), ST_y(geom), temp24_hour from
        alldata a JOIN stations t on
        (a.station = t.id) where day = %s
        """,
        (date.strftime("%Y-%m-%d"),),
    )
    for row in ccursor:
        if mode == NC_MODE:
            (x, y) = pyproj.transform(P4326, LCC, row[1], row[2])
            (gridxs, gridys, distances) = get_gp(xc, yc, x, y)
        else:
            # Note the hacky reversing of grid coords
            (gridys, gridxs, distances) = get_gp2(lons, lats, row[1], row[2])

        z0 = total[gridys[0], gridxs[0]]
        z1 = total[gridys[1], gridxs[1]]
        z2 = total[gridys[2], gridxs[2]]
        z3 = total[gridys[3], gridxs[3]]

        val = (
            z0 / distances[0]
            + z1 / distances[1]
            + z2 / distances[2]
            + z3 / distances[3]
        ) / (

def LLH_to_ECEF(self, lat, lon, alt):
        ecef, llh = pyproj.Proj(proj='geocent'), pyproj.Proj(proj='latlong')
        x, y, z = pyproj.transform(llh, ecef, lon, lat, alt, radians=False)
        return x/1000000, y/1000000, z/1000000