Python shapely.geometry 模块，shape() 实例源码

def coords_for(name):
    geocoder = Nominatim()
    location = geocoder.geocode(name, geometry='geojson')

    try:
        geometry = shape(location.raw['geojson'])

        # Coordinates have to be flipped in order to work in overpass
        if geometry.geom_type == 'Polygon':
            west, south, east, north = geometry.bounds
            return BoundingBox(south, west, north, east)
        elif geometry.geom_type == 'MultiPolygon':
            bboxs = (BoundingBox(*(g.bounds[0:2][::-1] + g.bounds[2:][::-1]))
                     for g in geometry)
            return bboxs
        elif geometry.geom_type == 'Point':
            south, north, west, east = (float(coordinate)
                                        for coordinate in
                                        location.raw['boundingbox'])
            return BoundingBox(south, west, north, east)

    except (KeyError, AttributeError):
        raise AttributeError(
            'No bounding box available for this location name.')

def colorize_raster(masks):
    ''' (H, W, 10) -> (H, W, 3)
    '''
    assert masks.shape[2] == 10
    palette = np.array([(180, 180, 180), (100, 100, 100),  # Buildings, Misc.
                        (6, 88, 179), (125, 194, 223),  # Road, Track
                        (55, 120, 27), (160, 219, 166),  # Trees, Crops
                        (209, 173, 116), (180, 117, 69),  # Waterway, Standing
                        (67, 109, 244), (39, 48, 215)], dtype=np.uint8)  # Car

    r = []
    for obj_type in range(10):
        c = palette[obj_type]
        result = np.stack([masks[:, :, obj_type]] * 3, axis=2)
        r.append(result * c)
    r = np.stack(r)
    r = np.max(r, axis=0)
    return r

def jaccard_raster(true_raster, pred_raster):
    assert true_raster.shape[2] == pred_raster.shape[2] == 10
    score = []
    for i in range(10):
        true = true_raster[:, :, i] != 0
        pred = pred_raster[:, :, i] != 0
        tp = np.sum(true * pred)
        fp = np.sum(pred) - tp
        fn = np.sum(true) - tp
        if tp == 0:
            jac = 0.
        else:
            jac = tp / float(fp + fn + tp)
        score.append((tp, fp, fn, jac))
    score = np.array(score)
    assert score.shape == (10, 4)
    return score

def query_location(self, query):
        def swap_coords(geojson):
            out = []
            for x in geojson:
                if isinstance(x, list):
                    out.append(swap_coords(x))
                else:
                    return geojson[1], geojson[0]
            return out

        nom = Nominatim()
        try:
            geo = nom.geocode(query=query, geometry='geojson', timeout=3).raw
            geojson = geo['geojson']
        except (AttributeError, KeyError):
            raise FailedQuery('Query for {} did not return results.'.format(query))
        self.log.info('Nominatim returned {} for {}'.format(geo['display_name'], query))
        geojson['coordinates'] = swap_coords(geojson['coordinates'])
        self.location = shape(geojson)

def rasterizeShapes(pshapes, geodict, all_touched=True):
    """
    Rastering a shape

    Args:
        pshapes: Sequence of orthographically projected shapes.
        goedict: Mapio geodictionary.
        all_touched: Turn pixel "on" if shape touches pixel, otherwise turn it
            on if the center of the pixel is contained within the shape. Note
            that the footprint of the shape is inflated and the amount of
            inflation depends on the pixel size if all_touched=True.

    Returns:
        Rasterio grid.
    """
    outshape = (geodict.ny, geodict.nx)
    txmin = geodict.xmin - geodict.dx / 2.0
    tymax = geodict.ymax + geodict.dy / 2.0
    transform = Affine.from_gdal(
        txmin, geodict.dx, 0.0, tymax, 0.0, -geodict.dy)
    img = rasterio.features.rasterize(
            pshapes, out_shape=outshape, fill=0.0, transform=transform,
            all_touched=all_touched, default_value=1.0)
    return img

def sampleYes(array, N):
    """Sample without replacement N points from an array of XY coordinates.

    Args:
        array: 2D numpy array of XY points.
        N: Integer number of points to sample without replacement from
            input array.

    Returns:
        Tuple of (sampled points, unsampled points).

    """

    # array is a Mx2 array of X,Y points
    m, n = array.shape
    allidx = np.arange(0, m)
    sampleidx = np.random.choice(allidx, size=N, replace=False)
    nosampleidx = np.setxor1d(allidx, sampleidx)
    sampled = array[sampleidx, :]
    notsampled = array[nosampleidx, :]
    return (sampled, notsampled)

