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

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

    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

    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))

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)

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

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

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)

	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

 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

 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

    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)

    def getHorizontalDistance(self, point):
        """
		Get horizontal distance (great circle distance, in km) between current point
		and provided point (point).
		"""
        g = Geod(ellps="sphere")
        forward_azimuth, back_azimuth, horizontal_distance = g.inv(
            self.longitude, self.latitude, point.longitude, point.latitude
        )
        return horizontal_distance * 1e-3  # 1e-3 is needed to convert from m to km

 def _distance(self, start, to):
     '''
     Return distance.
     '''
     q = Geod(ellps='WGS84')
     fa, ba, d =\
         q.inv(start['lon'],
               start['lat'],
               to['lon'],
               to['lat'])
     return d

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

    def cell_height(self):
        cell_height = np.full(self.shape, np.nan)
        geod = Geod(ellps='WGS84')
        for i in range(self.shape[0]):
            _az12, _az21, cell_length = geod.inv(
                self.longitude[i, 0],
                self.latitude[i, 0] + .25,
                self.longitude[i, 0],
                self.latitude[i, 0] - .25,
            )
            cell_height[i, :] = cell_length

        return cell_height

    def getPoint(self, horizontal_distance, vertical_distance, azimuth):
        """
		Get point with given horizontal, and vertical distances (in km, 
		vertical distance: positive-downward, negative-upward) 
		and azimuth (in degrees) from current point. 
		"""
        # TODO: check horizontal distance is positive
        g = Geod(ellps="sphere")
        longitude, latitude, back_azimuth = g.fwd(
            self.longitude, self.latitude, azimuth, horizontal_distance * 1e3
        )  # 1e3 is needed to convert from km to m
        depth = self.depth + vertical_distance
        return Point(longitude, latitude, depth)

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

    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

    def size (self):
        from pyproj import Geod

        g = Geod(ellps='WGS84')
        lon_min = self.lonw
        lon_max = self.lone
        lat_min = self.lats
        lat_max = self.latn
        lon = (lon_min+lon_max)/2
        lat = (lat_min+lat_max)/2

        sn = g.inv (lon, lat_min, lon, lat_max, radians=False)[2]
        we = g.inv (lon_min, lat, lon_max, lat, radians=False)[2]

        return (sn, we)