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

def bounded_segments(lines, bounding_box, cut_segment=True):
    """
    Extract the bounded segments from a list of lines
    :param lines: a list of LineString
    :param bounding_box: the bounding coordinates in (minx, miny, maxx, maxy)
           or Polygon instance
    :return: a list of bounded segments
    """
    if isinstance(bounding_box, Polygon):
        bbox = bounding_box
    else:
        bbox = box(bounding_box[0], bounding_box[1],
                   bounding_box[2], bounding_box[3])
    segments = []
    for line in lines:
        if line.intersects(bbox):
            if cut_segment:
                segments.append(line.intersection(bbox))
            else:
                segments.append(line)
    return segments

def subgraph_within_box(self, bounding_box):
        """
        Extract a subgraph bounded by a box.
        :param bounding_box: the bounding coordinates in
            (minx, miny, maxx, maxy) or a Polygon instance
        :return: a subgraph of nx.Graph
        """
        if isinstance(bounding_box, Polygon):
            bbox = bounding_box
        else:
            bbox = box(bounding_box[0], bounding_box[1],
                       bounding_box[2], bounding_box[3])
        nbunch = set()
        for edge in self.graph.edges():
            s, e = edge
            if bbox.intersects(LineString([self.node_xy[s], self.node_xy[e]])):
                nbunch.add(s)
                nbunch.add(e)
        return self.graph.subgraph(nbunch)

def get_geometry_cells(self, geometry, bounds=None):
        if bounds is None:
            bounds = self._get_geometry_bounds(geometry)
        minx, miny, maxx, maxy = bounds

        height, width = self.data.shape
        minx = max(minx, self.x)
        miny = max(miny, self.y)
        maxx = min(maxx, self.x + width)
        maxy = min(maxy, self.y + height)

        from shapely import prepared
        from shapely.geometry import box

        cells = np.zeros_like(self.data, dtype=np.bool)
        prep = prepared.prep(geometry)
        res = self.resolution
        for iy, y in enumerate(range(miny * res, maxy * res, res), start=miny - self.y):
            for ix, x in enumerate(range(minx * res, maxx * res, res), start=minx - self.x):
                if prep.intersects(box(x, y, x + res, y + res)):
                    cells[iy, ix] = True

        return cells

def __init__(self, width: int, height: int, xoff=0, yoff=0, zoff=0,
                 scale=1, buffer=0, background='#FFFFFF', center=True):
        self.width = width
        self.height = height
        self.minx = xoff
        self.miny = yoff
        self.base_z = zoff
        self.scale = scale
        self.orig_buffer = buffer
        self.buffer = int(math.ceil(buffer*self.scale))
        self.background = background

        self.maxx = self.minx + width / scale
        self.maxy = self.miny + height / scale

        # how many pixels around the image should be added and later cropped (otherwise rsvg does not blur correctly)
        self.buffer = int(math.ceil(buffer*self.scale))
        self.buffered_width = self.width + 2 * self.buffer
        self.buffered_height = self.height + 2 * self.buffer
        self.buffered_bbox = box(self.minx, self.miny, self.maxx, self.maxy).buffer(buffer, join_style=JOIN_STYLE.mitre)

        self.background_rgb = tuple(int(background[i:i + 2], 16)/255 for i in range(1, 6, 2))

def testDispatcherClassifierOneRule(self):
        # create polygons to test
        box1 = box(0, 0, 100, 100)
        box2 = box(0, 0, 10, 10)

        dispatcher = RuleBasedDispatcher([CatchAllRule()])
        dispatcher_classifier = DispatcherClassifier(dispatcher, [AreaClassifier(500)])
        # simple dispatch test
        cls, probability, dispatch, _ = dispatcher_classifier.dispatch_classify(None, box1)
        self.assertEqual(1, cls)
        self.assertEqual(1.0, probability)
        self.assertEqual(0, dispatch)
        classes, probas, dispatches, _ = dispatcher_classifier.dispatch_classify_batch(None, [box1, box2])
        self.assertEqual(1, classes[0])
        self.assertEqual(0, classes[1])
        self.assertEqual(1.0, probas[0])
        self.assertEqual(1.0, probas[1])
        self.assertEqual(0, dispatches[0])
        self.assertEqual(0, dispatches[1])