def convert2Prob(gdict, inventory, mustContainCenter=False):
    """Convert inventory shapefile to binary grid (with geodict of gdict) with 1 if any part of landslide was in a cell, 0 if not

    :param gdict: geodict, likely taken from model to compare inventory against
    :param inventory: full file path to shapefile of inventory, must be in geographic coordinates, WGS84
    :type inventory: string
    :returns rast: Grid2D object containing rasterized version of inventory where cell is 1. if any part of a landslide was in the cell
    """
    with fiona.open(inventory) as f:
        invshp = list(f.items())

    shapes = [shape(inv[1]['geometry']) for inv in invshp]

    # Rasterize with allTouch
    rast = Grid2D.rasterizeFromGeometry(shapes, gdict, fillValue=0., burnValue=1.0, mustContainCenter=mustContainCenter)

    return rast

def getProjectedShapes(shapes, xmin, xmax, ymin, ymax):
    """
    Take a sequence of geographic shapes and project them to a bounds-centered
    orthographic projection.

    Args:
        shapes: Sequence of shapes, as read in by fiona.collection().
        xmin: Eastern boundary of all shapes.
        xmax: Western boundary of all shapes.
        ymin: Southern boundary of all shapes.
        ymax: Northern boundary of all shapes.

    Returns:
       Tuple of
           - Input sequence of shapes, projected to orthographic
           - PyProj projection object used to transform input shapes
    """
    latmiddle = ymin + (ymax - ymin) / 2.0
    lonmiddle = xmin + (xmax - xmin) / 2.0
    projstr = ('+proj=ortho +datum=WGS84 +lat_0=%.4f +lon_0=%.4f '
               '+x_0=0.0 +y_0=0.0' % (latmiddle, lonmiddle))
    proj = pyproj.Proj(projparams=projstr)
    project = partial(
        pyproj.transform,
        pyproj.Proj(proj='latlong', datum='WGS84'),
        proj)

    pshapes = []
    for tshape in shapes:
        if tshape['geometry']['type'] == 'Polygon':
            pshapegeo = shape(tshape['geometry'])
        else:
            pshapegeo = shape(tshape['geometry'])
        pshape = transform(project, pshapegeo)
        pshapes.append(pshape)  # assuming here that these are simple polygons

    return (pshapes, proj)

def get_final_value(self, value, as_json=False):
        json_value = format_geojson(mapping(value))
        rounded_value = shape(json_value)

        shapely_logger.setLevel('ERROR')
        if rounded_value.is_valid:
            return json_value if as_json else rounded_value
        shapely_logger.setLevel('INFO')

        rounded_value = rounded_value.buffer(0)
        if not rounded_value.is_empty:
            value = rounded_value
        else:
            logging.debug('Fixing rounded geometry failed, saving it to the database without rounding.')

        return format_geojson(mapping(value), round=False) if as_json else value

def annotate_polygons(df, df_filter, ax, font_size, annotater):
        df[df_filter].apply(
            lambda x: ax.annotate(
                s=annotater(x),
                xy=x['shape'].centroid.coords[0],
                ha='center',
                va='center',
                color=Choroplether.GREY,
                fontsize=font_size,
                path_effects=[
                    patheffects.withStroke(linewidth=1, foreground="w")
                ]
            ),
            axis=1
        )
        return ax

def load_into_index(self, t, level, tile_dir):

    if tile_hierarchy.levels.has_key(level):
      file_name = tile_hierarchy.levels[level].GetFilename(t, level, os.environ['TILE_DIR'])

      if (file_name is not None and os.path.isfile(file_name)):
        # if the file has not be cached, we must load it up into the index.
        if t not in cached_tiles:

          with open(file_name) as f:
            geojson = json.load(f)

          for feature in geojson['features']:
            geom = shape(feature['geometry'])
            osmlr_id = long(feature['properties']['osmlr_id'])
            index.insert(osmlr_id, geom.bounds)

          # this is our set of tiles that have been loaded into the index, only load each tile once.
          cached_tiles.add(t)

  #parse the request because we dont get this for free!

