def getDistance(p1,p2,unit="m",p1_proj="EPSG:4326",p2_proj="EPSG:4326"):
    if p1_proj == 'aea':
        p1_aea = p1
    else:
        p1_aea = ops.transform(
            partial(
                pyproj.transform,
                pyproj.Proj(init=p1_proj),
                #pyproj.Proj(proj="aea",lat1=geometry.bounds[1],lat2=geometry.bounds[3])
                #use projection 'Albers Equal Conic Area for WA' to calcuate the area
                pyproj.Proj("+proj=aea +lat_1=-17.5 +lat_2=-31.5 +lat_0=0 +lon_0=121 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs ",lat1=p1.bounds[1],lat2=p1.bounds[3])
            ),
            p1
        )

    if p2_proj == 'aea':
        p2_aea = p2
    else:
        p2_aea = ops.transform(
            partial(
                pyproj.transform,
                pyproj.Proj(init=p2_proj),
                #pyproj.Proj(proj="aea",lat1=geometry.bounds[1],lat2=geometry.bounds[3])
                #use projection 'Albers Equal Conic Area for WA' to calcuate the area
                pyproj.Proj("+proj=aea +lat_1=-17.5 +lat_2=-31.5 +lat_0=0 +lon_0=121 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs ",lat1=p2.bounds[1],lat2=p2.bounds[3])
            ),
            p2
        )

    data = p1_aea.distance(p2_aea)
    if unit == "km" :
        return data / 1000.00 
    else:
        return data

    def update(self, transect):
        """
        Updates the container data to a profile that intersect the
        transect line.

        Returns nothing. Sets attributes as a side effect.

        Args:
            transect (LineString): A transect line.
        """
        Notice.info("Updating " + self.__class__.__name__)

        pad = self.settings['map_padding']
        aspect = 12/3.  # I was expecting it to be 8:3.

        # Calculate the map bounds and centre.
        bounds = transect.bounds
        llx, lly, urx, ury = bounds
        w, h = urx - llx, ury - lly

        x_adj, y_adj = 0, 0
        if h > (w/aspect):
            x_adj = ((aspect*h) - w) / 2.  # Aspect is hard-coded in uberplot
        else:
            y_adj = ((w/aspect) - h) / 2.

        utm_nad83 = pp.Proj("+init=EPSG:26920")
        ll_nad83 = pp.Proj("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")
        utm2lola = partial(pp.transform, utm_nad83, ll_nad83)

        ll = transform(utm2lola, Point(llx-pad-x_adj, lly-pad-y_adj))
        ur = transform(utm2lola, Point(urx+pad+x_adj, ury+pad+y_adj))
        self.ll, self.ur = ll, ur
        self.mid = Point(ll.x + 0.5*(ur.x-ll.x), ll.y + 0.5*(ur.y - ll.y))

        # Go over the layers and collect data.
        for layer, details in self.layers.items():
            path = details['file']
            print layer, path

            # Set up convenient params dictionary for plotting function.
            params = {k: v for k, v in details.items() if k != 'file'}
            self.layers[layer]['params'] = params

            # Get a list of shapes from the file.
            shapes = []
            fname, ext = os.path.splitext(os.path.basename(path))
            if ext.strip('.').lower() in self.settings['raster_extensions']:
                # TODO: Deal with rasters.
                pass
            elif ext.strip('.').lower() == 'shp':
                with fiona.open(path) as c:
                    for s in c:
                        shapes.append(shape(s['geometry']))
                        # name = s.get('name') or s.get('id') or None
                        # data = {name: shape(s['geometry'])}
                        # setattr(self, 'data', data)
                setattr(self, layer, shapes)
            else:
                pass