def draw_multisquare(image, position, size, color_out=255, color_in=255):
    """Draw a square with color 'color_out' and given size at a given position (x, y)
        Then draw four square of size (size/5) with color 'color_in' at:
            1) coord: (y + (size / 5), x + (size / 5))
            2) coord: (y + (size / 5), x + (3 * size / 5))
            3) coord: (y + (3 * size / 5), x + (size / 5))
            4) coord: (y + (3 * size / 5), x + (3 * size / 5))
    """
    x, y = position
    small_size = size / 5
    image = draw_poly(image, box(x, y, x + size, y + size), color=color_out)
    square1 = box(x + small_size, y + small_size, x + 2 * small_size, y + 2 * small_size)
    square2 = box(x + 3 * small_size, y + small_size, x + 4 * small_size, y + 2 * small_size)
    square3 = box(x + small_size, y + 3 * small_size, x + 2 * small_size, y + 4 * small_size)
    square4 = box(x + 3 * small_size, y + 3 * small_size, x + 4 * small_size, y + 4 * small_size)
    squares = [square1, square2, square3, square4]
    for square in squares:
        image = draw_poly(image, square, color=color_in)
    return image

def window_from_polygon(self, polygon, mask=False):
        """Build and return a window being the minimum bounding box around the passed polygon.
        At least a part of the polygon should fit the image

        Parameters
        ----------
        polygon: shapely.geometry.Polygon
            The polygon of which the minimum bounding window should be returned
        mask: boolean (optional, default: False)
            True for applying a mask to the image

        Returns
        -------
        window: ImageWindow
            The resulting image window

        Raises
        ------
        IndexError: if the polygon box offset is not inside the image
        """
        offset, width, height = Image.polygon_box(polygon)
        if not self._check_tile_offset(offset):
            raise IndexError("Offset {} is out of the image.".format(offset))
        return self.window(offset, width, height, polygon_mask=polygon if mask else None)

def tile_from_polygon(self, tile_builder, polygon, mask=False):
        """Build a tile being the minimum bounding box around the passed polygon

        Parameters
        ----------
        tile_builder: TileBuilder
            The builder for effectively building the tile
        polygon: shapely.geometry.Polygon
            The polygon of which the bounding tile should be returned
        mask: boolean (optional, default: False)
            True for applying the polygon as an alpha mask to the tile

        Returns
        -------
        tile: Tile
            The bounding tile

        Raises
        ------
        IndexError: if the offset is not inside the image
        TileExtractionException: if the tile cannot be extracted
        """
        offset, width, height = Image.polygon_box(polygon)
        return self.tile(tile_builder, offset, width, height, polygon_mask=polygon if mask else None)

def polygon_box(polygon):
        """From a shapely polygon, return the information about the polygon bounding box.
        These information are offset (x, y), width and height.

        Parameters
        ----------
        polygon: shapely.geometry.Polygon
            The polygon of which the bounding box should be computed

        Returns
        -------
        offset: tuple (int, int)
            The offset of the polygon bounding box
        width: int
            The bounding box width
        height
            The bounding box heigth
        """
        minx, miny, maxx, maxy = polygon.bounds
        fminx, fminy = int(math.floor(minx)), int(math.floor(miny))
        cmaxx, cmaxy = int(math.ceil(maxx)), int(math.ceil(maxy))
        offset = (fminx, fminy)
        width = cmaxx - fminx
        height = cmaxy - fminy
        return offset, width, height

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 __init__(self, access_token=os.environ.get("DG_MAPS_API_TOKEN"),
                 url="https://api.mapbox.com/v4/digitalglobe.nal0g75k/{z}/{x}/{y}.png",
                 zoom=22, bounds=None):
        self.zoom_level = zoom
        self._token = access_token
        self._name = "image-{}".format(str(uuid.uuid4()))
        self._url_template = url + "?access_token={token}"

        _first_tile = mercantile.Tile(z=self.zoom_level, x=0, y=0)
        _last_tile = mercantile.Tile(z=self.zoom_level, x=180, y=-85.05)
        g = box(*mercantile.xy_bounds(_first_tile)).union(box(*mercantile.xy_bounds(_last_tile)))

        self._full_bounds = g.bounds

        # TODO: populate rest of fields automatically
        self._tile_size = 256
        self._nbands = 3
        self._dtype = "uint8"
        self.bounds = self._expand_bounds(bounds)
        self._chunks = tuple([self._nbands] + [self._tile_size, self._tile_size])