def drawBuildings(polygons, displayRatio=1.0):
    fig = plt.figure()
    ax = fig.add_subplot(111)   
    for feature in polygons['features']:
        if displayRatio < 1.0:
            if random.random() >= displayRatio:
                continue
        polygon = shape(feature['geometry'])
        patch = PolygonPatch(polygon, facecolor='#FD7400', edgecolor='#FD7400', alpha=0.5, zorder=1)
        ax.add_patch(patch)

    bounds = getBuildingsBoundaries(polygons)
    minx, miny, maxx, maxy = bounds
    ax.set_xlim(minx,maxx)
    ax.set_ylim(miny,maxy)

    plt.show()

def update(self, bounds, data):


        ileft = int(round(np.floor((bounds[0] - self.bounds[0]) * 100) / (self.pixel_size * 100)))
        itop = int(round(np.floor((self.bounds[3] - bounds[3]) * 100) / (self.pixel_size * 100)))

        iright = ileft + data.shape[1]
        ibottom = itop + data.shape[0]

        try:

            # add worker surf as slice to sum_mean_surf
            self.grid[itop:ibottom, ileft:iright] += data

        except:
            print "#####"
            print "master: shape ({0}) grid shape ({1}) bounds ({2})".format(self.shape, self.grid.shape, self.bounds)
            print "subset: shape({0} left,right,top,bot ({1})".format(self.grid[itop:ibottom, ileft:iright].shape, (ileft, iright, itop, ibottom))
            print "data: shape ({0}) bounds ({1})".format(data.shape, bounds)
            print "#####"

            raise

def get_shape_within(self, shp, polys):
        """Find shape in set of shapes which another given shape is within.

        Args:
            shp (shape): shape object
            polys (List[shape]): list of shapes
        Returns:
            If shape is found in polys which shp is within, return shape.
            If not shape is found, return None.
        """
        if not hasattr(shp, 'geom_type'):
            raise Exception("CoreMSR [get_shape_within] : invalid shp given")

        if not isinstance(polys, list):
            raise Exception("CoreMSR [get_shape_within] : invalid polys given")

        for poly in polys:
            tmp_poly = shape(poly)
            if shp.within(tmp_poly):
                return tmp_poly

        print "CoreMSR [get_shape_within] : shp not within any given poly"

        return "None"

def is_in_grid(self, shp):
        """Check if arbitrary polygon is within grid bounding box.

        Args:
            shp (shape):
        Returns:
            Bool whether shp is in grid box.

        Depends on self.grid_box (shapely prep type) being defined
        in environment.
        """
        if not hasattr(shp, 'geom_type'):
            raise Exception("CoreMSR [is_in_grid] : invalid shp given")

        if not isinstance(self.grid_box, type(prep(Point(0,0)))):
            raise Exception("CoreMSR [is_in_grid] : invalid prep_adm0 "
                            "found")

        return self.grid_box.contains(shp)

def limit_geom_chars(geom, limit=8000000, step=0.0001):
    """Limit chars in geom geojson string by simplification

    Default limit for character count is 8000000
        (8MB - MongoDB max document size is 16MB)
    Default step for simplication is 0.0001
    """
    geom = shape(geom)
    curr_geom = geom
    ix = 0
    while len(str(curr_geom)) > limit:
        ix = ix + 1
        curr_geom = geom.simplify(step * ix)

    simplify_tolerance = step * ix
    return curr_geom.__geo_interface__, simplify_tolerance

def tmp_master_process(self, worker_result):
    task, geom, surf, bounds = worker_result

    if geom != "None":

        active_data.set_value(task, 'geom_val', geom)

        master_grid.update(bounds, surf)
        # ileft = (bounds[0] - core.bounds[0]) / core.pixel_size
        # itop = (core.bounds[3] - bounds[3]) / core.pixel_size

        # iright = ileft + surf.shape[0]
        # ibottom = itop + surf.shape[1]


        # # add worker surf as slice to sum_mean_surf
        # sum_mean_surf[ileft:iright, itop:ibottom] += surf


        # mstack.append_stack(surf)

        # if mstack.get_stack_size() > 1:
           # print "reducing stack"
        #    mstack.reduce_stack()

