def centroid_planar(longlat_pairs):
    """
    Centroid by projecting onto a plane. Inaccurate over long distances.
    """
    #
    # projections prep
    from pyproj import Proj, transform
    p_geodetic = Proj(proj='latlong')
    p_flat = Proj('+init=EPSG:27700')  # uk OS

    #
    # convert
    N = len(longlat_pairs)

    cartesians = np.zeros([N, 2])  # N-high
    for indx, pair in enumerate(longlat_pairs):
        lng, lat = map(float, pair)
        x, y = transform(p_geodetic, p_flat, lng, lat)
        cartesians[indx,:] = (x, y)

    x, y = np.average(cartesians, axis=0)

    lng, lat = transform(p_flat, p_geodetic, x, y)

    return (lng, lat)

	def test_indices(self):
		"""docstring for test_indices"""
		srs = 'epsg:25832'
		bbox = self.bbox
		self.assertEqual(self.ts.indices(srs=srs, bbox=bbox, resolution=1), (0,4,0,8))
		# lower left
		bbox = Bbox(500000, 5000000, 500999, 5000487)
		self.assertEqual(self.ts.indices(srs=srs, bbox=bbox, resolution=1), (2,4,0,4))
		# single cell
		bbox = Bbox(500010, 5000900, 500020, 5000910)
		self.assertEqual(self.ts.indices(srs=srs, bbox=bbox, resolution=1), (0,1,0,1))
		# outside
		bbox = Bbox(200000, 1000000, 300000, 2000000)
		self.assertEqual(self.ts.indices(srs=srs, bbox=bbox, resolution=1), None)
		# trunc
		bbox = Bbox(200000, 1000000, 700000, 7000000)
		self.assertEqual(self.ts.indices(srs=srs, bbox=bbox, resolution=1), (0,4,0,8))
		# reproject
		srs = 'epsg:4326'		
		utm = proj.Proj(init='epsg:25832')
		geo = proj.Proj(init='epsg:4326')
		min_x, min_y = proj.transform(utm, geo, 500010, 5000900)
		max_x, max_y = proj.transform(utm, geo, 500020, 5000910)
		bbox = Bbox(min_x, min_y, max_x, max_y)
		print "Bbox: %f, %f, %f, %f, srs: %s" % (bbox.min_x, bbox.min_y, bbox.max_x, bbox.max_y, srs)
		self.assertEqual(self.ts.indices(srs=srs, bbox=bbox, resolution=1), (0,1,0,1))

def get_map_geo():
    #przychodza w EPSG:3857
    if request.method == "GET":
        print "Heeelo"
        p3857 = pyproj.Proj(init="EPSG:3857")
        p4326 = pyproj.Proj(init="EPSG:4326")


        _long_min = request.args.get('long_min', None)
        _lat_min = request.args.get('lat_min', None)
        #_geo_diff = request.args.get('geo_diff', geo_diff)
        _long_max = request.args.get('long_max', None)
        _lat_max = request.args.get('lat_max', None)
        logger.debug(_long_min, _lat_min,)
        _long_min, _lat_min = pyproj.transform(p3857, p4326, _long_min, _lat_min)
        #logger.debug(str(_long_min), str(_lat_min), str(_long_max), str(_lat_max))
        _long_max, _lat_max = pyproj.transform(p3857, p4326, _long_max, _lat_max)
        print "bede preparowal"
        _wms_params = prepare_params_to_request_4326(_lat_min, _long_min, _lat_max, _long_max)
        print "spreparowalem " + str(_wms_params)
        response = get_map_image_in_b64(server_topo, _wms_params, "topo.jpg")
        print "przygotowuje odpowiedz"
        return Response(status=200)
    else:
        return Response(status=400)

