[docs]classGeoContext(object):""" Specifies spatial parameters to use when loading a raster from the Descartes Labs catalog. Two Images loaded with the same GeoContext will result in images with the same shape (in pixels), covering the same spatial extent, regardless of the dimensions or projection of the original data. Specifically, a fully-defined GeoContext specifies: * geometry to use as a cutline (WGS84), and/or bounds * resolution (m) or a shape defining the extent in pixels * EPSG code of the output coordinate reference system * whether to align pixels to the output CRS (see docstring for `AOI.align_pixels` for more information) GeoContexts are immutable. """__slots__=("_all_touched",)# slots *suffixed* with an underscore will be ignored by `__eq__` and `__repr__`.# a double-underscore prefix would be more conventional, but that actually breaks as a slot name.def__init__(self,all_touched=False):""" Parameters ---------- all_touched: bool, default False If True, this ensures that any source pixel which intersects the AOI GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). """self._all_touched=bool(all_touched)def__getstate__(self):return{attr:getattr(self,attr)forsinself.__class__.__mro__forattringetattr(s,"__slots__",tuple())}def__setstate__(self,state):forattr,valinstate.items():setattr(self,attr,val)@propertydefall_touched(self):""" bool: If True, this ensures that any source pixel which intersects the GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). """returnself._all_touched@propertydefraster_params(self):""" dict: The properties of this GeoContext, as keyword arguments to use for `Raster.ndarray` or `Raster.raster`. """raster_params={}ifself.all_touched:raster_params["cutline_all_touched"]=Truereturnraster_paramsdef__eq__(self,other):""" Two GeoContexts are equal only if they are the same type, and every property is equal. """ifnotisinstance(other,self.__class__):returnFalseforattrinself.__slots__:ifgetattr(self,attr)!=getattr(other,attr):returnFalsereturnTruedef__repr__(self):classname=self.__class__.__name__delim=",\n"+" "*(len(classname)+1)props=delim.join("{}={}".format(attr.lstrip("_"),reprlib.repr(getattr(self,attr)))forsinself.__class__.__mro__forattringetattr(s,"__slots__",tuple()))return"{}({})".format(classname,props)
[docs]classAOI(GeoContext):""" A GeoContext that clips imagery to a geometry, and/or to square bounds, with any output resolution and CRS. Examples -------- .. code-block:: python cutline_aoi = dl.geo.AOI(my_geometry, resolution=40) aoi_with_cutline_disabled = cutline_aoi.assign(geometry=None) no_cutline_aoi = dl.geo.AOI(geometry=None, resolution=15, bounds=(-40, 35, -39, 36)) aoi_without_auto_bounds = dl.geo.AOI(geometry=my_geometry, resolution=15, bounds=(-40, 35, -39, 36)) aoi_with_specific_pixel_dimensions = dl.geo.AOI(geometry=my_geometry, shape=(200, 400)) """__slots__=("_geometry","_resolution","_crs","_align_pixels","_bounds","_bounds_crs","_shape",)def__init__(self,geometry=None,resolution=None,crs=None,align_pixels=None,bounds=None,bounds_crs="EPSG:4326",shape=None,all_touched=False,):""" Parameters ---------- geometry: GeoJSON-like dict, object with ``__geo_interface__``; optional When searching, filter for elements which intersect this geometry. When rastering, clip imagery to this geometry. Coordinates must be WGS84 (lat-lon). If :const:`None`, imagery will just be clipped to :py:attr:`~descarteslabs.common.gecontext.AOI.bounds`. resolution: float, optional Distance, in native units of the CRS, that the edge of each pixel represents on the ground. Do not assume this to always be either degrees or meters. Can only specify one of `resolution` and `shape`. crs: str, optional Coordinate Reference System into which imagery will be projected, expressed as an EPSG code (like :const:`EPSG:4326`), a PROJ.4 definition, or an OGC CRS Well-Known Text string. align_pixels: bool, optional, default True if resolution is not None If :const:`True`, this ensures that, in different images rasterized with this same AOI GeoContext, pixels ``(i, j)`` correspond to the same area in space. This is accomplished by snapping the coordinates of the origin (top-left corner of top-left pixel) to a non-fractional interval of `resolution`. Note that in cases where `shape` has been specified, this may lead to the resulting image being one pixel larger in each dimension, so the the entire bounds is included. If `align_pixels` is :const:`False`, when using imagery with different native resolutions and/or projections, pixels at the same indices can be misaligned by a fraction of `resolution` (i.e. correspond to *slighly* different coordinates in space). However, this requires warping of the original image, which can be undesireable when you want to work with the original data in its native resolution and projection. bounds: 4-tuple, optional Clip imagery to these ``(min_x, min_y, max_x, max_y)`` bounds, expressed in :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds_crs` (which defaults to WGS84 lat-lon). :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds` are automatically computed from `geometry` if not specified. Otherwise, :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds` are required. bounds_crs: str, optional, default "EPSG:4326" The Coordinate Reference System of the :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds`, given as an EPSG code (like :const:`EPSG:4326`), a PROJ.4 definition, or an OGC CRS Well-Known Text string. shape: 2-tuple, optional ``(rows, columns)``, in pixels, the output raster should fit within; the longer side of the raster will be min(shape). Can only specify one of `resolution` and `shape`. Note that when `align_pixels` is :const:`True`, the actual resulting raster may be one pixel larger in each direction. all_touched: bool, default False If True, this ensures that any source pixel which intersects the AOI GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). """super(AOI,self).__init__(all_touched=all_touched)# If no bounds were given, use the bounds of the geometryifboundsisNoneandgeometryisnotNone:bounds="update"self._assign(geometry,resolution,crs,align_pixels,bounds,bounds_crs,shape,"unchanged",)self._validate()@propertydefgeometry(self):""" shapely geometry: Clip imagery to this geometry Coordinates must be WGS84 (lat-lon). If :const:`None`, imagery will just be clipped to :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds`. """returnself._geometry@propertydefresolution(self):""" float: Distance, in units of the CRS, that the edge of each pixel represents on the ground. """returnself._resolution@propertydefcrs(self):""" str: Coordinate reference system into which imagery will be projected, expressed as an EPSG code (like :const:`EPSG:4326`), a PROJ.4 definition, or an OGC CRS Well-Known Text string. """returnself._crs@propertydefalign_pixels(self):""" bool: If True, this ensures that, in different images rasterized with this same AOI GeoContext, pixels ``(i, j)`` correspond to the same area in space. This is accomplished by snapping the coordinates of the origin (top-left corner of top-left pixel) to a non-fractional interval of `resolution`. Note that in cases where `shape` has been specified, this may lead to the resulting image being one pixel larger in each dimension, so the the entire bounds is included. If `align_pixels` is False, when using imagery with different native resolutions and/or projections, pixels at the same indicies can be misaligned by a fraction of ``resolution`` (i.e. correspond to *slighly* different coordinates in space). However, this requires warping of the original image, which can be undesireable when you want to work with the original data in its native resolution and projection. """ifself._align_pixelsisNone:returnself._resolutionisnotNoneelse:returnself._align_pixels@propertydefbounds(self):""" tuple: Clip imagery to these ``(min_x, min_y, max_x, max_y)`` bounds, expressed in the coordinate reference system in :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds_crs`. """returnself._bounds@propertydefbounds_crs(self):""" str: The coordinate reference system of the :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds`, given as an EPSG code (like :const:`EPSG:4326`), a PROJ.4 definition, or an OGC CRS Well-Known Text string. """returnself._bounds_crs@propertydefshape(self):""" tuple: ``(rows, columns)``, in pixels, the output raster should fit within; the longer side of the raster will be min(shape). """returnself._shape@propertydefraster_params(self):""" dict: The properties of this `AOI`, as keyword arguments to use for :class:`~descarteslabs.client.services.raster.raster.Raster.ndarray` or :class:`~descarteslabs.client.services.raster.raster.Raster.raster`. Raises ValueError if :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds`, `crs`, :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds_crs`, `resolution`, or `align_pixels` is :const:`None`. """# Ensure that there can be no ambiguity: every parameter must be specified,# so every raster call using this context will return spatially equivalent dataifself._boundsisNone:raiseValueError("AOI must have bounds specified")ifself._bounds_crsisNone:raiseValueError("AOI must have bounds_crs specified")ifself._crsisNone:raiseValueError("AOI must have CRS specified")ifself._resolutionisNoneandself._shapeisNone:raiseValueError("AOI must have one of resolution or shape specified")# align_pixels will always be True or False based on resolution# all_touched doesn't affect the spatial equivalencecutline=(self._geometry.__geo_interface__ifself._geometryisnotNoneelseNone)dimensions=((self._shape[1],self._shape[0])ifself._shapeisnotNoneelseNone)return{**super().raster_params,"cutline":cutline,"resolution":self._resolution,"srs":self._crs,"bounds_srs":self._bounds_crs,"align_pixels":self.align_pixels,"bounds":self._bounds,"dimensions":dimensions,}@propertydef__geo_interface__(self):""" dict: :py:attr:`~descarteslabs.common.geo.geocontext.AOI.geometry` as a GeoJSON Geometry dict, otherwise :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds` as a GeoJSON Polygon dict if :py:attr:`~descarteslabs.common.geo.geocontext.AOI.geometry` is :const:`None` and :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds_crs` is :const:`EPSG:4326`, otherwise raises :exc:`RuntimeError`. """ifself._geometryisnotNone:returnself._geometry.__geo_interface__elifself._boundsisnotNoneandis_wgs84_crs(self._bounds_crs):returnpolygon_from_bounds(self._bounds)else:raiseRuntimeError("AOI GeoContext must have a geometry set, or bounds set and a WGS84 `bounds_crs`, ""to have a __geo_interface__")
[docs]defassign(self,geometry="unchanged",resolution="unchanged",crs="unchanged",align_pixels="unchanged",bounds="unchanged",bounds_crs="unchanged",shape="unchanged",all_touched="unchanged",):""" Return a copy of the AOI with the given values assigned. Note ---- If you are assigning a new geometry and want bounds to updated as well, use ``bounds="update"``. This will also change :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds_crs` to :const:`EPSG:4326`, since the geometry's coordinates are in WGS84 decimal degrees, so the new bounds determined from those coordinates must be in that CRS as well. If you assign :py:attr:`~descarteslabs.common.geo.geocontext.AOI.geometry` without changing :py:attr:`~descarteslabs.common.geo.geocontext.AOI.bounds`, the new AOI GeoContext will produce rasters with the same shape and covering the same spatial area as the old one, just with pixels masked out that fall outside your new geometry. Returns ------- new : `AOI` """new=copy.deepcopy(self)new._assign(geometry,resolution,crs,align_pixels,bounds,bounds_crs,shape,all_touched,)new._validate()returnnew
def_validate(self):# validate shapeifself._shapeisnotNone:ifnotisinstance(self._shape,(list,tuple))orlen(self._shape)!=2:raiseTypeError("Shape must be a tuple of (rows, columns) in pixels")# validate resolutionifself._resolutionisnotNone:ifnotisinstance(self._resolution,(int,float)):raiseTypeError("Resolution must be an int or float, got type '{}'".format(type(self._resolution).__name__))ifself._resolution<=0:raiseValueError("Resolution must be greater than zero")# can't set both resolution and shapeifself._resolutionisnotNoneandself._shapeisnotNone:raiseValueError("Cannot set both resolution and shape")# test that bounds are saneifself._boundsisnotNone:shapely_support.check_valid_bounds(self._bounds)# rough check that bounds values actually make sense for bounds_crsifself._bounds_crsisnotNoneandself._boundsisnotNone:is_geographic,lon_wrap=is_geographic_crs(self._bounds_crs,with_lon_wrap=True)ifis_geographic:# some whole-globe products are funky around the dateline. Try# to allow up to a 1/2 pixel slop there. This will generally only# occur with AOIs created automatically from Image properties.ifself._resolutionandself._crsandis_geographic_crs(self._crs):tol=self._resolution/2elifself._shapeisnotNone:tol=(max((self._bounds[2]-self._bounds[0])/self._shape[1],(self._bounds[3]-self._bounds[1])/self._shape[0],)/2)else:tol=0.001ifnotvalid_latlon_bounds(self._bounds,tol,lon_wrap=lon_wrap):raiseValueError("Bounds must be in lat-lon coordinates, ""but the given bounds are outside [-90, 90] for y or [-180, 180] for x.")else:ifvalid_latlon_bounds(self._bounds):# Warn that bounds are probably in the wrong CRS.# But we can't be sure without a proper tool for working with CRSs,# since bounds that look like valid lat-lon coords# *could* be valid in a different CRS, though unlikely.warnings.warn("You might have the wrong `bounds_crs` set.\n""Bounds appear to be in lat-lon decimal degrees, but the `bounds_crs` ""does not seem to be a geographic coordinate reference system ""(i.e. its units are not degrees, but meters, feet, etc.).\n\n""If this is unexpected, set `bounds_crs='EPSG:4326'`.")# check that bounds and geometry actually intersect (if bounds in wgs84)if(self._geometryisnotNoneandself._boundsisnotNoneandis_wgs84_crs(self._bounds_crs)):bounds_shp=shapely.geometry.box(*self._bounds)ifnotbounds_shp.intersects(self._geometry):raiseValueError("Geometry and bounds do not intersect. This would result in all data being masked. ""If you're assigning new geometry, assign new bounds as well ""(use `bounds='update'` to use the bounds of the new geometry).")# Helpful warning about a common mistake: resolution < width# The CRS of bounds and CRS of resolution must be the same to compare between those values# This most often happens when switching from a projected to a geodetic CRS (i.e. UTM to WGS84)# and not updating the (units of the) resolution accordingly, so you now have, say,# 30 decimal degrees as your resolution. Probably not what you meant.# TODO: better way of checking equivalence between CRSs than string equalityif(notself._all_touchedandself._crsisnotNoneandself._resolutionisnotNoneandself._boundsisnotNoneandself._bounds_crs==self._crs):crs_width=self._bounds[2]-self._bounds[0]crs_height=self._bounds[3]-self._bounds[1]msg=("Output raster's {dim} ({dim_len:.4f}) is smaller than its resolution ""({res:.4f}), meaning it would be less than one pixel {dim_adj}.\n""Remember that resolution is specified in units of the output CRS, ""which are not necessarily meters.")ifis_geographic_crs(self._crs):msg+="\nSince your CRS is in lat-lon coordinates, resolution must be given in decimal degrees."msg+=("\nIf you are intending to raster an area smaller than the source imagery resolution, then you""should set an appropriate value of resolution, shape, or all_touched=True on the supplied AOI"" to signal your intentions.")ifcrs_width<self._resolution:raiseValueError(msg.format(dim="width",dim_len=crs_width,res=self._resolution,dim_adj="wide",))ifcrs_height<self._resolution:raiseValueError(msg.format(dim="height",dim_len=crs_height,res=self._resolution,dim_adj="tall",))def_assign(self,geometry,resolution,crs,align_pixels,bounds,bounds_crs,shape,all_touched,):# we use "unchanged" as a sentinel value, because None is a valid thing to set attributes to.ifgeometryisnotNoneandgeometry!="unchanged":geometry=shapely_support.geometry_like_to_shapely(geometry)ifboundsisnotNoneandbounds!="unchanged":ifbounds=="update":ifbounds_crsnotin(None,"unchanged","EPSG:4326"):raiseValueError("Can't compute bounds from a geometry while also explicitly setting a `bounds_crs`.\n\n""To resolve: don't set `bounds_crs`. It will be set to 'EPSG:4326' for you. ""(Though you can do so explicitly if you'd like.)\n\n""Explanation: the coordinates in a geometry are latitudes and longitudes ""in decimal degrees, defined in the WGS84 coordinate reference system ""(referred to by the code EPSG:4326). When we infer `bounds` from a `geometry`, ""those bounds will be in the same coordinate reference system as the geometry---i.e., WGS84. ""Therefore, setting `bounds_crs` to anything besides 'EPSG:4326' doesn't make sense.")bounds_crs="EPSG:4326"ifgeometryisnotNoneandgeometry!="unchanged":bounds=geometry.boundselse:raiseValueError("A geometry must be given with which to update the bounds")else:bounds=tuple(bounds)ifgeometry!="unchanged":self._geometry=geometryifresolution!="unchanged":# To avoid breaking existing code, avoid a conflict with shape.# getattr() to handle pre-init cases.if(getattr(self,"_resolution",None)isNoneandgetattr(self,"_shape",None)isnotNone):self._shape=Noneself._resolution=resolutionifcrs!="unchanged":self._crs=crsifalign_pixels!="unchanged":self._align_pixels=align_pixelsifbounds!="unchanged":self._bounds=boundsifbounds_crs!="unchanged":self._bounds_crs=bounds_crsifshape!="unchanged":self._shape=shapeifall_touched!="unchanged":self._all_touched=bool(all_touched)
[docs]classDLTile(GeoContext):""" A GeoContext that clips and projects imagery to a single DLTile. DLTiles allow you to define a grid of arbitrary spacing, resolution, and overlap that can cover the globe. DLTiles are always in a UTM projection. Example ------- >>> import descarteslabs as dl >>> from descarteslabs.geo import DLTile >>> tile = DLTile.from_latlon( ... lat=35.691, ... lon=-105.944, ... tilesize=512, ... resolution=10, ... pad=0 ... ) >>> product = dl.catalog.Product.get("usgs:landsat:oli-tirs:c2:l2:v0") # doctest: +SKIP >>> images = product.images().intersects(tile).collect() # doctest: +SKIP >>> images # doctest: +SKIP ImageCollection of 558 images * Dates: Mar 18, 2013 to Sep 14, 2023 * Products: usgs:landsat:oli-tirs:c2:l2:v0: 558 >>> images.geocontext # doctest: +SKIP DLTile(key='512:0:10.0:13:-17:771', resolution=10.0, tilesize=512, pad=0, crs='EPSG:32613', bounds=(412960.0, 3947520.0, 418080.0, 3952640.0), bounds_crs='EPSG:32613', geometry=<POLYGON ((-1....962 35.71...>, zone=13, ti=-17, tj=771, geotrans=(412960.0, 10.0, 0.0, 3952640.0, 0.0, -10.0), proj4='+proj=utm +z...s=m +no_defs ', wkt='PROJCS["WGS ...SG","32613"]]', all_touched=False) """__slots__=("_key","_resolution","_tilesize","_pad","_crs","_bounds","_bounds_crs","_geometry","_zone","_ti","_tj","_geotrans","_proj4","_wkt",)def__init__(self,dltile_dict,all_touched=False):""" Constructs a DLTile from a parameter dictionary. It is preferred to use the :meth:`DLTile.from_latlon, :meth:`DLTile.from_shape`, or :meth:`DLTile.from_key` class methods to construct a DLTile GeoContext. Parameters ---------- dltile_dict: Dict[Str, Any] Dictionary for the tile. all_touched: bool, default False If True, this ensures that any source pixel which intersects the AOI GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). """super(DLTile,self).__init__(all_touched=all_touched)ifisinstance(dltile_dict["geometry"],shapely.geometry.polygon.Polygon):self._geometry=dltile_dict["geometry"]else:self._geometry=shapely.geometry.shape(dltile_dict["geometry"])properties=dltile_dict["properties"]self._key=properties["key"]self._resolution=properties["resolution"]self._tilesize=properties["tilesize"]self._pad=properties["pad"]self._crs=properties["cs_code"]self._bounds=tuple(properties["outputBounds"])self._bounds_crs=properties["cs_code"]self._zone=properties["zone"]self._ti=properties["ti"]self._tj=properties["tj"]# these properties may not be presentself._geotrans=properties.get("geotrans",None)self._proj4=properties.get("proj4",None)self._wkt=properties.get("wkt",None)
[docs]@classmethoddeffrom_latlon(cls,lat,lon,resolution,tilesize,pad,all_touched=False):""" Return a DLTile GeoContext that covers a latitude/longitude. Where the point falls within the tile will vary, depending on the point and tiling parameters. Parameters ---------- lat : float Latitude (WGS84) lon : float Longitude (WGS84) resolution : float Distance, in meters, that the edge of each pixel represents on the ground tilesize : int Length of each side of the tile, in pixels pad : int Number of extra pixels by which each side of the tile is buffered. This determines the number of pixels by which two tiles overlap. all_touched: bool, default False If True, this ensures that any source pixel which intersects the AOI GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). Returns ------- tile : DLTile Example ------- >>> from descarteslabs.geo import DLTile >>> # make a tile with total size 100, centered on lat, lon >>> # total tilesize == tilesize + 2 * pad >>> params = { ... "lat": 30.0131, ... "lon": 31.2089, ... "resolution": 10, ... "tilesize": 2, ... "pad": 49, ... } >>> tile = DLTile.from_latlon(**params) >>> tile.key '2:49:10.0:36:-8637:166079' >>> tile.geometry.centroid.xy # doctest: +SKIP (array('d', [31.20899205942612]), array('d', [30.013121672688087])) """grid=Grid(resolution=resolution,tilesize=tilesize,pad=pad)tile=grid.tile_from_lonlat(lat=lat,lon=lon)returncls(tile.geocontext,all_touched=all_touched)
[docs]@classmethoddeffrom_shape(cls,shape,resolution,tilesize,pad,keys_only=False,all_touched=False):""" Return a list of DLTiles that intersect the given geometry. Parameters ---------- shape : GeoJSON-like A GeoJSON dict, or object with a ``__geo_interface__``. Must be in :const:`EPSG:4326` (WGS84 lat-lon) projection. resolution : float Distance, in meters, that the edge of each pixel represents on the ground. tilesize : int Length of each side of the tile, in pixels. pad : int Number of extra pixels by which each side of the tile is buffered. This determines the number of pixels by which two tiles overlap. keys_only : bool, default False Whether to return DLTile objects or only DLTile keys. Set to True when returning a large number of tiles and you do not need the full objects. all_touched: bool, default False If True, this ensures that any source pixel which intersects the AOI GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). Returns ------- tiles : List[DLTile] or List[Str] Example ------- >>> from descarteslabs.geo import DLTile >>> shape = { ... "type":"Feature", ... "geometry":{ ... "type":"Polygon", ... "coordinates":[[ ... [-122.51140471760839,37.77130087547876], ... [-122.45475646845254,37.77475476721895], ... [-122.45303985468301,37.76657207194229], ... [-122.51057242081689,37.763446782666094], ... [-122.51140471760839,37.77130087547876]]] ... },"properties": None ... } >>> tiles = DLTile.from_shape( ... shape=shape, ... resolution=1, ... tilesize=500, ... pad=0, ... ) >>> len(tiles) 31 """grid=Grid(resolution=resolution,tilesize=tilesize,pad=pad)ifgrid._estimate_ntiles_from_shape(shape)>50000:warnings.warn("DLTile.from_shape will return a large number of tiles. ""Consider using DLTile.iter_from_shape instead.")tiles=grid.tiles_from_shape(shape=shape,keys_only=keys_only)ifkeys_only:result=[tilefortileintiles]else:result=[cls(tile.geocontext,all_touched=all_touched)fortileintiles]returnresult
[docs]@classmethoddefiter_from_shape(cls,shape,resolution,tilesize,pad,keys_only=False,all_touched=False):""" Return a iterator for DLTiles that intersect the given geometry. Parameters ---------- shape : GeoJSON-like A GeoJSON dict, or object with a ``__geo_interface__``. Must be in :const:`EPSG:4326` (WGS84 lat-lon) projection. resolution : float Distance, in meters, that the edge of each pixel represents on the ground. tilesize : int Length of each side of the tile, in pixels. pad : int Number of extra pixels by which each side of the tile is buffered. This determines the number of pixels by which two tiles overlap. keys_only : bool, default False Whether to return DLTile objects or only DLTile keys. Set to True when returning a large number of tiles and you do not need the full objects. all_touched: bool, default False If True, this ensures that any source pixel which intersects the AOI GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). Returns ------- Iterator of DLTiles or str Example ------- >>> from descarteslabs.geo import DLTile >>> shape = { ... "type":"Feature", ... "geometry":{ ... "type":"Polygon", ... "coordinates":[[ ... [-122.51140471760839,37.77130087547876], ... [-122.45475646845254,37.77475476721895], ... [-122.45303985468301,37.76657207194229], ... [-122.51057242081689,37.763446782666094], ... [-122.51140471760839,37.77130087547876]]] ... },"properties": None ... } >>> gen = DLTile.from_shape( ... shape=shape, ... resolution=1, ... tilesize=500, ... pad=0, ... keys_only=True ... ) >>> tiles = [tile for tile in gen] # doctest: +SKIP >>> tiles[0] # doctest: +SKIP '500:0:1.0:10:94:8359' """grid=Grid(resolution=resolution,tilesize=tilesize,pad=pad)tiles=grid.tiles_from_shape(shape=shape,keys_only=keys_only)fortileintiles:ifkeys_only:yieldtileelse:yieldcls(tile.geocontext,all_touched=all_touched)
[docs]@classmethoddeffrom_key(cls,dltile_key,all_touched=False):""" Return a DLTile GeoContext from a DLTile key. Parameters ---------- dltile_key : str DLTile key, e.g. '128:16:960.0:15:-1:37' all_touched: bool, default False If True, this ensures that any source pixel which intersects the AOI GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). Returns ------- tile: DLTile Example ------- >>> from descarteslabs.geo import DLTile >>> tile = DLTile.from_key("2048:16:30.0:15:3:80") >>> tile # doctest: +SKIP DLTile(key='2048:16:30.0:15:3:80', resolution=30.0, tilesize=2048, pad=16, crs='EPSG:32615', bounds=(683840.0, 4914720.0, 746240.0, 4977120.0), bounds_crs='EPSG:32615', geometry=<shapely.geom...>, zone=15, ti=3, tj=80, geotrans=[ ... """tile=Tile.from_key(dltile_key)returncls(tile.geocontext,all_touched=all_touched)
[docs]defsubtile(self,subdivide,resolution=None,pad=None,keys_only=False):""" Return an iterator for new DLTiles that subdivide this tile. The DLtile will be sub-divided into subdivide^2 total sub-tiles each with a side length of tile_size / subdivide. The resulting sub-tile size must be an integer. Each sub-tile will by default inherit the same resolution and pad as the orginal tile. Parameters ---------- subdivide : int The value to subdivide the tile. The total number of sub-tiles will be the square of this value. This value must evenly divide the original tilesize. resolution : None, float A new resolution for the sub-tiles. None defaults to the original DLTile resolution. The new resolution must evenly divide the the original tilesize divided by the subdivide ratio. pad : None, int A new pad value for the sub-tiles. None defaults to the original DLTile pad value. keys_only : bool, default False Whether to return DLTile objects or only DLTile keys. Set to True when returning a large number of tiles and you do not need the full objects. Returns ------- Iterator over DLTiles or str Example: ------- >>> from descarteslabs.geo import DLTile >>> tile = DLTile.from_key("2048:0:30.0:15:3:80") >>> tiles = [tile for tile in tile.subtile(8)] >>> len(tiles) 64 >>> tiles[0].tilesize 256 """subtiles=Tile.from_key(self.key).subtile(subdivide=subdivide,new_resolution=resolution,new_pad=pad,)fortileinsubtiles:ifkeys_only:yieldtile.keyelse:yieldDLTile(tile.geocontext,all_touched=self.all_touched)
[docs]defrowcol_to_latlon(self,row,col):""" Convert pixel coordinates to lat, lon coordinates Parameters ---------- row : int or List[int] Pixel row coordinate or coordinates col : int or List[int] Pixel column coordinate or coordinates Returns ------- coords : List[Tuple[float], Tuple[float]] List with the first element the latitude values and the second element longitude values Example ------- >>> from descarteslabs.geo import DLTile >>> tile = DLTile.from_key("2048:0:30.0:15:3:80") >>> tile.rowcol_to_latlon(row=56, col=1111) [(44.894653081367544,), (-90.24334206726267,)] """lonlat=Tile.from_key(self.key).rowcol_to_lonlat(row=row,col=col)lonlat=lonlat.tolist()ifisinstance(lonlat[0],(int,float)):result=[(lonlat[1],),(lonlat[0],)]else:result=list(zip(*lonlat))result[0],result[1]=result[1],result[0]returnresult
[docs]deflatlon_to_rowcol(self,lat,lon):""" Convert lat, lon coordinates to pixel coordinates Parameters ---------- lat: float or List[float] Latitude coordinate or coordinates lon: float or List[float] Longitude coordinate or coordinates Returns ------- coords: List[Tuple[int] Tuple[int]] Tuple with the first element the row values and the second element column values Example ------- >>> from descarteslabs.geo import DLTile >>> tile = DLTile.from_key("2048:0:30.0:15:3:80") >>> tile.latlon_to_rowcol(lat=44.8, lon=-90.2) [(403,), (1237,)] """rowcol=Tile.from_key(self.key).lonlat_to_rowcol(lat=lat,lon=lon)rowcol=rowcol.tolist()ifisinstance(rowcol[0],(int,float)):result=[(rowcol[0],),(rowcol[1],)]else:result=list(zip(*rowcol))returnresult
[docs]defassign(self,pad="unchanged",all_touched="unchanged"):""" Return a copy of the DLTile with the pad and/or all_touched value modified. Parameters ---------- pad : int, default "unchanged" New pad value all_touched : bool, default "unchanged" New all_touched value Returns ------- tile : DLTile Example: -------- >>> from descarteslabs.geo import DLTile >>> tile = DLTile.from_key("2048:16:30.0:15:3:80") >>> tile.pad 16 >>> tile = tile.assign(123) >>> tile.pad 123 """tile=Tile.from_key(self.key)ifpad!="unchanged":tile=tile.assign(pad=pad)ifall_touched=="unchanged":all_touched=self.all_touchedreturnDLTile(tile.geocontext,all_touched=all_touched)
@propertydefkey(self):""" str: The DLTile's key, which encodes the tiling parameters, and which number in the grid this tile is. """returnself._key@propertydefresolution(self):"""float: Distance, in meters, that the edge of each pixel represents on the ground"""returnself._resolution@propertydeftilesize(self):""" int: Length of each side of the tile, in pixels. Note that the total number of pixels along each side of an image is ``tile_size + 2 * padding`` """returnself._tilesize@propertydeftile_extent(self):""" int: total extent of geocontext length in pixels, including pad. Size is ``tile_size + 2 * pad``. """returnself._tilesize+2*self._pad@propertydefpad(self):""" int: Number of extra pixels by which each side of the tile is buffered. This determines the number of pixels by which two tiles overlap. """returnself._pad@propertydefcrs(self):""" str: Coordinate reference system into which imagery will be projected. For DLTiles, this is always a UTM projection, given as an EPSG code. """returnself._crs@propertydefbounds(self):""" tuple: The ``(min_x, min_y, max_x, max_y)`` of the area covered by this DLTile, in the UTM coordinate reference system given in :py:attr:`~descarteslabs.common.geo.geocontext.DLTile.bounds_crs`. """returnself._bounds@propertydefbounds_crs(self):""" str: The coordinate reference system of the :py:attr:`~descarteslabs.common.geo.geocontext.DLTile.bounds`, given as an EPSG code (like :const:`EPSG:32615`). A DLTile's CRS is always UTM. """returnself._bounds_crs@propertydefgeometry(self):""" shapely.geometry.Polygon: The polygon covered by this DLTile in WGS84 (lat-lon) coordinates """returnself._geometry@propertydefzone(self):"""int: The UTM zone of this tile"""returnself._zone@propertydefti(self):"""int: The y-index of this tile in its grid"""returnself._ti@propertydeftj(self):"""int: The x-index of this tile in its grid"""returnself._tj@propertydefraster_params(self):""" dict: The properties of this DLTile, as keyword arguments to use for `Raster.ndarray` or `Raster.raster`. """return{**super().raster_params,"dltile":self._key,# QUESTION: shouldn't align_pixels be True?# based on the GDAL documentation for `-tap`, seems like that should be true# to ensure that pixels of images with different resolutions/projections# are aligned with the same dltile. otherwise, pixel (0,0) in 1 image could be at# different coordinates than the other"align_pixels":False,}@propertydefgeotrans(self):""" tuple: The 6-tuple GDAL geotrans for this DLTile in the shape ``(a, b, c, d, e, f)`` where | a is the top left pixel's x-coordinate | b is the west-east pixel resolution | c is the row rotation, always 0 for DLTiles | d is the top left pixel's y-coordinate | e is the column rotation, always 0 for DLTiles | f is the north-south pixel resolution, always a negative value """ifself._geotransisNone:returnNonereturntuple(self._geotrans)@propertydefproj4(self):"""str: PROJ.4 definition for this DLTile's coordinate reference system"""returnself._proj4@propertydefwkt(self):"""str: OGC Well-Known Text definition for this DLTile's coordinate reference system"""returnself._wkt@propertydef__geo_interface__(self):"""dict: :py:attr:`~descarteslabs.common.geo.geocontext.DLTile.geometry` as a GeoJSON Polygon"""returnself._geometry.__geo_interface__
[docs]classXYZTile(GeoContext):""" A GeoContext for XYZ tiles, such as those used in web maps. The tiles are always 256x256 pixels, in the spherical Mercator or "Web Mercator" coordinate reference system (:const:`EPSG:3857`). """__slots__=("_x","_y","_z")def__init__(self,x,y,z,all_touched=False):""" Parameters ---------- x: int X-index of the tile (increases going east) y: int Y-index of the tile (increases going south) z: int Zoom level of the tile all_touched: bool, default False If True, this ensures that any source pixel which intersects the AOI GeoContext contributes to the raster result. Normally this mode is not enabled, and its use is strongly discouraged. However, it can be useful when the AOI is smaller than a source pixel, which under many situations will return no result at all (i.e. entirely masked). """self._x=xself._y=yself._z=zsuper(XYZTile,self).__init__(all_touched=all_touched)@propertydefx(self):"int: X-index of the tile (increases going east)"returnself._x@propertydefy(self):"int: Y-index of the tile (increases going south)"returnself._y@propertydefz(self):"int: Zoom level of the tile"returnself._z
[docs]defparent(self):"The parent XYZTile enclosing this one"returnself.__class__(*mercantile.parent(self._x,self._y,self._z))
[docs]defchildren(self):"List of child XYZTiles contained within this one"return[self.__class__(*t)fortinmercantile.children(self._x,self._y,self._z)]
@propertydefgeometry(self):""" shapely.geometry.Polygon: The polygon covered by this XYZTile in :const:`WGS84` (lat-lon) coordinates """returnshapely.geometry.box(*mercantile.bounds(self._x,self._y,self._z))@propertydefbounds(self):""" tuple: The ``(min_x, min_y, max_x, max_y)`` of the area covered by this XYZTile, in spherical Mercator coordinates (EPSG:3857). """returntuple(mercantile.xy_bounds(self._x,self._y,self._z))@propertydefcrs(self):""" str: Coordinate reference system into which common.geo will be projected. Always :const:`EPSG:3857` (spherical Mercator, aka "Web Mercator") """return"EPSG:3857"@propertydefbounds_crs(self):""" str: The coordinate reference system of the :py:attr:`~descarteslabs.common.geo.geocontext.XYZTile.bounds`. Always :const:`EPSG:3857` (spherical Mercator, aka "Web Mercator") """return"EPSG:3857"@propertydeftilesize(self):""" int: Length of each side of the tile, in pixels. Always 256. """return256@propertydefresolution(self):""" float: Distance, in meters, that the edge of each pixel represents in the spherical Mercator ("Web Mercator", EPSG:3857) projection. """num_tiles=1<<self.zreturnEARTH_CIRCUMFERENCE_WGS84/num_tiles/self.tilesize@propertydef__geo_interface__(self):"dict: :py:attr:`~descarteslabs.common.geo.geocontext.XYZTile.geometry` as a GeoJSON Polygon"returnself.geometry.__geo_interface__@propertydefraster_params(self):""" dict: The properties of this XYZTile, as keyword arguments to use for `Raster.ndarray` or `Raster.raster`. """return{**super().raster_params,"bounds":self.bounds,"srs":self.crs,"bounds_srs":self.bounds_crs,"align_pixels":False,"resolution":self.resolution,}