def load_nodes(self, nodes_file, id_field):

        # Now we load the layer as geographic objects for the actual analysis
        source = fiona.open(nodes_file, 'r')

        azims = []
        dirs = []
        ids = []
        self.nodes = []
        print 'Creating nodes spatial index'
        for feature in source:
            geom = shape(feature['geometry'])
            i_d = int(feature["properties"][id_field])
            self.idx_nodes.insert(i_d, geom.bounds)
            x, y = geom.centroid.coords.xy
            self.nodes.append([i_d, float(x[0]), float(y[0])])
        self.nodes = pd.DataFrame(self.nodes, columns = ['ID', 'X', 'Y']).set_index(['ID'])
        print '     NODES spatial index created successfully'

def setArray(self, incomingArray, copy=False):
        """
        You can use the self.array directly but if you want to copy from one array
        into a raster we suggest you do it this way
        :param incomingArray:
        :return:
        """
        masked = isinstance(self.array, np.ma.MaskedArray)
        if copy:
            if masked:
                self.array = np.ma.copy(incomingArray)
            else:
                self.array = np.ma.masked_invalid(incomingArray, copy=True)
        else:
            if masked:
                self.array = incomingArray
            else:
                self.array = np.ma.masked_invalid(incomingArray)

        self.rows = self.array.shape[0]
        self.cols = self.array.shape[1]
        self.min = np.nanmin(self.array)
        self.max = np.nanmax(self.array)

def PrintArr(arr):
    """
    Print an ASCII representation of the array with an up-down flip if the
    the cell height is negative.

    Int this scenario:
        - '-' means masked
        - '_' means nodata
        - '#' means a number
        - '0' means 0
    :param arr:
    """
    print "\n"
    for row in range(arr.shape[0]):
        rowStr = ""
        for col in range(arr[row].shape[0]):
            rowStr += str(arr[row][col])
        print "{0}:: {1}".format(row, rowStr)
    print "\n"

def get_data_polygon(rasterfile):

    import rasterio
    from rasterio.features import shapes
    from shapely.geometry import shape

    r = Raster(rasterfile)
    array = np.array(r.array.mask*1, dtype=np.int16)

#    with rasterio.drivers():
    with rasterio.open(rasterfile) as src:
        image = src.read(1)
        mask_array = image != src.nodata
        results = ({'properties': {'raster_val': v}, 'geometry': s}
                    for i, (s, v) in enumerate(shapes(array, mask=mask_array, transform=src.affine)))

    geoms = list(results)
    polygons = [shape(geom['geometry']) for geom in geoms]

    return polygons

def read_roads(roads_shp, records, buffer_size):
    """Reads shapefile and extracts roads and projection
    :param roads_shp: Path to the shapefile containing roads
    :param records: List of shapely geometries representing record points
    :param buffer_size: Number of units to buffer record for checking if road should be kept
    """
    # Create a spatial index for record buffers to efficiently find intersecting roads
    record_buffers_index = rtree.index.Index()
    for idx, record in enumerate(records):
        record_point = record['point']
        record_buffer_bounds = record_point.buffer(buffer_size).bounds
        record_buffers_index.insert(idx, record_buffer_bounds)

    shp_file = fiona.open(roads_shp)
    roads = []
    logger.info('Number of total roads in shapefile: {:,}'.format(len(shp_file)))
    for road in shp_file:
        road_shp = shape(road['geometry'])
        if should_keep_road(road, road_shp, record_buffers_index):
            roads.append(road_shp)

    return (roads, shp_file.bounds)

def is_point_on_land(coords, shape_file=None):
    """Checks if a coords are over land or over water. This is done useing
    a shape file of world boundaries and looping over all Polygons to see
    if the point is in any of it.

    """

    import fiona
    from shapely.geometry import Point, shape

    if not shape_file:
        shape_file='/home/data/mapdata/other/tm_world/TM_WORLD_BORDERS-0.3.shp'

    lon, lat=coords
    pt=Point(lon, lat)
    try:
        fc=fiona.open(shape_file)
    except:
        pass
    result=False
    for feature in fc:
        if shape(feature['geometry']).contains(pt):
          result=True
    fc.close()
    return result

def __getattribute__(self, name):
        fn = object.__getattribute__(self, name)
        if(isinstance(fn, types.MethodType) and
           any(name in C.__dict__ for C in self.__class__.__mro__)):
            @wraps(fn)
            def wrapped(*args, **kwargs):
                result = fn(*args, **kwargs)
                if isinstance(result, da.Array) and len(result.shape) in [2,3]:
                    copy = super(DaskImage, self.__class__).__new__(self.__class__,
                                                                    result.dask, result.name, result.chunks,
                                                                    result.dtype, result.shape)
                    copy.__dict__.update(self.__dict__)
                    try:
                        copy.__dict__.update(result.__dict__)
                    except AttributeError:
                        # this means result was an object with __slots__
                        pass
                    return copy
                return result
            return wrapped
        else:
            return fn