def get_mbb_and_srs(layerDict):
    """Deprecated, please use get_crs and get_mbb"""
    #This function is unsatisfactory, it should actually be two different functions, one for the MBB and one for the SRS. Let's do that.
    for name in layerDict:
        with fiona.open(layerDict[name]['file']) as s:
            if 'minx' not in locals():
                minx, miny, maxx, maxy = s.bounds
                CRS = s.crs
                #very unclean way to get properly formatted SRS
                vector = ogr.GetDriverByName('ESRI Shapefile')
                src_ds = vector.Open(layerDict[name]['file'])
                src_lyr = src_ds.GetLayer(0)
                srs = osr.SpatialReference()
                srs.ImportFromWkt(src_lyr.GetSpatialRef().__str__())
            elif s.crs == CRS:
                x1, y1, x2, y2 = s.bounds
                if x1 < minx: minx = x1
                if y1 < miny: miny = y1
                if x2 > maxx: maxx = x2
                if y2 > maxy: maxy = y2
            else: #raise NameError('Layer: ' + str(name) + ' is of wrong CRS: ' + str(s.crs) + 'desired CRS: ' + str(CRS))
                x1, y1, x2, y2 = s.bounds
                proj1=pyproj.Proj(s.crs)
                proj2=pyproj.Proj(CRS,preserve_units=True)
                x1, y1 = pyproj.transform(proj1,proj2,x1,y1)
                x2, y2 = pyproj.transform(proj1,proj2,x2,y2)
                if x1 < minx: minx = x1
                if y1 < miny: miny = y1
                if x2 > maxx: maxx = x2
                if y2 > maxy: maxy = y2
    return (minx,miny,maxx,maxy),srs

def get_extrema(data,**kwargs):
	transform = False
	for key,value in kwargs.iteritems():
		if key == 'transform':
			transform = value
	north = data.top
	south = data.bottom
	west = data.left
	east = data.right

	extrema = {'n':north,'s':south,'w':west,'e':east}

	# newcorners
	ur_point = [extrema['e'],extrema['n']]
	ll_point = [extrema['w'],extrema['s']]


	if transform == True:
		p1 = pyproj.Proj(init='epsg:'+str(3857))
		p2 = pyproj.Proj(init='epsg:'+str(4326))

	 	ur_point = pyproj.transform(p1,p2,ur_point[0],ur_point[1])
	 	ll_point = pyproj.transform(p1,p2,ll_point[0],ll_point[1])

	 	print ll_point
	extrema = {'n':ur_point[1],'s':ll_point[1],'w':ll_point[0],'e':ur_point[0]}
	print extrema
	return extrema

    def getBoundsByAspect(bbList,aspect,projSrc,projDest):
        ''' Return bounds adjusted for aspect ratio '''
        bb = bbList
        bl = pyproj.transform(projSrc,projDest,bb[0],bb[1])
        tr = pyproj.transform(projSrc,projDest,bb[2],bb[3])

        projW = tr[0]-bl[0]
        projH = tr[1]-bl[1]
        projcx = (tr[0]+bl[0])/2
        projcy = (tr[1]+bl[1])/2

        pbb=[0,0,0,0]

        if projW/projH > aspect:
            pbb[0] = bl[0]
            pbb[2] = tr[0]
            pbb[1] = projcy-projW*0.5/aspect
            pbb[3] = projcy+projW*0.5/aspect
        else:
            pbb[1] = bl[1]
            pbb[3] = tr[1]
            pbb[0] = projcx-projH*0.5*aspect
            pbb[2] = projcx+projH*0.5*aspect

        nbl = pyproj.transform(projDest,projSrc,pbb[0],pbb[1])
        ntr = pyproj.transform(projDest,projSrc,pbb[2],pbb[3])

        return [nbl[0],nbl[1],ntr[0],ntr[1]]