def getDistance(p1,p2,unit="m",p1_proj="EPSG:4326",p2_proj="EPSG:4326"):
    if p1_proj == 'aea':
        p1_aea = p1
    else:
        p1_aea = ops.transform(
            partial(
                pyproj.transform,
                pyproj.Proj(init=p1_proj),
                #pyproj.Proj(proj="aea",lat1=geometry.bounds[1],lat2=geometry.bounds[3])
                #use projection 'Albers Equal Conic Area for WA' to calcuate the area
                proj_aea(p1)
            ),
            p1
        )

    if p2_proj == 'aea':
        p2_aea = p2
    else:
        p2_aea = ops.transform(
            partial(
                pyproj.transform,
                pyproj.Proj(init=p2_proj),
                #pyproj.Proj(proj="aea",lat1=geometry.bounds[1],lat2=geometry.bounds[3])
                #use projection 'Albers Equal Conic Area for WA' to calcuate the area
                proj_aea(p2)
            ),
            p2
        )

    data = p1_aea.distance(p2_aea)
    if unit == "km" :
        return data / 1000.00 
    else:
        return data

    def test_drop_label(self):
        from shapely.ops import transform
        from tilequeue.tile import reproject_mercator_to_lnglat
        import math
        import dsl

        thresholds = {
            8:    200000000,
            9:    100000000,
            10:    10000000,
            11:     4000000,
            12:      750000,
            13:      100000,
            14:       50000,
            15:       10000,
        }

        for zoom in range(8, 16):
            area = thresholds.get(zoom)
            radius = math.sqrt(area / math.pi)

            coord = 2 ** (zoom - 1)

            # larger feature should retain name
            shape = dsl.tile_centre_shape(
                zoom, coord, coord).buffer(radius * 1.1)
            shape_lnglat = transform(
                    reproject_mercator_to_lnglat, shape)

            self.generate_fixtures(
                dsl.way(1, shape_lnglat, {
                    'natural': 'water',
                    'name': 'Foo',
                }),
            )

            self.assert_has_feature(
                zoom, coord, coord, 'water', {
                    'kind': 'water',
                    'name': 'Foo',
                })

            # smaller shape should drop it
            shape = dsl.tile_centre_shape(
                zoom, coord, coord).buffer(radius / 1.1)
            shape_lnglat = transform(
                    reproject_mercator_to_lnglat, shape)

            self.generate_fixtures(
                dsl.way(1, shape_lnglat, {
                    'natural': 'water',
                    'name': 'Foo',
                }),
            )

            self.assert_has_feature(
                zoom, coord, coord, 'water', {
                    'kind': 'water',
                    'name': type(None),
                })

def doPolygonize():
  blocks = polygonize(lines)
  writeBlocks(blocks, args[0] + '-blocks.geojson')

  blocks = polygonize(lines)
  bounds = Polygon([
    [minlng, minlat],
    [minlng, maxlat],
    [maxlng, maxlat],
    [maxlng, minlat],
    [minlng, minlat]
  ])
  # Geometry transform function based on pyproj.transform
  project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:3785'),
    pyproj.Proj(init='EPSG:4326'))
  print bounds
  print transform(project, bounds)

  print 'finding holes'
  for index, block in enumerate(blocks):
    if index % 1000 == 0:
      print "diff'd  %s" % (index)
    if not block.is_valid:
      print explain_validity(block)
      print transform(project, block)
    else:
      bounds = bounds.difference(block)
  print bounds

def fit_in_tile(z, x, y, shape):
    """
    Fit shape into the tile. Shape should be a Shapely geometry or WKT string
    with coordinates between 0 and 1. This unit square is then remapped into
    the tile z/x/y.
    """

    from ModestMaps.Core import Coordinate
    from shapely.ops import transform
    from shapely.wkt import loads as wkt_loads
    from tilequeue.tile import coord_to_mercator_bounds
    from tilequeue.tile import reproject_mercator_to_lnglat

    bounds = coord_to_mercator_bounds(Coordinate(zoom=z, column=x, row=y))

    if isinstance(shape, (str, unicode)):
        shape = wkt_loads(shape)

    # check shape fits within unit square, so we can transform it to fit
    # within the tile.
    assert shape.bounds[0] >= 0
    assert shape.bounds[1] >= 0
    assert shape.bounds[2] <= 1
    assert shape.bounds[3] <= 1

    def _transform(x, y, *unused_coords):
        return (
            x * (bounds[2] - bounds[0]) + bounds[0],
            y * (bounds[3] - bounds[1]) + bounds[1],
        )

    merc_shape = transform(_transform, shape)
    return transform(reproject_mercator_to_lnglat, merc_shape)

def apply_map_projection(polygon, meta_data, projection):
    m = Basemap(llcrnrlon=-50.,llcrnrlat=-50.,urcrnrlon=340.,urcrnrlat=65.,\
                resolution='h',projection=projection,
                lat_0=40.,lon_0=-20.,lat_ts=20.)
    polygon = ops.transform(m, polygon)
    meta_data.POINTS = [ops.transform(m, p) for p in meta_data.POINTS]
    meta_data.LONG = [p.coords[0][0] for p in meta_data.POINTS]
    meta_data.LAT = [p.coords[0][1] for p in meta_data.POINTS]

    return polygon, meta_data