def _parse_geoms(self, **kwargs):
        """ Finds supported geometry types, parses them and returns the bbox """
        bbox = kwargs.get('bbox', None)
        wkt_geom = kwargs.get('wkt', None)
        geojson = kwargs.get('geojson', None)
        if bbox is not None:
            g = box(*bbox)
        elif wkt_geom is not None:
            g = wkt.loads(wkt_geom)
        elif geojson is not None:
            g = shape(geojson)
        else:
            return None
        if self.proj is None:
            return g
        else:
            return self._reproject(g, from_proj=kwargs.get('from_proj', 'EPSG:4326'))

def test_coords_input(self):
        geoms = []
        with fiona.open(self.shp) as source:
            for r in source:
                s = shape(r['geometry'])
                geoms.append(s)
        s2g.ShapeGraph(geoms[0:30], to_graph=True, resolution=0.01,
                       properties=['osm_id'])

def setUp(self):
        shp = os.path.join(os.path.dirname(__file__), '../data/campus.shp')
        self.geoms = []
        with fiona.open(shp) as source:
            for r in source:
                s = shape(r['geometry'])
                self.geoms.append(s)

def get_spectral_data(img_id, h, w, bands=['A', 'M', 'P']):
    res = []
    for waveband in bands:
        image_path = '{}/{}_{}.tif'.format(sixteen_band_path, img_id, waveband)
        image = tiff.imread(image_path)
        if len(image.shape) == 2:  # for panchromatic band
            image.shape = (1,) + image.shape
        image = image.transpose((1, 2, 0))
        image = cv2.resize(image, (w, h), interpolation=cv2.INTER_LANCZOS4)
        if len(image.shape) == 2:  # for panchromatic band
            image.shape += (1,)
        res.append(image)
    image = np.concatenate(res, axis=2)
    image = image.astype(np.float32)
    return image

def polygonize(im):
    assert len(im.shape) == 2
    shapes = features.shapes(im)
    polygons = []
    for i, shape in enumerate(shapes):
#        if i % 10000 == 0:
#            print(i)
        if shape[1] == 0:
            continue
        polygons.append(geometry.shape(shape[0]))
#    polygons = [geometry.shape(shape[0]) for shape in shapes if shape[1] > 0]
    mp = geometry.MultiPolygon(polygons)
    return mp

def unsoft(max_idx, num_class=11):
    max_idx = max_idx.squeeze()
    assert len(max_idx.shape) == 2
    preds = np.zeros((num_class - 1, max_idx.shape[0], max_idx.shape[1]))
    for i in range(1, num_class):
        preds[i - 1] = (max_idx == i)
    preds = preds.astype(np.uint8)
    return preds

def polygonize_sk(mask, level):
    contours = measure.find_contours(mask, level)
    polys = []
    for contour in contours:
        if contour.shape[0] < 4:
            continue
        poly = Polygon(shell=contour[:, [1, 0]])
        polys.append(poly)
    polys = MultiPolygon(polys)
    return polys

def get_state_from_coordinates(coordnates):

    #coordinates = np.array([(34, -118), (40.7, -74)])
    sf = shapefile.Reader('./data/us_states_st99/st99_d00')

    #sf = shapefile.Reader("./data/states/cb_2015_us_state_20m")
    shapes = sf.shapes()
    #shapes[i].points
    fields = sf.fields
    records = sf.records()
    state_polygons = defaultdict(list)
    for i, record in enumerate(records):
        state = record[5]
        points = shapes[i].points
        poly = shape(shapes[i])
        state_polygons[state].append(poly)

    coor_state = {}
    for i in range(coordnates.shape[0]):
        lat, lon = coordnates[i]
        for state, polies in state_polygons.iteritems():
            for poly in polies:
                point = Point(lon, lat)
                if poly.contains(point):
                    coor_state[(lat, lon)] = state.lower()
    return coor_state