def tricontourf_canvas(topology, datasetnc, request):
    """
    topology - netcdf topology object
    dataset - netcdf dataset object
    request - original http request
    """
    import wms_handler

    xmin, ymin, xmax, ymax = wms_handler.get_bbox(request)

    #proj = get_pyproj(request)
    #lonmin, latmin = proj(xmin, ymin, inverse=True)
    #lonmax, latmax = porj(xmax, ymax, inverse=True)
    CRS = get_pyproj(request)
    lonmin, latmin = pyproj.transform(CRS, EPSG4326, xmin, ymin)
    lonmax, latmax = pyproj.transform(CRS, EPSG4326, xmax, ymax)
    #logger.info("lonmin, latmin: {0} {1}".format(lonmin, latmin))
    #logger.info("lonmax, latmax: {0} {1}".format(lonmax, latmax))


    #compute triangular subset
    lon = topology.nodes[:,0]
    lat = topology.nodes[:,1]
    latlon_sub_idx = get_lat_lon_subset_idx(topology.nodes[:,0],
                                            topology[:,1],
                                            lonmin,
                                            latmin,
                                            lonmax,
                                            latmax)

    nv_sub_idx = get_nv_subset_idx(topology.faces[:], sub_idx)

 def get_buscar_ruta(self, **kwargs):
     values = json.loads(kwargs['rutas_wp'])
     google = pyproj.Proj("+init=EPSG:3857")
     gps = pyproj.Proj("+init=EPSG:4326")
     start = pyproj.transform(gps, google, values['start']['lng'], values['start']['lat'])
     finish = pyproj.transform(gps, google, values['end']['lng'], values['end']['lat'])
     sql = """SELECT id
         FROM movilidad_sostenible_oferta o
         WHERE ST_DWithin(o.shape, ST_SetSRID(ST_MakePoint(%s, %s), 900913), 500) AND
         ST_DWithin(o.shape, ST_SetSRID(ST_MakePoint(%s, %s), 900913), 500)
     """ # Busca rutas que esten a 500 metros de los puntos origen y destino seleccionados
     params = (start[0], start[1], finish[0], finish[1])
     request.env.cr.execute(sql, params)
     res = request.env.cr.fetchall()
     #print values['start']['lng'], values['start']['lat']
     #print values['end']['lng'], values['end']['lat']
     #print sql
     #print params
     #print res
     rutas_ids = [ i[0] for i in res ]
     rutas_model = request.env['movilidad_sostenible.oferta']
     rutas = rutas_model.browse(rutas_ids)
     #print rutas
     values.update(kwargs=kwargs.items())
     return http.request.render('movilidad_sostenible.lista_rutas_ofertar', {
         'lista_ofertas': rutas,
     })

    def _get_projected_extent(self, src_srs, target_srs, edge_points=9):
        """
        Densifies the edges with edge_points points between corners, and projects all of them.
        Returns the outer bounds of the projected coords.
        Geographic latitudes must be bounded to  -85.0511 <= y <= 85.0511 prior to calling this or things will fail!
        """

        srcProj = Proj(init=src_srs) if src_srs.strip().lower().startswith('epsg:') else Proj(str(src_srs))
        targetProj = Proj(init=target_srs) if ':' in target_srs else Proj(str(target_srs))
        if edge_points < 2:  # cannot densify (calculations below will fail), just use corners
            x_values, y_values = transform(srcProj, targetProj, [self.xmin, self.xmax], [self.ymin, self.ymax])
            return x_values[0], y_values[0], x_values[1], y_values[1]

        samples = range(0, edge_points)
        x_diff, y_diff = self.get_dimensions()
        xstep = x_diff/(edge_points-1)
        ystep = y_diff/(edge_points-1)
        x_values = []
        y_values = []
        for i, j in product(samples, samples):
            x_values.append(self.xmin + xstep * i)
            y_values.append(self.ymin + ystep * j)
        # TODO: check for bidrectional consistency, as is done in ncserve BoundingBox.project() method
        x_values, y_values = transform(srcProj, targetProj, x_values, y_values)
        return min(x_values), min(y_values), max(x_values), max(y_values)