def network_copy_offset(DG,distance=0.0):
    '''copies a network and shifts the subcatchment and stream coordinates by a distance
    
    used for testing conflation simpily using only one network
    Can be used to provide an understanding of the sensitivity of the network to positional error in features'''
    DG2 = DG.copy()
    if distance != 0:
        for f,t,k,data in DG2.edges_iter(data=True,keys=True):
            data['subCatch'] = transform(lambda x, y, z=None: (x+distance, y+distance), data['subCatch'])
            data['stream'] =  transform(lambda x, y, z=None: (x+distance, y+distance), data['stream'])
    return DG2

    def __getitem__(self, geometry):
        if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
            g = shape(geometry)
            if g.disjoint(shape(self)):
                raise ValueError("AOI does not intersect image: {} not in {}".format(g.bounds, self.bounds))
            bounds = ops.transform(self.__geo_transform__.rev, g).bounds
            result, xmin, ymin = self._slice_padded(bounds)
        else:
            if len(geometry) == 1:
                assert geometry[0] == Ellipsis
                return self

            elif len(geometry) == 2:
                arg0, arg1 = geometry
                if isinstance(arg1, slice):
                    assert arg0 == Ellipsis
                    return self[:, :, arg1.start:arg1.stop]
                elif arg1 == Ellipsis:
                    return self[arg0, :, :]

            elif len(geometry) == 3:
                try:
                    nbands, ysize, xsize = self.shape
                except:
                    ysize, xsize = self.shape
                band_idx, y_idx, x_idx = geometry
                if y_idx == Ellipsis:
                    y_idx = slice(0, ysize)
                if x_idx == Ellipsis:
                    x_idx = slice(0, xsize)
                if not(isinstance(y_idx, slice) and isinstance(x_idx, slice)):
                    di = DaskImage(self)
                    return di.__getitem__(geometry)
                xmin, ymin, xmax, ymax = x_idx.start, y_idx.start, x_idx.stop, y_idx.stop
                xmin = 0 if xmin is None else xmin
                ymin = 0 if ymin is None else ymin
                xmax = xsize if xmax is None else xmax
                ymax = ysize if ymax is None else ymax
                if ymin > ysize and xmin > xsize:
                    raise IndexError("Index completely out of image bounds")

                g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
                result = super(GeoDaskImage, self).__getitem__(geometry)

            else:
                return super(GeoDaskImage, self).__getitem__(geometry)

        gi = mapping(g)
        gt = self.__geo_transform__ + (xmin, ymin)
        image = super(GeoDaskImage, self.__class__).__new__(self.__class__, result, __geo_interface__ = gi, __geo_transform__ = gt)
        return image

def draw_preview_plot(polygon):
    ax, fig = new_fig()
    m = Basemap(llcrnrlon=-20., llcrnrlat=-50.,
                urcrnrlon=330., urcrnrlat=65.,
                rsphere=(6378137.00, 6356752.3142),
                resolution='h', projection='merc',
                lat_0=40., lon_0=-20., lat_ts=20.)
    m.drawcoastlines()
    m.fillcontinents()

    ops.transform(m, polygon)
    plot_mpp(ax, polygon, fc="red")

    finalize_plot(fig, ax)

def project(geom, projection1, projection2):
    """Reproject a shapely geometry object to new coordinate system

    Parameters
    ----------
    geom: shapely geometry object
    projection1: string
        Proj4 string specifying source projection
    projection2: string
        Proj4 string specifying destination projection
    """
    projection1 = str(projection1)
    projection2 = str(projection2)


    # define projections
    pr1 = pyproj.Proj(projection1, errcheck=True, preserve_units=True)
    pr2 = pyproj.Proj(projection2, errcheck=True, preserve_units=True)

    # projection function
    # (see http://toblerity.org/shapely/shapely.html#module-shapely.ops)
    project = partial(pyproj.transform, pr1, pr2)

    # do the transformation!
    return transform(project, geom)

def getClassBalance(pshapes,bounds,proj):
    """
    Get native class balance of projected shapes, assuming a rectangular bounding box.
    :param pshapes:
      Sequence of projected shapely shapes
    :param bounds:
      Desired bounding box, in decimal degrees.
    :param proj:
      PyProj object defining orthographic projection of shapes.
    :returns:
      Float fraction of hazard polygons (area of hazard polygons/total area of bbox) 
    """
    xmin,ymin,xmax,ymax = bounds
    bpoly = Polygon([(xmin,ymax),
                     (xmax,ymax),
                     (xmax,ymin),
                     (xmin,ymin)])
    project = partial(
        pyproj.transform,
        pyproj.Proj(proj='latlong', datum='WGS84'),
        proj)
    bpolyproj = transform(project,bpoly)
    totalarea = bpolyproj.area
    polyarea = 0
    for pshape in pshapes:
        polyarea += pshape.area

    return polyarea/totalarea