def subset_s2(path, file, aoi):
    # Read Data________________________________________________________________
    print("SUBSET: Read Product...")
    sentinel = ProductIO.readProduct(path+file)
    print("SUBSET: Done reading!")
    # Get Band Names and print info
    name = sentinel.getName()
    print("SUBSET: Image ID:        %s" % name)
    band_names = sentinel.getBandNames()
    print("SUBSET: Bands:       %s" % (list(band_names)))
    # Preprocessing ___________________________________________________________
    # Resampling
    parameters = HashMap()
    parameters.put('targetResolution', 10)
    print("SUBSET: resample target resolution: 10m")
    product_resample = snappy.GPF.createProduct('Resample', parameters, sentinel)
    # Geojson to wkt
    with open(aoi) as f:
            gjson = json.load(f)       
    for feature in gjson['features']:
        polygon = (feature['geometry'])
    str_poly = json.dumps(polygon)
    gjson_poly = geojson.loads(str_poly)
    poly_shp = shape(gjson_poly)
    wkt = poly_shp.wkt
    # Subset
    geom = WKTReader().read(wkt)
    op = SubsetOp()
    op.setSourceProduct(product_resample)
    op.setGeoRegion(geom)
    product_sub = op.getTargetProduct()            
    # Write Data_______________________________________________________            
    print("SUBSET: Writing subset.")
    subset = path + name + '_subset_'
    ProductIO.writeProduct(product_sub, subset, "BEAM-DIMAP")    
    print("SUBSET: Done and saved in %s" % path)

def getProjectedShapes(shapes, xmin, xmax, ymin, ymax):
    """
    Take a sequence of geographic shapes and project them to a bounds-centered
    orthographic projection.

    Args:
        shapes: Sequence of shapes, as read in by fiona.collection().
        xmin: Eastern boundary of all shapes.
        xmax: Western boundary of all shapes.
        ymin: Southern boundary of all shapes.
        ymax: Northern boundary of all shapes.

    Returns:
       Tuple of
           - Input sequence of shapes, projected to orthographic
           - PyProj projection object used to transform input shapes
    """
    latmiddle = ymin + (ymax - ymin) / 2.0
    lonmiddle = xmin + (xmax - xmin) / 2.0
    projstr = ('+proj=ortho +datum=WGS84 +lat_0=%.4f +lon_0=%.4f '
               '+x_0=0.0 +y_0=0.0' % (latmiddle, lonmiddle))
    proj = pyproj.Proj(projparams=projstr)
    project = partial(
        pyproj.transform,
        pyproj.Proj(proj='latlong', datum='WGS84'),
        proj)

    pshapes = []
    for tshape in shapes:
        if tshape['geometry']['type'] == 'Polygon':
            pshapegeo = shape(tshape['geometry'])
        else:
            pshapegeo = shape(tshape['geometry'])
        pshape = transform(project, pshapegeo)
        pshapes.append(pshape)  # assuming here that these are simple polygons

    return (pshapes, proj)

def sampleNoPoints(sampleimg, N, xvar, yvar):
    '''
    Sample from our "no" sample grid, where that grid contains 1s where no
    pixels should be sampled from, and 0s where they should not.

    Args:
        sampleimg: Grid of shape (len(yvar),len(xvar)) where 1's represent
            pixels from which "no" values can be sampled.
        N: Number of pixels to sample (without replacement from sampleimg.
        xvar: Numpy array of centers of columns of sampling grid.
        yvar: Numpy array of centers of rows of sampling grid.

    Returns:
        Tuple of
          - nopoints sequence of x,y tuples representing "no" samples.
          - Modified version of input sampleimg, with nopoints pixels set to 0.

    '''

    # get N points from sampleimg without replacement
    # avoid nosampleidx indices
    # return an sequence of X,Y tuples from those indices
    npixels = len(xvar) * len(yvar)
    allidx = np.arange(0, npixels)
    sampleimg = sampleimg.flatten()  # flatten out the input image
    sampleidx = np.random.choice(allidx[sampleimg == 1], size=N, replace=False)
    sampleidx.sort()
    sampleimg[sampleidx] = 0
    sampleimg.shape = (len(yvar), len(xvar))
    sampley, samplex = np.unravel_index(sampleidx, sampleimg.shape)
    xp = xvar[samplex]
    yp = yvar[sampley]
    nopoints = list(zip(xp, yp))
    return (nopoints, sampleimg, sampleidx)