def get_kml_poly(trans):
    outProj = Proj(init='EPSG:4326')
    inProj = Proj(init='EPSG:3857')

    kml = dict()
    kml['lonUL'], kml['latUL'] = transform(inProj, outProj, trans['lonUL'], trans['latUL'])
    kml['lonUR'], kml['latUR'] = transform(inProj, outProj, trans['lonUR'], trans['latUR'])
    kml['lonLL'], kml['latLL'] = transform(inProj, outProj, trans['lonLL'], trans['latLL'])
    kml['lonLR'], kml['latLR'] = transform(inProj, outProj, trans['lonLR'], trans['latLR'])

    return """
    <Polygon>
      <extrude>1</extrude>
      <altitudeMode>clampToGround</altitudeMode>
      <outerBoundaryIs>
        <LinearRing>
          <coordinates>
            {lonUL},{latUL},0
            {lonUR},{latUR},0
            {lonLR},{latLR},0
            {lonLL},{latLL},0
            {lonUL},{latUL},0
          </coordinates>
        </LinearRing>
      </outerBoundaryIs>
    </Polygon>
    """.format(**kml)

    def extent(self):
        """Return extent of the image as a dict with keys minx, maxx,
        miny, maxy, in Google projection."""
        data = gdal.Open(self.filename)
        pixelsx = data.RasterXSize
        pixelsy = data.RasterYSize

        startx, dxx, dyx, starty, dxy, dyy = data.GetGeoTransform()

        endx = startx + (pixelsx - 1) * dxx
        endy = starty - (pixelsy - 1) * dyy

        startx, starty = transform(
            RD_PROJECTION, GOOGLE_PROJECTION, startx, starty)
        endx, endy = transform(
            RD_PROJECTION, GOOGLE_PROJECTION, endx, endy)

        if startx > endx:
            startx, endx = endx, startx
        if starty > endy:
            starty, endy = endy, starty

        return {
            'minx': startx,
            'maxx': endx,
            'miny': starty,
            'maxy': endy,
            }

def check_grid_consistency(grid1,grid2,proj1,proj2):
    ''' Given two 2D (mesh)grids and their respective projections,
        perform a number of checks to verify their consistency '''

    def print_diffs(xt,yt,x2,y2):
        print 'Total accumulated longitude mismatch:',(xt-x2).sum()
        print 'Total accumulated latitude mismatch:',(yt-y2).sum()
        print 'Average longitudinal grid mismatch:',(xt-x2).mean()
        print 'Average latitudinal mismatch:',(yt-y2).mean()
        print 'Maximum longitudinal grid mismatch:',(xt-x2).max(),np.unravel_index((xt-x2).argmax(),xt.shape)
        print 'Maximum latitudinal mismatch:',(yt-y2).max(),np.unravel_index((yt-y2).argmax(),yt.shape),'\n'
        return

    for grid in ['mass_grid','u_grid','v_grid']:
        print 'Performing transformation on the %s'%grid

        x1,y1 = grid1['mass_grid']
        x2,y2 = grid2['mass_grid']

        # Transform proj1 to proj 2
        print 'Transform proj1 to proj2 (probably LCC to lat/lon)'
        xt,yt = pyproj.transform(proj1,proj2,x1,y1)
        print_diffs(xt,yt,x2,y2)

        # Transform proj2 to proj 1
        print 'Transform proj2 to proj1 (probably lat/lon to LCC)'
        xt,yt = pyproj.transform(proj2,proj1,x2,y2)
        print_diffs(xt,yt,x1,y1)

    return

def reproject_cutline_gmerc(src_proj, points):
    if points:
        points = densify_linestring(points)
        points1 = zip(*pyproj.transform(src_proj, proj_gmerc, *zip(*points)))
        points2 = zip(*pyproj.transform(src_proj, proj_gmerc_180, *zip(*points)))
        return points1 if calc_area(points1) <= calc_area(points2) else points2
    return []