def aoi(self, **kwargs):
        """ Subsets the Image by the given bounds

        kwargs:
            bbox: optional. A bounding box array [minx, miny, maxx, maxy]
            wkt: optional. A WKT geometry string
            geojson: optional. A GeoJSON geometry dictionary

        Returns:
            image (ndarray): an image instance
        """
        g = self._parse_geoms(**kwargs)
        if g is None:
            return self
        else:
            return self[g]

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 __getitem__(self, geometry):
        if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
            image = GeoImage.__getitem__(self, geometry)
            return image
        else:
            result = DaskImage.__getitem__(self, geometry)
            image = super(GeoDaskWrapper, self.__class__).__new__(self.__class__,
                                                            result.dask, result.name, result.chunks,
                                                            result.dtype, result.shape)

            if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
                xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
                xmin = 0 if xmin is None else xmin
                ymin = 0 if ymin is None else ymin
                xmax = self.shape[2] if xmax is None else xmax
                ymax = self.shape[1] if ymax is None else ymax

                g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
                image.__geo_interface__ = mapping(g)
                image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
            else:
                image.__geo_interface__ = self.__geo_interface__
                image.__geo_transform__ = self.__geo_transform__
            return image

def test_subgraph_within_box(self):
        bounding_box = box(121.428387, 31.027371, 121.430863, 31.030227)
        a = time.time()
        subgraph = self.sg.subgraph_within_box(bounding_box)
        print(time.time() - a)
        if self.show_plots:
            plt.figure()
            nx.draw(subgraph, pos=self.sg.node_xy, node_size=50)
            plt.show()

def test_lines_within_box(self):
        bounding_box = box(121.428387, 31.027371, 121.430863, 31.030227)
        lines = self.sg.lines_within_box(bounding_box)
        for line in lines:
            assert line.is_major

def box_overlay(a, b):
    """Checking overlay by bounds (minx, miny, maxx, maxy)
    """
    # for i, j in product([0, 2], [1, 3]):
    #     px, py = (b1[i], b1[j])
    #     if b2[0] <= px <= b2[2] and b2[1] <= py <= b2[3]:
    #         return True
    # return False
    bbox1 = box(a[0], a[1], a[2], a[3])
    bbox2 = box(b[0], b[1], b[2], b[3])
    return bbox1.intersects(bbox2)

def test_query1(self):
        intersection = box(8348915.46680623, 543988.943201519, 8348915.4669, 543988.943201520)
        queried = query(self.uri, self.layer_name, 11, intersection)

        self.assertEqual(queried.to_numpy_rdd().first()[0], SpatialKey(1450, 996))

def test_query3(self):
        intersection = box(8348915.46680623, 543988.943201519, 8348915.4669, 543988.943201520).wkb
        queried = query(self.uri, self.layer_name,
                        11, intersection)

        self.assertEqual(queried.to_numpy_rdd().first()[0], SpatialKey(1450, 996))

def test_query_partitions(self):
        intersection = box(8348915.46680623, 543988.943201519, 8348915.4669, 543988.943201520)
        queried = query(self.uri, self.layer_name,
                        11, intersection, num_partitions=2)

        self.assertEqual(queried.to_numpy_rdd().first()[0], SpatialKey(1450, 996))

def test_query_crs(self):
        intersection = box(74.99958369653905, 4.8808219582513095, 74.99958369738141, 4.880821958251324)
        queried = query(self.uri, self.layer_name, 11, intersection,
                        query_proj=4326)

        self.assertEqual(queried.to_numpy_rdd().first()[0], SpatialKey(1450, 996))

def to_polygon(self):
        """Converts this instance to a Shapely Polygon.

        The resulting Polygon will be in the shape of a box.

        Returns:
            ``shapely.geometry.Polygon``
        """

        return box(*self)