def getProjectedPatch(polygon, m, edgecolor, facecolor, lw=1., zorder=10):
    polyjson = mapping(polygon)
    tlist = []
    for sequence in polyjson['coordinates']:
        lon, lat = list(zip(*sequence))
        x, y = m(lon, lat)
        tlist.append(tuple(zip(x, y)))
    polyjson['coordinates'] = tuple(tlist)
    ppolygon = shape(polyjson)
    patch = PolygonPatch(ppolygon, facecolor=facecolor, edgecolor=edgecolor,
                         zorder=zorder, linewidth=lw, fill=True, visible=True)
    return patch

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def from_db_value(self, value, expression, connection, context):
        if value is None:
            return value
        return shape(json.loads(value))

def to_python(self, value):
        if value is None or value == '':
            return None
        try:
            geometry = shape(json.loads(value))
        except Exception:
            raise ValidationError(_('Invalid GeoJSON.'))
        self._validate_geomtype(geometry)
        try:
            geometry = clean_geometry(geometry)
        except Exception:
            raise ValidationError(_('Could not clean geometry.'))
        self._validate_geomtype(geometry)
        return geometry

def calc_area_of_polygons(clons, clats):
    """
    Calculate the area of lat-lon polygons.

    We project sphere onto a flat surface using an equal area projection
    and then calculate the area of flat polygon.

    This is slow we should do some caching to avoid recomputing.
    """

    areas = np.zeros(clons.shape[1:])
    areas[:] = np.NAN

    for j in range(areas.shape[0]):
        for i in range(areas.shape[1]):

            lats = clats[:, j, i]
            lons = clons[:, j, i]

            lat_centre = lats[0] + abs(lats[2] - lats[1]) / 2
            lon_centre = lons[0] + abs(lons[1] - lons[0]) / 2

            pa = pyproj.Proj(proj_str.format(lat_centre, lon_centre))
            x, y = pa(lons, lats)

            cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
            areas[j, i] = shape(cop).area

    assert(np.sum(areas) is not np.NAN)
    assert(np.min(areas) > 0)

    return areas

def add_poly_patches(geo_df, basemap):
        geo_df['shape'] = Choroplether.make_shapes(geo_df, basemap)
        geo_df['patches'] = Choroplether.make_default_polygon_patches(
            geo_df['shape']
        )
        return geo_df

def make_shapes(geo_df, bmap):
        return geo_df['geometry'].map(lambda x: transform(bmap, shape(x)))

def rasterize_geom(geom, shape, affine, all_touched=False):
    """
    Parameters
    ----------
    geom: GeoJSON geometry
    shape: desired shape
    affine: desired transform
    all_touched: rasterization strategy

    Returns
    -------
    ndarray: boolean
    """
    geoms = [(geom, 1)]
    rv_array = features.rasterize(
        geoms,
        out_shape=shape,
        transform=affine,
        fill=0,
        dtype='uint8',
        all_touched=all_touched)

    return rv_array.astype(bool)


# https://stackoverflow.com/questions/8090229/
#   resize-with-averaging-or-rebin-a-numpy-2d-array/8090605#8090605

def rebin_sum(a, shape, dtype):
    sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]
    return a.reshape(sh).sum(-1, dtype=dtype).sum(1, dtype=dtype)

def rasterize_pctcover_geom(geom, shape, affine, scale=None, all_touched=False):
    """
    Parameters
    ----------
    geom: GeoJSON geometry
    shape: desired shape
    affine: desired transform
    scale: scale at which to generate percent cover estimate

    Returns
    -------
    ndarray: float32
    """
    if scale is None:
        scale = 10

    min_dtype = min_scalar_type(scale ** 2)

    new_affine = Affine(affine[0]/scale, 0, affine[2],
                        0, affine[4]/scale, affine[5])

    new_shape = (shape[0] * scale, shape[1] * scale)

    rv_array = rasterize_geom(geom, new_shape, new_affine, all_touched=all_touched)

    rv_array = rebin_sum(rv_array, shape, min_dtype)

    return rv_array.astype('float32') / (scale**2)

def iter_latlong(self, indices=None):
        """Iterate over dataset as Shapely geometries in WGS 84 CRS."""
        _ = partial(project, from_proj=self.crs, to_proj='')
        if indices is None:
            for index, feature in enumerate(self):
                yield (index, _(shape(feature['geometry'])))
        else:
            for index in indices:
                yield (index, _(shape(self[index]['geometry'])))