def get_proj_dict(dem_dataset, epsg):
    dem_transform = {'upper_left_x': dem_dataset.GetGeoTransform()[0],
                     'x_res_m' : dem_dataset.GetGeoTransform()[1],
                     'x_rotation': dem_dataset.GetGeoTransform()[2],
                     'upper_left_y': dem_dataset.GetGeoTransform()[3],
                     'y_rotation': dem_dataset.GetGeoTransform()[4],
                     'y_res_m': dem_dataset.GetGeoTransform()[5],
                     'n_cols': dem_dataset.RasterXSize,
                     'n_rows': dem_dataset.RasterYSize}
    
    dem_transform['y_res_m'] *= -1
    dem_transform['east_min'] = dem_transform['upper_left_x']

    dem_transform['east_max'] = (dem_transform['x_res_m'] 
                                 * dem_transform['n_cols']
                                 + dem_transform['east_min'])
    
    dem_transform['north_max'] = dem_transform['upper_left_y']
    
    dem_transform['north_min'] = (-1 * dem_transform['y_res_m'] 
                                  * dem_transform['n_rows']
                                  + dem_transform['north_max'])

    wgs84 = pj.Proj(init='epsg:4326')
    utm = pj.Proj(init='epsg:{}'.format(epsg))
    
    dem_transform['lon_min'], dem_transform['lat_min'] = pj.transform(utm, wgs84,
                                                    dem_transform['east_min'], 
                                                    dem_transform['north_min'])
    
    dem_transform['lon_max'], dem_transform['lat_max'] = pj.transform(utm, wgs84,
                                                    dem_transform['east_max'], 
                                                    dem_transform['north_max'])
    return dem_transform

def load_grid(filename, proj_lat_long, proj_meters):
    print "loading:", filename
    grid_params = {}
    execfile(filename, grid_params)
    lon_min = grid_params['lon_min']
    lon_max = grid_params['lon_max']
    lat_min = grid_params['lat_min']
    lat_max = grid_params['lat_max']

    # project from lat-long
    ( x_min, y_min ) = pyproj.transform( proj_lat_long, proj_meters, lon_min, lat_min )
    ( x_max, y_max ) = pyproj.transform( proj_lat_long, proj_meters, lon_max, lat_max )

    grid = np.zeros( (grid_params['num_z'],
                      grid_params['num_lat'],
                      grid_params['num_lon'])
                     , np.double)
    return(grid,
           x_min,
           x_max,
           y_min,
           y_max,
           grid_params['z_min'],
           grid_params['z_max'],
           )

    def saveKML(self, kmlFile):
        """
        Saves a KML template to use with google earth.  Assumes x/y coordinates
        are lat/long, and creates an overlay to display the heatmap within Google
        Earth.

        kmlFile ->  output filename for the KML.
        """
        if self.img is None:
            raise Exception("Must first run heatmap() to generate image file.")

        tilePath = os.path.splitext(kmlFile)[0] + ".png"
        self.img.save(tilePath)

        if self.override:
            ((west, south), (east, north)) = self.area
        else:
            ((west, south), (east, north)) = self._ranges()

        #convert overlay BBOX if required
        if use_pyproj and self.srcepsg is not None and self.srcepsg != 'EPSG:4326':
          source = pyproj.Proj(init=self.srcepsg)
          dest = pyproj.Proj(init='EPSG:4326')
          (east,south) = pyproj.transform(source,dest,east,south)
          (west,north) = pyproj.transform(source,dest,west,north)

        bytes = self.KML % (tilePath, north, south, east, west)
        fh = open(kmlFile, "w")
        fh.write(bytes)
        fh.close()

    def project(self, target_projection, edge_points=9):
        """
        target_projection must be a pyproj projection.
        Densifies the edges with edge_points points between corners, and projects all of them.
        Returns the outer bounds of the projected coords.

        Note: beware projection issues when going to projections that don't fully encapsulate the world domain
        (e.g., Web Mercator has singularities above and below latitudes of ~ 85).
        """

        assert self.projection and isinstance(target_projection, Proj)

        if target_projection.srs == self.projection.srs:
            return self.clone()

        if edge_points < 2:
            # use corners only
            x_values, y_values = transform(self.projection, target_projection, [self.xmin, self.xmax], [self.ymin, self.ymax])
            return BBox((x_values[0], y_values[0], x_values[1], y_values[1]), projection=target_projection)

        samples = range(0, edge_points)
        xstep = float(self.xmax-self.xmin)/(edge_points-1)
        ystep = float(self.ymax-self.ymin)/(edge_points-1)
        x_values = []
        y_values = []
        for i, j in product(samples, samples):
            x_values.append(self.xmin + xstep * i)
            y_values.append(self.ymin + ystep * j)
        # TODO: check for bidrectional consistency, as is done in ncserve BoundingBox.project() method
        x_values, y_values = transform(self.projection, target_projection, x_values, y_values)
        return BBox((min(x_values), min(y_values), max(x_values), max(y_values)), target_projection)

    def _limit(self, lon, lat, target_cs): # TODO: lat long boundaries are not rectangular
        data_proj = Proj("+init=EPSG:4326")  # WGS84
        target_proj = Proj(target_cs)

        # Find bounding box in ERA projection
        bbox = self.bounding_box
        bb_proj = transform(target_proj, data_proj, bbox[0], bbox[1])
        lon_min, lon_max = min(bb_proj[0]), max(bb_proj[0])
        lat_min, lat_max = min(bb_proj[1]), max(bb_proj[1])
        #print(lon_min,lon_max,lat_min,lat_max)

        # Limit data
        lon_upper = lon >= lon_min
        lon_lower = lon <= lon_max

        lat_upper = lat >= lat_min
        lat_lower = lat <= lat_max

        lon_inds = np.nonzero(lon_upper == lon_lower)[0]
        lat_inds = np.nonzero(lat_upper == lat_lower)[0]
        # Masks
        lon_mask = lon_upper == lon_lower
        lat_mask = lat_upper == lat_lower
        
        #print (lon_inds,lat_inds)
        #print (lon[lon_inds],lat[lat_inds])

        if lon_inds.size == 0:
            raise ERAInterimDataRepositoryError("Bounding box longitudes don't intersect with dataset.")
        if lat_inds.size == 0:
            raise ERAInterimDataRepositoryError("Bounding box latitudes don't intersect with dataset.")

        x, y = transform(data_proj, target_proj, *np.meshgrid(lon[lon_inds], lat[lat_inds]))

        return x, y, (lon_mask, lat_mask), (lon_inds, lat_inds)