def test_ccw(self):
        b = geometry.box(0, 0, 1, 1, ccw=True)
        self.assertEqual(b.exterior.coords[0], (1.0, 0.0))
        self.assertEqual(b.exterior.coords[1], (1.0, 1.0))

def test_ccw_default(self):
        b = geometry.box(0, 0, 1, 1)
        self.assertEqual(b.exterior.coords[0], (1.0, 0.0))
        self.assertEqual(b.exterior.coords[1], (1.0, 1.0))

def test_cw(self):
        b = geometry.box(0, 0, 1, 1, ccw=False)
        self.assertEqual(b.exterior.coords[0], (0.0, 0.0))
        self.assertEqual(b.exterior.coords[1], (0.0, 1.0))

def parse_geo_box(geo_box_str):
    """
    parses [-90,-180 TO 90,180] to a shapely.geometry.box
    :param geo_box_str:
    :return:
    """

    from_point_str, to_point_str = parse_solr_geo_range_as_pair(geo_box_str)
    from_point = parse_lat_lon(from_point_str)
    to_point = parse_lat_lon(to_point_str)
    rectangle = box(from_point[0], from_point[1], to_point[0], to_point[1])
    return rectangle

def aa_bbox(self):
        """Return axis-aligned bounding box."""
        return self.fov.bounds

def __init__(self, aa_bbox):
        """Init."""
        super().__init__(box(aa_bbox))

def bbox(self):
        return box(self.minx-1, self.miny-1, self.maxx+1, self.maxy+1)

def testMerger(self):
        # build chunks for the polygons
        tile_box = box(0, 0, 512, 256)  # a box having the same dimension as the tile
        circle = Point(600, 360)
        circle = circle.buffer(250)

        circle_part1 = tile_box.intersection(circle)
        circle_part2 = translate(tile_box, xoff=512).intersection(circle)
        circle_part3 = translate(tile_box, yoff=256).intersection(circle)
        circle_part4 = translate(tile_box, xoff=512, yoff=256).intersection(circle)
        circle_part5 = translate(tile_box, yoff=512).intersection(circle)
        circle_part6 = translate(tile_box, xoff=512, yoff=512).intersection(circle)

        # create topology
        fake_image = FakeImage(1024, 768, 3)
        fake_builder = FakeTileBuilder()
        topology = fake_image.tile_topology(fake_builder, 512, 256)

        tile1 = topology.tile(1)
        tile2 = topology.tile(2)
        tile3 = topology.tile(3)
        tile4 = topology.tile(4)
        tile5 = topology.tile(5)
        tile6 = topology.tile(6)

        tiles = [tile1.identifier, tile2.identifier, tile3.identifier, tile4.identifier, tile5.identifier, tile6.identifier]
        tile_polygons = [[circle_part1], [circle_part2], [circle_part3], [circle_part4], [circle_part5], [circle_part6]]

        polygons = SemanticMerger(5).merge(tiles, tile_polygons, topology)
        self.assertEqual(len(polygons), 1, "Number of found polygon")

        # use recall and false discovery rate to evaluate the error on the surface
        tpr = circle.difference(polygons[0]).area / circle.area
        fdr = polygons[0].difference(circle).area / polygons[0].area
        self.assertLessEqual(tpr, 0.002, "Recall is low for circle area")
        self.assertLessEqual(fdr, 0.002, "False discovery rate is low for circle area")

def testRuleBasedDispatcherNoLabels(self):
        # prepare data for test
        box1 = box(0, 0, 100, 100)
        box2 = box(0, 0, 10, 10)

        dispatcher = RuleBasedDispatcher([CatchAllRule()])
        self.assertEqual(dispatcher.dispatch(None, box1), 0)
        dispatch_batch = dispatcher.dispatch_batch(None, [box1, box2])
        assert_array_equal(dispatch_batch, [0, 0])
        labels, dispatch_map = dispatcher.dispatch_map(None, [box1, box2])
        assert_array_equal(labels, dispatch_batch)
        assert_array_equal(dispatch_batch, dispatch_map)