def area_from_lon_lat_poly(geometry):
    """
    Compute the area in km^2 of a shapely geometry, whose points are in longitude and latitude.

    Parameters
    ----------
    geometry: shapely geometry
        Points must be in longitude and latitude.

    Returns
    -------
    area:  float
        Area in km^2.

    """

    import pyproj
    from shapely.ops import transform
    from functools import partial


    project = partial(
        pyproj.transform,
        pyproj.Proj(init='epsg:4326'), # Source: Lon-Lat
        pyproj.Proj(proj='aea')) # Target: Albers Equal Area Conical https://en.wikipedia.org/wiki/Albers_projection

    new_geometry = transform(project, geometry)

    #default area is in m^2
    return new_geometry.area/1e6

 def test_multipolygon(self):
     g = geometry.MultiPoint([(0, 1), (0, 4)]).buffer(1.0)
     h = transform(lambda x, y, z=None: (x+1.0, y+1.0), g)
     self.assertEqual(h.geom_type, 'MultiPolygon')
     self.assertAlmostEqual(g.area, h.area)
     self.assertAlmostEqual(h.centroid.x, 1.0)
     self.assertAlmostEqual(h.centroid.y, 3.5)

 def test_polygon(self):
     g = geometry.Point(0, 1).buffer(1.0)
     h = transform(lambda x, y, z=None: (x+1.0, y+1.0), g)
     self.assertEqual(h.geom_type, 'Polygon')
     self.assertAlmostEqual(g.area, h.area)
     self.assertAlmostEqual(h.centroid.x, 1.0)
     self.assertAlmostEqual(h.centroid.y, 2.0)

def geo_geodetic_to_cartesian(geo_wkt):
    """
    Transforms a wkt string of geodetic coordinates to cartesian coordinates
    :param geo_wkt: wkt string of geodetic coordinates
    :returns : wkt string object of cartesian coordinates
    - - - Transformation information - - -
    'EPSG:4326': LatLon with WGS84 datum used by GPS units and Google Earth
    'aea': 'Albers Equal Area is a conic, equal area map projection that uses
    two standard parallels. Although scale and shape are not preserved,
    distortion is minimal between the standard parallels.
    """

    geo = loads(geo_wkt)
    lon_vec = geo.xy[0]
    lat_vec = geo.xy[1]
    geo_transformed = transform(
        partial(pyproj.transform,
                pyproj.Proj(init='EPSG:4326'),
                pyproj.Proj(proj='aea',
                            lon=lon_vec,
                            lat=lat_vec)
            ),
        geo)

    return geo_transformed.wkt

def projectdf(df, projection1, projection2):
    """Reproject a dataframe's geometry column to new coordinate system

    Parameters
    ----------
    df: dataframe
        Contains "geometry" column of shapely geometries

    projection1: string
        Proj4 string specifying source projection
    projection2: string
        Proj4 string specifying destination projection
    """
    projection1 = str(projection1)
    projection2 = str(projection2)


    # define projections
    pr1 = pyproj.Proj(projection1, errcheck=True, preserve_units=True)
    pr2 = pyproj.Proj(projection2, errcheck=True, preserve_units=True)

    # projection function
    # (see http://toblerity.org/shapely/shapely.html#module-shapely.ops)
    project = partial(pyproj.transform, pr1, pr2)

    # do the transformation!
    newgeo = [transform(project, g) for g in df.geometry]

    return newgeo

def get_record_buffers(road_srid, records_csv_uuid):
    """Returns a list of all records, buffered by the configured MATCH_TOLERANCE
    :param road_srid: EPSG ID of roads shapefile projection
    :param records_csv_uuid: UUID of the BlackSpotRecordsFile CSV
    """
    road_projection = {'init': 'epsg:{}'.format(road_srid)}
    record_projection = {'init': RECORD_PROJECTION}
    project = partial(
        pyproj.transform,
        pyproj.Proj(record_projection),
        pyproj.Proj(road_projection)
    )
    record_buffers = []
    with BlackSpotRecordsFile.objects.get(uuid=records_csv_uuid).csv as records_csv:
        records_csv.open('rb')
        try:
            for row in csv.DictReader(records_csv):
                try:
                    record_point = transform(project, Point(float(row[RECORD_COL_X]),
                                                            float(row[RECORD_COL_Y])))
                    record_buffers.append(record_point.buffer(MATCH_TOLERANCE))
                except RuntimeError:
                    continue
        finally:
            records_csv.close()

    return record_buffers

def convert_coordinates(x, y, transformation):
    """"""

    src_pt = Point(x, y)
    dst_pt = ops.transform(transformation, src_pt)

    return dst_pt.x, dst_pt.y

def transform(g, from_srs, to_srs):
    project = partial(
        pyproj.transform,
        parse_projection(from_srs),
        parse_projection(to_srs))

    return ops.transform(project, g)