def get_image_tile(raster, x, y, z):
    try:
        bound = bounds(x, y, z)

        with rasterio.open(raster) as src:
            x_res, y_res = src.transform[0], src.transform[4]
            p1 = Proj({'init': 'epsg:4326'})
            p2 = Proj(**src.crs)

            # project tile boundaries from lat/lng to source CRS
            tile_ul_proj = transform(p1, p2, bound.west, bound.north)
            tile_lr_proj = transform(p1, p2, bound.east, bound.south)
            # get origin point from the TIF
            tif_ul_proj = (src.bounds.left, src.bounds.top)

            # use the above information to calculate the pixel indices of the window
            top = int((tile_ul_proj[1] - tif_ul_proj[1]) / y_res)
            left = int((tile_ul_proj[0] - tif_ul_proj[0]) / x_res)
            bottom = int((tile_lr_proj[1] - tif_ul_proj[1]) / y_res)
            right = int((tile_lr_proj[0] - tif_ul_proj[0]) / x_res)

            window = ((top, bottom), (left, right))

            # read the first three bands (assumed RGB) of the TIF into an array
            data = np.empty(shape=(3, 256, 256)).astype(src.profile['dtype'])
            for k in (1, 2, 3):
                src.read(k, window=window, out=data[k - 1], boundless=True)

            return Image.fromarray(np.moveaxis(data, 0, -1), mode='RGB')
    except Exception as err:
        raise TileNotFoundError(err)

    def create_circle(self, source, dist=1000, n_segments=20):
        """
        Create a circle around source in a given distance
        and with the given number of segments

        Parameters
        ----------
        source : Point-instance

        dist : float
            the distance around the point

        n_segments : int
            the number of segments of the (nearly) circle

        """
        if source.x is None:
            source_x, source_y = transform(self.p1, self.p2,
                                           source.lon, source.lat)
        else:
            source_x, source_y = source.x, source.y
        angel = np.linspace(0, np.pi*2, n_segments)
        x = source_x + dist * np.cos(angel)
        y = source_y + dist * np.sin(angel)
        lon, lat = transform(self.p2, self.p1, x, y)
        destinations = np.vstack([lon, lat]).T

        return destinations