def testRuleBasedDispatcher(self):
        # prepare data for test
        box1 = box(0, 0, 100, 100)
        box2 = box(0, 0, 10, 10)

        dispatcher = RuleBasedDispatcher([CatchAllRule()], ["catchall"])
        self.assertEqual(dispatcher.dispatch(None, box1), "catchall")
        dispatch_batch = dispatcher.dispatch_batch(None, [box1, box2])
        assert_array_equal(dispatch_batch, ["catchall", "catchall"])
        labels, dispatch_map = dispatcher.dispatch_map(None, [box1, box2])
        assert_array_equal(labels, dispatch_batch)
        assert_array_equal(dispatch_map, [0, 0])

def testCustomDispatcher(self):
        # prepare data for test
        box1 = box(0, 0, 500, 500)
        box2 = box(0, 0, 10, 10)
        box3 = box(0, 0, 1000, 1000)

        dispatcher = CustomDispatcher()
        self.assertEqual(dispatcher.dispatch(None, box1), "BIG")
        dispatch_batch = dispatcher.dispatch_batch(None, [box1, box2, box3])
        assert_array_equal(dispatch_batch, ["BIG", "SMALL", "BIG"])
        labels, dispatch_map = dispatcher.dispatch_map(None, [box1, box2, box3])
        assert_array_equal(labels, dispatch_batch)
        assert_array_equal(dispatch_map, [1, 0, 1])

def _init_polygon_mask(self, polygon_mask):
        """Sets the polygon mask taking into account parent and passed mask"""
        parent_pmask = self._parent_polygon_mask(do_translate=True)

        if polygon_mask is not None and parent_pmask is not None:
            self._polygon_mask = polygon_mask.intersection(parent_pmask)
        elif polygon_mask is not None:
            self._polygon_mask = polygon_mask
        elif parent_pmask is not None:
            self._polygon_mask = box(0, 0, self.width, self.height).intersection(parent_pmask)
        else:
            self._polygon_mask = None

def __init__(self, min_x, min_y, max_x, max_y):
     self.minx = min_x
     self.miny = min_y
     self.maxx = max_x
     self.maxy = max_y

  #The bounding boxes do NOT intersect if the other bounding box is
  #entirely LEFT, BELOW, RIGHT, or ABOVE bounding box.

def TileId(self, y, x):
    if (y < self.bbox.miny or x < self.bbox.minx or
        y > self.bbox.maxy or x > self.bbox.maxx):
      return -1

    #Find the tileid by finding the latitude row and longitude column
    return (self.Row(y) * self.ncolumns) + self.Col(x)

  # Get the bounding box of the specified tile.

def BottomNeighbor(self, tileid):
    return tileid if (tileid < self.ncolumns) else (tileid - self.ncolumns)

  # Get the list of tiles that lie within the specified bounding box.
  # The method finds the center tile and spirals out by finding neighbors
  # and recursively checking if tile is inside and checking/adding
  # neighboring tiles

def get_projects_geojson(search_bbox_dto: ProjectSearchBBoxDTO) -> geojson.FeatureCollection:
        """  search for projects meeting criteria provided return as a geojson feature collection"""

        # make a polygon from provided bounding box
        polygon = ProjectSearchService._make_4326_polygon_from_bbox(search_bbox_dto.bbox, search_bbox_dto.input_srid)

        # validate the bbox area is less than or equal to the max area allowed to prevent
        # abuse of the api or performance issues from large requests
        if not ProjectSearchService.validate_bbox_area(polygon):
            raise BBoxTooBigError('Requested bounding box is too large')

        # get projects intersecting the polygon for created by the author_id
        intersecting_projects = ProjectSearchService._get_intersecting_projects(polygon, search_bbox_dto.project_author)

        # allow an empty feature collection to be returned if no intersecting features found, since this is primarily
        # for returning data to show on a map
        features = []
        for project in intersecting_projects:
            try:
                localDTO = ProjectInfo.get_dto_for_locale(project.id, search_bbox_dto.preferred_locale,
                                                          project.default_locale)
            except Exception as e:
                pass

            properties = {
                "projectId": project.id,
                "projectStatus": ProjectStatus(project.status).name,
                "projectName": localDTO.name
            }
            feature = geojson.Feature(geometry=geojson.loads(project.geometry), properties=properties)
            features.append(feature)

        return geojson.FeatureCollection(features)

def _make_4326_polygon_from_bbox(bbox: list, srid: int) -> Polygon:
        """ make a shapely Polygon in SRID 4326 from bbox and srid"""
        try:
            polygon = box(bbox[0], bbox[1], bbox[2], bbox[3])
            if not srid == 4326:
                geometry = shape.from_shape(polygon, srid)
                geom_4326 = db.engine.execute(ST_Transform(geometry, 4326)).scalar()
                polygon = shape.to_shape(geom_4326)
        except Exception as e:
            raise ProjectSearchServiceError(f'error making polygon: {e}')
        return polygon

def initialize_grid(self, geom):
        """
        """

        bounds = geom.bounds

        (minx, miny, maxx, maxy) = bounds

        (minx, miny, maxx, maxy) = (
            round(np.floor(minx * self.psi)) / self.psi,
            round(np.floor(miny * self.psi)) / self.psi,
            round(np.ceil(maxx * self.psi)) / self.psi,
            round(np.ceil(maxy * self.psi)) / self.psi)

        clean_bounds = (minx, miny, maxx, maxy)

        top_left_lon = minx
        top_left_lat = maxy
        affine = Affine(self.pixel_size, 0, top_left_lon,
                        0, -self.pixel_size, top_left_lat)


        # base_rasterize, base_bounds = self.rasterize_geom(geom)
        # self.shape = base_rasterize.shape

        nrows = int(np.ceil( (maxy - miny) / self.pixel_size ))
        ncols = int(np.ceil( (maxx - minx) / self.pixel_size ))
        self.shape = (nrows, ncols)


        self.bounds = clean_bounds
        self.affine = affine
        self.topleft = (top_left_lon, top_left_lat)
        self.grid_box = prep(box(*self.bounds))


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

def _expand_bounds(self, bounds):
        if bounds is None:
            return bounds
        min_tile_x, min_tile_y, max_tile_x, max_tile_y = self._tile_coords(bounds)

        ul = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=min_tile_x, y=max_tile_y)))
        lr = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=max_tile_x, y=min_tile_y)))

        return ul.union(lr).bounds

def _tile_coords(self, bounds):
        """ Convert tile coords mins/maxs to lng/lat bounds """
        tfm = partial(pyproj.transform,
                      pyproj.Proj(init="epsg:3857"),
                      pyproj.Proj(init="epsg:4326"))
        bounds = ops.transform(tfm, box(*bounds)).bounds
        params = list(bounds) + [[self.zoom_level]]
        tile_coords = [(tile.x, tile.y) for tile in mercantile.tiles(*params)]
        xtiles, ytiles = zip(*tile_coords)
        minx = min(xtiles)
        maxx = max(xtiles)
        miny = min(ytiles)
        maxy = max(ytiles)
        return minx, miny, maxx, maxy

def __new__(cls, access_token=os.environ.get("DG_MAPS_API_TOKEN"),
                url="https://api.mapbox.com/v4/digitalglobe.nal0g75k/{z}/{x}/{y}.png",
                zoom=22, **kwargs):
        _tms_meta = TmsMeta(access_token=access_token, url=url, zoom=zoom, bounds=kwargs.get("bounds"))
        self = super(TmsImage, cls).create(_tms_meta)
        self._base_args = {"access_token": access_token, "url": url, "zoom": zoom}
        self._tms_meta = _tms_meta
        self.__geo_interface__ = mapping(box(*_tms_meta.bounds))
        self.__geo_transform__ = _tms_meta.__geo_transform__
        g = self._parse_geoms(**kwargs)
        if g is not None:
            return self[g]
        else:
            return self

def __getitem__(self, geometry):
        if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
            if self._tms_meta._bounds is None:
                return self.aoi(geojson=mapping(geometry), from_proj=self.proj)
            image = GeoImage.__getitem__(self, geometry)
            image._tms_meta = self._tms_meta
            return image
        else:
            result = super(TmsImage, self).__getitem__(geometry)
            image = super(TmsImage, self.__class__).__new__(self.__class__,
                                                            result.dask, result.name, result.chunks,
                                                            result.dtype, result.shape)

            if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
                xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
                xmin = 0 if xmin is None else xmin
                ymin = 0 if ymin is None else ymin
                xmax = self.shape[2] if xmax is None else xmax
                ymax = self.shape[1] if ymax is None else ymax

                g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
                image.__geo_interface__ = mapping(g)
                image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
            else:
                image.__geo_interface__ = self.__geo_interface__
                image.__geo_transform__ = self.__geo_transform__
            image._tms_meta = self._tms_meta
            return image

def bounds(self):
        """ Access the spatial bounding box of the image 

        Returns:
            bounds (list): list of bounds in image projected coordinates (minx, miny, maxx, maxy)
        """
        return shape(self).bounds

def _transpix(self, geometry, gsd, dem, proj):
        xmin, ymin, xmax, ymax = geometry.bounds
        x = np.linspace(xmin, xmax, num=int((xmax-xmin)/gsd))
        y = np.linspace(ymax, ymin, num=int((ymax-ymin)/gsd))
        xv, yv = np.meshgrid(x, y, indexing='xy')

        if self.proj is None:
            from_proj = "EPSG:4326"
        else:
            from_proj = self.proj

        itfm = partial(pyproj.transform, pyproj.Proj(init=proj), pyproj.Proj(init=from_proj))

        xv, yv = itfm(xv, yv) # if that works

        if isinstance(dem, GeoImage):
            g = box(xv.min(), yv.min(), xv.max(), yv.max())
            try:
                dem = dem[g].compute(get=dask.get) # read(quiet=True)
            except AssertionError:
                dem = 0 # guessing this is indexing by a 0 width geometry.

        if isinstance(dem, np.ndarray):
            dem = tf.resize(np.squeeze(dem), xv.shape, preserve_range=True, order=1, mode="edge")

        return self.__geo_transform__.rev(xv, yv, z=dem, _type=np.float32)[::-1]

def __contains__(self, g):
        geometry = ops.transform(self.__geo_transform__.rev, g)
        img_bounds = box(0, 0, *self.shape[2:0:-1])
        return img_bounds.contains(geometry)

def _image_by_type(cls, cat_id, **kwargs):
        vectors = Vectors()
        aoi = wkt.dumps(box(-180, -90, 180, 90))
        query = "item_type:GBDXCatalogRecord AND attributes.catalogID:{}".format(cat_id)
        query += " AND NOT item_type:IDAHOImage AND NOT item_type:DigitalGlobeAcquisition"
        result = vectors.query(aoi, query=query, count=1)
        if len(result) == 0:
            raise Exception('Could not find a catalog entry for the given id: {}'.format(cat_id))
        else:
            return cls._image_class(cat_id, result[0], **kwargs)

def _build_standard_products(idaho_id, bbox, proj):
        wkt = box(*bbox).wkt
        dem = ipe.GeospatialCrop(ipe.IdahoRead(bucketName="idaho-dems", imageId=idaho_id, objectStore="S3"), geospatialWKT=str(wkt))
        if proj is not "EPSG:4326":
            dem = ipe.Reproject(dem, **reproject_params(proj))
        return {
            "dem": dem
        }

def __getitem__(self, geometry):
        if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
            image = GeoImage.__getitem__(self, geometry)
            image._ipe_op = self._ipe_op
            return image
        else:
            result = super(IpeImage, self).__getitem__(geometry)
            image = super(IpeImage, self.__class__).__new__(self.__class__,
                                                            result.dask, result.name, result.chunks,
                                                            result.dtype, result.shape)

            if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
                xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
                xmin = 0 if xmin is None else xmin
                ymin = 0 if ymin is None else ymin
                xmax = self.shape[2] if xmax is None else xmax
                ymax = self.shape[1] if ymax is None else ymax

                g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
                image.__geo_interface__ = mapping(g)
                image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
            else:
                image.__geo_interface__ = self.__geo_interface__
                image.__geo_transform__ = self.__geo_transform__
            image._ipe_op = self._ipe_op
            return image