зеркало из https://github.com/microsoft/landcover.git
Removed the concept of "shapes" on the server side
Fixed the DataLoaderUSALayer and added an appropriate dataset entry in datasets.json
This commit is contained in:
Caleb Robinson
Родитель
0ef7352bda
Коммит
1bf6e4c576
24
server.py
24
server.py
|
|
@ -7,6 +7,8 @@ import base64
|
|||
import json
|
||||
import logging
|
||||
import os
|
||||
os.environ["CURL_CA_BUNDLE"] = "/etc/ssl/certs/ca-certificates.crt"
|
||||
|
||||
import sys
|
||||
import time
|
||||
|
||||
|
|
@ -19,9 +21,9 @@ import rasterio.warp
|
|||
|
||||
LOGGER = logging.getLogger("server")
|
||||
|
||||
from web_tool.DataLoader import warp_data_to_3857, crop_data_by_extent, crop_data_by_geometry
|
||||
from web_tool.Datasets import load_datasets, get_area_from_geometry
|
||||
DATASETS = load_datasets()
|
||||
from web_tool.DataLoader import warp_data_to_3857, crop_data_by_extent, crop_data_by_geometry, get_area_from_geometry
|
||||
from web_tool.Datasets import load_datasets
|
||||
DATALOADERS = load_datasets()
|
||||
|
||||
from web_tool.Utils import setup_logging, get_random_string, class_prediction_to_img
|
||||
from web_tool import ROOT_DIR
|
||||
|
|
@ -207,10 +209,10 @@ def pred_patch():
|
|||
name_list = [item["name"] for item in class_list]
|
||||
color_list = [item["color"] for item in class_list]
|
||||
|
||||
if dataset not in DATASETS:
|
||||
if dataset not in DATALOADERS:
|
||||
raise ValueError("Dataset doesn't seem to be valid, do the datasets in js/tile_layers.js correspond to those in TileLayers.py")
|
||||
else:
|
||||
current_data_loader = DATASETS[dataset]["data_loader"]
|
||||
current_data_loader = DATALOADERS[dataset]
|
||||
|
||||
input_raster = current_data_loader.get_data_from_extent(extent)
|
||||
current_session.latest_input_raster = input_raster
|
||||
|
|
@ -255,19 +257,19 @@ def pred_tile():
|
|||
zone_layer_name = data["zoneLayerName"]
|
||||
model_idx = data["modelIdx"]
|
||||
|
||||
if dataset not in DATASETS:
|
||||
if dataset not in DATALOADERS:
|
||||
raise ValueError("Dataset doesn't seem to be valid, do the datasets in js/tile_layers.js correspond to those in TileLayers.py")
|
||||
else:
|
||||
current_data_loader = DATASETS[dataset]["data_loader"]
|
||||
current_data_loader = DATALOADERS[dataset]
|
||||
|
||||
try:
|
||||
input_raster = current_data_loader.get_data_from_shape(geom["geometry"])
|
||||
input_raster = current_data_loader.get_data_from_geometry(geom["geometry"])
|
||||
shape_area = get_area_from_geometry(geom["geometry"])
|
||||
except NotImplementedError as e: # Example of how to handle errors from the rest of the server
|
||||
bottle.response.status = 400
|
||||
return json.dumps({"error": "Cannot currently download imagery with 'Basemap' based datasets"})
|
||||
|
||||
output_raster = current_session.pred_tile(input_raster)
|
||||
output_raster = current_session.pred_tile(input_raster)
|
||||
if output_raster.shape[2] > len(color_list):
|
||||
LOGGER.warning("The number of output channels is larger than the given color list, cropping output to number of colors (you probably don't want this to happen")
|
||||
output_raster.data = output_raster.data[:,:,:len(color_list)]
|
||||
|
|
@ -349,10 +351,10 @@ def get_input():
|
|||
extent = data["extent"]
|
||||
dataset = data["dataset"]
|
||||
|
||||
if dataset not in DATASETS:
|
||||
if dataset not in DATALOADERS:
|
||||
raise ValueError("Dataset doesn't seem to be valid, please check Datasets.py")
|
||||
else:
|
||||
current_data_loader = DATASETS[dataset]["data_loader"]
|
||||
current_data_loader = DATALOADERS[dataset]
|
||||
|
||||
input_raster = current_data_loader.get_data_from_extent(extent)
|
||||
warped_output_raster = warp_data_to_3857(input_raster) # warp image to 3857
|
||||
|
|
|
|||
|
|
@ -23,6 +23,7 @@ import rasterio.merge
|
|||
import rtree
|
||||
|
||||
import mercantile
|
||||
import utm
|
||||
|
||||
import cv2
|
||||
import pickle
|
||||
|
|
@ -30,6 +31,8 @@ import pickle
|
|||
from . import ROOT_DIR
|
||||
from .DataLoaderAbstract import DataLoader
|
||||
|
||||
NAIP_BLOB_ROOT = 'https://naipblobs.blob.core.windows.net/naip'
|
||||
|
||||
|
||||
class InMemoryRaster(object):
|
||||
|
||||
|
|
@ -87,7 +90,6 @@ def extent_to_transformed_geom(extent, dst_crs):
|
|||
return fiona.transform.transform_geom(src_crs, dst_crs, geom)
|
||||
|
||||
|
||||
|
||||
def warp_data_to_3857(input_raster):
|
||||
"""Warps an input raster to EPSG:3857
|
||||
|
||||
|
|
@ -185,32 +187,54 @@ def crop_data_by_geometry(input_raster, geometry, geometry_crs):
|
|||
return InMemoryRaster(dst_img, input_raster.crs, dst_transform, (left, bottom, right, top))
|
||||
|
||||
|
||||
def get_area_from_geometry(geom, src_crs="epsg:4326"):
|
||||
"""Semi-accurately calculates the area for an input GeoJSON shape in km^2 by reprojecting it into a local UTM coordinate system.
|
||||
|
||||
Args:
|
||||
geom (dict): A polygon (or multipolygon) in GeoJSON format
|
||||
src_crs (str, optional): The CRS of `geom`. Defaults to "epsg:4326".
|
||||
|
||||
Raises:
|
||||
ValueError: This will be thrown if geom isn't formatted correctly, or is not a Polygon or MultiPolygon type
|
||||
|
||||
Returns:
|
||||
area (float): The area of `geom` in km^2
|
||||
"""
|
||||
|
||||
# get one of the coordinates
|
||||
try:
|
||||
if geom["type"] == "Polygon":
|
||||
lon, lat = geom["coordinates"][0][0]
|
||||
elif geom["type"] == "MultiPolygon":
|
||||
lon, lat = geom["coordinates"][0][0][0]
|
||||
else:
|
||||
raise ValueError("Polygons and MultiPolygons only")
|
||||
except:
|
||||
raise ValueError("Input shape is not in the correct format")
|
||||
|
||||
zone_number = utm.latlon_to_zone_number(lat, lon)
|
||||
hemisphere = "+north" if lat > 0 else "+south"
|
||||
dest_crs = "+proj=utm +zone=%d %s +datum=WGS84 +units=m +no_defs" % (zone_number, hemisphere)
|
||||
projected_geom = fiona.transform.transform_geom(src_crs, dest_crs, geom)
|
||||
area = shapely.geometry.shape(projected_geom).area / 1000000.0 # we calculate the area in square meters then convert to square kilometers
|
||||
return area
|
||||
|
||||
|
||||
# ------------------------------------------------------
|
||||
# DataLoader for arbitrary GeoTIFFs
|
||||
# ------------------------------------------------------
|
||||
class DataLoaderCustom(DataLoader):
|
||||
|
||||
def __init__(self, data_fn, shapes, padding):
|
||||
def __init__(self, padding, **kwargs):
|
||||
"""A `DataLoader` object made for single raster datasources (single .tif files, .vrt files, etc.). This provides functionality for extracting data
|
||||
from different shapes and calculating the area of shapes.
|
||||
|
||||
Args:
|
||||
data_fn (str): Path (relative to the root project directory) to a raster file that GDAL can read
|
||||
shapes (dict): A dictionary of with key:value pairs in the format `<layer name>: {"geoms": [...], "areas": [...], "crs": <str>}`. These describe the different polygon layers from
|
||||
the "datasets.json" file.
|
||||
padding (float): Amount of padding in terms of units of the CRS of the raster source pointed to by `data_fn`
|
||||
padding (float): Amount of padding in terms of units of the CRS of the raster source pointed to by `data_fn`.
|
||||
**kwargs: Should contain a "path" key that points to the location of the datasource (e.g. the .tif or .vrt file to use)
|
||||
"""
|
||||
self.data_fn = data_fn
|
||||
self._shapes = shapes
|
||||
self._padding = padding
|
||||
|
||||
@property
|
||||
def shapes(self):
|
||||
return self._shapes
|
||||
|
||||
@shapes.setter
|
||||
def shapes(self, value):
|
||||
self._shapes = value
|
||||
self.data_fn = kwargs["path"]
|
||||
|
||||
@property
|
||||
def padding(self):
|
||||
|
|
@ -221,15 +245,6 @@ class DataLoaderCustom(DataLoader):
|
|||
self._padding = value
|
||||
|
||||
def get_data_from_extent(self, extent):
|
||||
"""Returns the data from the class' raster source corresponding to a *buffered* version of the input extent.
|
||||
Buffering is done by `self.padding` number of units (in terms of the source coordinate system).
|
||||
|
||||
Args:
|
||||
extent (dict): A geographic extent formatted as a dictionary with the following keys: xmin, xmax, ymin, ymax, crs
|
||||
|
||||
Returns:
|
||||
output_raster (InMemoryRaster): A raster cropped to a *buffered* version of the input extent.
|
||||
"""
|
||||
with rasterio.open(self.data_fn, "r") as f:
|
||||
src_crs = f.crs.to_string()
|
||||
transformed_geom = extent_to_transformed_geom(extent, src_crs)
|
||||
|
|
@ -244,61 +259,16 @@ class DataLoaderCustom(DataLoader):
|
|||
src_image = np.rollaxis(src_image, 0, 3)
|
||||
return InMemoryRaster(src_image, src_crs, src_transform, buffed_geom.bounds)
|
||||
|
||||
def get_area_from_shape_by_extent(self, extent, shape_layer):
|
||||
"""Returns the area from the shape that _contains_ the centroid of the given extent.
|
||||
|
||||
Args:
|
||||
extent (dict): A geographic extent formatted as a dictionary with the following keys: xmin, xmax, ymin, ymax, crs
|
||||
shape_layer (str): The name of a key in `self._shapes` to look through
|
||||
|
||||
Returns:
|
||||
area (float): The area of the shape (or a ValueError is raised if no shape is found)
|
||||
"""
|
||||
i, shape = self.get_shape_by_extent(extent, shape_layer)
|
||||
return self.shapes[shape_layer]["areas"][i]
|
||||
|
||||
def get_shape_by_extent(self, extent, shape_layer):
|
||||
"""Returns the index and shape where shape _contains_ the centroid of the given extent.
|
||||
Specifically, this will iterate through `self._shapes[shape_layer]["geoms"]` to find which shape contains the centroid of the extent.
|
||||
|
||||
Args:
|
||||
extent (dict): A geographic extent formatted as a dictionary with the following keys: xmin, xmax, ymin, ymax, crs
|
||||
shape_layer (str): A key in `self._shapes` to look through
|
||||
|
||||
Raises:
|
||||
ValueError: Will be raised if no shape within `self._shapes[shape_layer]` contains the given extent
|
||||
|
||||
Returns:
|
||||
(index (int), shape (shapely.geometry.shape)): A tuple containing the index of the found shape and the corresponding shapely shape object
|
||||
"""
|
||||
transformed_geom = extent_to_transformed_geom(extent, self.shapes[shape_layer]["crs"])
|
||||
transformed_shape = shapely.geometry.shape(transformed_geom)
|
||||
mask_geom = None
|
||||
for i, shape in enumerate(self.shapes[shape_layer]["geoms"]):
|
||||
if shape.contains(transformed_shape.centroid):
|
||||
return i, shape
|
||||
raise ValueError("No shape contains the centroid")
|
||||
|
||||
def get_data_from_shape(self, shape):
|
||||
"""Returns the data from the class' raster source corresponding to the input `shape` without any buffering applied. Note that `shape` is expected
|
||||
to be a GeoJSON polygon (as a dictionary) in the EPSG:4326 coordinate system.
|
||||
|
||||
Args:
|
||||
shape (dict): A polygon in GeoJSON format describing the boundary to crop the input raster to
|
||||
|
||||
Returns:
|
||||
output_raster (InMemoryRaster): A raster cropped to the outline of `shape`
|
||||
"""
|
||||
|
||||
def get_data_from_geometry(self, geometry):
|
||||
#TODO: Figure out what happens if we call this with a geometry that doesn't intersect the data source.
|
||||
f = rasterio.open(self.data_fn, "r")
|
||||
src_profile = f.profile
|
||||
src_crs = f.crs.to_string()
|
||||
transformed_mask_geom = fiona.transform.transform_geom("epsg:4326", src_crs, shape)
|
||||
transformed_mask_geom = fiona.transform.transform_geom("epsg:4326", src_crs, geometry)
|
||||
src_image, src_transform = rasterio.mask.mask(f, [transformed_mask_geom], crop=True, all_touched=True, pad=False)
|
||||
f.close()
|
||||
|
||||
src_image = np.rollaxis(src_image, 0, 3)
|
||||
#return src_image, src_profile, src_transform, shapely.geometry.shape(transformed_mask_geom).bounds, src_crs
|
||||
return InMemoryRaster(src_image, src_crs, src_transform, shapely.geometry.shape(transformed_mask_geom).bounds)
|
||||
|
||||
|
||||
|
|
@ -309,115 +279,82 @@ class NAIPTileIndex(object):
|
|||
TILES = None
|
||||
|
||||
@staticmethod
|
||||
def lookup(extent):
|
||||
def lookup(geom):
|
||||
if NAIPTileIndex.TILES is None:
|
||||
assert all([os.path.exists(fn) for fn in [
|
||||
ROOT_DIR + "/data/tile_index.dat",
|
||||
ROOT_DIR + "/data/tile_index.idx",
|
||||
ROOT_DIR + "/data/tiles.p"
|
||||
"data/tile_index/tile_index.dat",
|
||||
"data/tile_index/tile_index.idx",
|
||||
"data/tile_index/tiles.p"
|
||||
]]), "You do not have the correct files, did you setup the project correctly"
|
||||
NAIPTileIndex.TILES = pickle.load(open(ROOT_DIR + "/data/tiles.p", "rb"))
|
||||
return NAIPTileIndex.lookup_naip_tile_by_geom(extent)
|
||||
NAIPTileIndex.TILES = pickle.load(open("data/tile_index/tiles.p", "rb"))
|
||||
return NAIPTileIndex.lookup_naip_tile_by_geom(geom)
|
||||
|
||||
@staticmethod
|
||||
def lookup_naip_tile_by_geom(extent):
|
||||
tile_index = rtree.index.Index(ROOT_DIR + "/data/tile_index")
|
||||
def lookup_naip_tile_by_geom(geom):
|
||||
tile_index = rtree.index.Index("data/tile_index/tile_index")
|
||||
|
||||
geom = extent_to_transformed_geom(extent, "EPSG:4269")
|
||||
minx, miny, maxx, maxy = shapely.geometry.shape(geom).bounds
|
||||
geom = shapely.geometry.mapping(shapely.geometry.box(minx, miny, maxx, maxy, ccw=True))
|
||||
|
||||
geom = shapely.geometry.shape(geom)
|
||||
|
||||
intersected_indices = list(tile_index.intersection(geom.bounds))
|
||||
for idx in intersected_indices:
|
||||
intersected_fn = NAIPTileIndex.TILES[idx][0]
|
||||
intersected_geom = NAIPTileIndex.TILES[idx][1]
|
||||
if intersected_geom.contains(geom):
|
||||
print("Found %d intersections, returning at %s" % (len(intersected_indices), intersected_fn))
|
||||
tile_index.close()
|
||||
return intersected_fn
|
||||
|
||||
tile_index.close()
|
||||
if len(intersected_indices) > 0:
|
||||
raise ValueError("Error, there are overlaps with tile index, but no tile completely contains selection")
|
||||
else:
|
||||
raise ValueError("No tile intersections")
|
||||
|
||||
class USALayerGeoDataTypes(Enum):
|
||||
NAIP = 1
|
||||
NLCD = 2
|
||||
LANDSAT_LEAFON = 3
|
||||
LANDSAT_LEAFOFF = 4
|
||||
BUILDINGS = 5
|
||||
LANDCOVER = 6
|
||||
|
||||
class DataLoaderUSALayer(DataLoader):
|
||||
|
||||
@property
|
||||
def shapes(self):
|
||||
return self._shapes
|
||||
@shapes.setter
|
||||
def shapes(self, value):
|
||||
self._shapes = value
|
||||
def __init__(self, padding):
|
||||
self._padding = padding
|
||||
|
||||
@property
|
||||
def padding(self):
|
||||
return self._padding
|
||||
|
||||
@padding.setter
|
||||
def padding(self, value):
|
||||
self._padding = value
|
||||
|
||||
def __init__(self, shapes, padding):
|
||||
self._shapes = shapes
|
||||
self._padding = padding
|
||||
def get_data_from_extent(self, extent):
|
||||
query_geom = extent_to_transformed_geom(extent, "epsg:4326")
|
||||
naip_fn = NAIPTileIndex.lookup(query_geom)
|
||||
|
||||
def get_fn_by_geo_data_type(self, naip_fn, geo_data_type):
|
||||
fn = None
|
||||
with rasterio.open(NAIP_BLOB_ROOT + "/" + naip_fn) as f:
|
||||
src_crs = f.crs.to_string()
|
||||
transformed_geom = extent_to_transformed_geom(extent, src_crs)
|
||||
transformed_geom = shapely.geometry.shape(transformed_geom)
|
||||
|
||||
if geo_data_type == USALayerGeoDataTypes.NAIP:
|
||||
fn = naip_fn
|
||||
elif geo_data_type == USALayerGeoDataTypes.NLCD:
|
||||
fn = naip_fn.replace("/esri-naip/", "/resampled-nlcd/")[:-4] + "_nlcd.tif"
|
||||
elif geo_data_type == USALayerGeoDataTypes.LANDSAT_LEAFON:
|
||||
fn = naip_fn.replace("/esri-naip/data/v1/", "/resampled-landsat8/data/leaf_on/")[:-4] + "_landsat.tif"
|
||||
elif geo_data_type == USALayerGeoDataTypes.LANDSAT_LEAFOFF:
|
||||
fn = naip_fn.replace("/esri-naip/data/v1/", "/resampled-landsat8/data/leaf_off/")[:-4] + "_landsat.tif"
|
||||
elif geo_data_type == USALayerGeoDataTypes.BUILDINGS:
|
||||
fn = naip_fn.replace("/esri-naip/", "/resampled-buildings/")[:-4] + "_building.tif"
|
||||
elif geo_data_type == USALayerGeoDataTypes.LANDCOVER:
|
||||
fn = naip_fn.replace("/esri-naip/", "/resampled-lc/")[:-4] + "_lc.tif"
|
||||
else:
|
||||
raise ValueError("GeoDataType not recognized")
|
||||
buffed_geom = transformed_geom.buffer(self.padding)
|
||||
buffed_geojson = shapely.geometry.mapping(buffed_geom)
|
||||
|
||||
return fn
|
||||
|
||||
def get_shape_by_extent(self, extent, shape_layer):
|
||||
transformed_geom = extent_to_transformed_geom(extent, shapes_crs)
|
||||
transformed_shape = shapely.geometry.shape(transformed_geom)
|
||||
mask_geom = None
|
||||
for i, shape in enumerate(shapes):
|
||||
if shape.contains(transformed_shape.centroid):
|
||||
return i, shape
|
||||
raise ValueError("No shape contains the centroid")
|
||||
src_image, src_transform = rasterio.mask.mask(f, [buffed_geojson], crop=True, all_touched=True, pad=False)
|
||||
|
||||
def get_data_from_extent(self, extent, geo_data_type=USALayerGeoDataTypes.NAIP):
|
||||
naip_fn = NAIPTileIndex.lookup(extent)
|
||||
fn = self.get_fn_by_geo_data_type(naip_fn, geo_data_type)
|
||||
src_image = np.rollaxis(src_image, 0, 3)
|
||||
return InMemoryRaster(src_image, src_crs, src_transform, buffed_geom.bounds)
|
||||
|
||||
f = rasterio.open(fn, "r")
|
||||
src_crs = f.crs
|
||||
transformed_geom = extent_to_transformed_geom(extent, f.crs.to_dict())
|
||||
transformed_geom = shapely.geometry.shape(transformed_geom)
|
||||
buffed_geom = transformed_geom.buffer(self.padding)
|
||||
geom = shapely.geometry.mapping(shapely.geometry.box(*buffed_geom.bounds))
|
||||
src_image, src_transform = rasterio.mask.mask(f, [geom], crop=True)
|
||||
f.close()
|
||||
def get_data_from_geometry(self, geometry):
|
||||
naip_fn = NAIPTileIndex.lookup(geometry)
|
||||
|
||||
return src_image, src_crs, src_transform, buffed_geom.bounds
|
||||
with rasterio.open(NAIP_BLOB_ROOT + "/" + naip_fn) as f:
|
||||
src_profile = f.profile
|
||||
src_crs = f.crs.to_string()
|
||||
transformed_mask_geom = fiona.transform.transform_geom("epsg:4326", src_crs, geometry)
|
||||
src_image, src_transform = rasterio.mask.mask(f, [transformed_mask_geom], crop=True, all_touched=True, pad=False)
|
||||
|
||||
def get_area_from_shape_by_extent(self, extent, shape_layer):
|
||||
raise NotImplementedError()
|
||||
src_image = np.rollaxis(src_image, 0, 3)
|
||||
return InMemoryRaster(src_image, src_crs, src_transform, shapely.geometry.shape(transformed_mask_geom).bounds)
|
||||
|
||||
def get_data_from_shape(self, shape):
|
||||
raise NotImplementedError()
|
||||
|
||||
# ------------------------------------------------------
|
||||
# DataLoader for loading RGB data from arbitrary basemaps
|
||||
|
|
@ -451,7 +388,7 @@ class DataLoaderBasemap(DataLoader):
|
|||
arr = np.asarray(bytearray(req.read()), dtype=np.uint8)
|
||||
img = cv2.imdecode(arr, -1)
|
||||
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
|
||||
return img
|
||||
return img
|
||||
|
||||
def get_tile_as_virtual_raster(self, tile):
|
||||
'''NOTE: Here "tile" refers to a mercantile "Tile" object.'''
|
||||
|
|
|
|||
|
|
@ -1,29 +1,48 @@
|
|||
from abc import ABC, abstractmethod
|
||||
import abc
|
||||
|
||||
class DataLoader(ABC):
|
||||
class DataLoader(abc.ABC):
|
||||
|
||||
@property
|
||||
@abstractmethod
|
||||
@abc.abstractmethod
|
||||
def padding(self):
|
||||
"""A float value that describes the amount of padding (in terms of the class' data source CRS) to apply on input shapes in `get_data_from_extent`.
|
||||
"""
|
||||
pass
|
||||
|
||||
@property
|
||||
@abstractmethod
|
||||
def shapes(self):
|
||||
pass
|
||||
@abc.abstractmethod
|
||||
def __init__(self, padding, **kwargs):
|
||||
"""A `DataLoader` object should be able to query a source of data by polygons / extents and return the result in a reasonable amount of time. This functionality is abstracted
|
||||
as different sources of data will need to be queried using different interfaces. For example, local raster data sources (".tif", ".vrt", etc.) can be simply accessed, while global data
|
||||
sources provided by a basemap will need more effort to access.
|
||||
|
||||
@abstractmethod
|
||||
def get_shape_by_extent(self, extent, shape_layer):
|
||||
Args:
|
||||
padding (float): Amount of padding in terms of units of the CRS of the raster source pointed to by `data_fn` to apply during `get_data_from_extent`.
|
||||
**kwargs: Key, value pairs created from the contents of this implementation's "dataLayer" key in datasets.json.
|
||||
"""
|
||||
raise NotImplementedError()
|
||||
|
||||
@abstractmethod
|
||||
|
||||
@abc.abstractmethod
|
||||
def get_data_from_extent(self, extent):
|
||||
"""Returns the data from the class' data source corresponding to a *buffered* version of the input extent.
|
||||
Buffering is done by `self.padding` number of units (in terms of the source coordinate system).
|
||||
|
||||
Args:
|
||||
extent (dict): A geographic extent formatted as a dictionary with the following keys: xmin, xmax, ymin, ymax, crs
|
||||
|
||||
Returns:
|
||||
output_raster (InMemoryRaster): A raster cropped to a *buffered* version of the input extent.
|
||||
"""
|
||||
raise NotImplementedError()
|
||||
|
||||
@abstractmethod
|
||||
def get_area_from_shape_by_extent(self, extent, shape_layer):
|
||||
raise NotImplementedError()
|
||||
@abc.abstractmethod
|
||||
def get_data_from_geometry(self, geometry):
|
||||
"""Returns the data from the class' raster source corresponding to the input `geometry` without any buffering applied. Note that `geometry` is expected
|
||||
to be a GeoJSON polygon (as a dictionary) in the EPSG:4326 coordinate system.
|
||||
|
||||
@abstractmethod
|
||||
def get_data_from_shape(self, shape):
|
||||
Args:
|
||||
geometry (dict): A polygon in GeoJSON format describing the boundary to crop the input raster to
|
||||
|
||||
Returns:
|
||||
output_raster (InMemoryRaster): A raster cropped to the outline of `geometry`
|
||||
"""
|
||||
raise NotImplementedError()
|
||||
|
|
@ -13,114 +13,58 @@ LOGGER = logging.getLogger("server")
|
|||
from . import ROOT_DIR
|
||||
from .DataLoader import DataLoaderCustom, DataLoaderUSALayer, DataLoaderBasemap
|
||||
|
||||
|
||||
def get_area_from_geometry(geom, src_crs="epsg:4326"):
|
||||
if geom["type"] == "Polygon":
|
||||
lon, lat = geom["coordinates"][0][0]
|
||||
elif geom["type"] == "MultiPolygon":
|
||||
lon, lat = geom["coordinates"][0][0][0]
|
||||
else:
|
||||
raise ValueError("Polygons and MultiPolygons only")
|
||||
|
||||
zone_number = utm.latlon_to_zone_number(lat, lon)
|
||||
hemisphere = "+north" if lat > 0 else "+south"
|
||||
dest_crs = "+proj=utm +zone=%d %s +datum=WGS84 +units=m +no_defs" % (zone_number, hemisphere)
|
||||
projected_geom = fiona.transform.transform_geom(src_crs, dest_crs, geom)
|
||||
area = shapely.geometry.shape(projected_geom).area / 1000000.0 # we calculate the area in square meters then convert to square kilometers
|
||||
return area
|
||||
|
||||
def _load_geojson_as_list(fn):
|
||||
''' Takes a geojson file as input and outputs a list of shapely `shape` objects in that file and their corresponding areas in km^2.
|
||||
|
||||
We calculate area here by re-projecting the shape into its local UTM zone, converting it to a shapely `shape`, then using the `.area` property.
|
||||
'''
|
||||
shapes = []
|
||||
areas = []
|
||||
crs = None
|
||||
with fiona.open(fn) as f:
|
||||
src_crs = f.crs
|
||||
for row in f:
|
||||
geom = row["geometry"]
|
||||
|
||||
area = get_area_from_geometry(geom, src_crs)
|
||||
areas.append(area)
|
||||
|
||||
shape = shapely.geometry.shape(geom)
|
||||
shapes.append(shape)
|
||||
return shapes, areas, src_crs
|
||||
|
||||
|
||||
def _load_dataset(dataset):
|
||||
# Step 1: load the shape layers
|
||||
shape_layers = {}
|
||||
if dataset["shapeLayers"] is not None:
|
||||
for shape_layer in dataset["shapeLayers"]:
|
||||
fn = shape_layer["shapesFn"]
|
||||
if os.path.exists(fn):
|
||||
shapes, areas, crs = _load_geojson_as_list(fn)
|
||||
|
||||
shape_layer["geoms"] = shapes
|
||||
shape_layer["areas"] = areas
|
||||
shape_layer["crs"] = crs["init"] # TODO: will this break with fiona version; I think `.crs` will turn into a PyProj object
|
||||
shape_layers[shape_layer["name"]] = shape_layer
|
||||
else:
|
||||
LOGGER.warning("Step 1 failed in loading dataset {}".format(dataset["metadata"]["displayName"]))
|
||||
# TODO: check that the displayName field is present
|
||||
return False # TODO: maybe we should make these errors more descriptive (explain why we can't load a dataset)
|
||||
|
||||
# Step 2: make sure the dataLayer exists
|
||||
# Step 1: make sure the dataLayer exists
|
||||
if dataset["dataLayer"]["type"] == "CUSTOM":
|
||||
fn = dataset["dataLayer"]["path"]
|
||||
if not os.path.exists(fn):
|
||||
LOGGER.warning("Step 2 failed in loading dataset {}".format(dataset["metadata"]["displayName"]))
|
||||
return False # TODO: maybe we should make these errors more descriptive (explain why we can't load a dataset)
|
||||
|
||||
# Step 3: setup the appropriate DatasetLoader
|
||||
# Step 2: setup the appropriate DatasetLoader
|
||||
if dataset["dataLayer"]["type"] == "CUSTOM":
|
||||
data_loader = DataLoaderCustom(dataset["dataLayer"]["path"], shape_layers, dataset["dataLayer"]["padding"])
|
||||
data_loader = DataLoaderCustom(**dataset["dataLayer"])
|
||||
elif dataset["dataLayer"]["type"] == "USA_LAYER":
|
||||
data_loader = DataLoaderUSALayer(shape_layers, dataset["dataLayer"]["padding"])
|
||||
data_loader = DataLoaderUSALayer(dataset["dataLayer"]["padding"])
|
||||
elif dataset["dataLayer"]["type"] == "BASEMAP":
|
||||
data_loader = DataLoaderBasemap(dataset["dataLayer"]["path"], dataset["dataLayer"]["padding"])
|
||||
data_loader = DataLoaderBasemap(dataset["dataLayer"]["padding"])
|
||||
else:
|
||||
LOGGER.warning("Step 3 failed in loading dataset {}".format(dataset["metadata"]["displayName"]))
|
||||
return False # TODO: maybe we should make these errors more descriptive (explain why we can't load a dataset)
|
||||
|
||||
return {
|
||||
"data_loader": data_loader,
|
||||
"shape_layers": shape_layers,
|
||||
}
|
||||
return data_loader
|
||||
|
||||
def load_datasets():
|
||||
"""Returns a dictionary of key:value where keys are dataset names (from "datasets.json" / "datasets.mine.json") and values are instances of classes that extend DataLoaderAbstract
|
||||
"""
|
||||
datasets = dict()
|
||||
|
||||
dataset_json = json.load(open(os.path.join(ROOT_DIR, "datasets.json"),"r"))
|
||||
for key, dataset in dataset_json.items():
|
||||
dataset_object = _load_dataset(dataset)
|
||||
data_loader = _load_dataset(dataset)
|
||||
|
||||
if dataset_object is False:
|
||||
if data_loader is False:
|
||||
LOGGER.warning("Files are missing, we will not be able to serve the following dataset: '%s'" % (key))
|
||||
else:
|
||||
datasets[key] = dataset_object
|
||||
datasets[key] = data_loader
|
||||
|
||||
if os.path.exists(os.path.join(ROOT_DIR, "datasets.mine.json")):
|
||||
dataset_json = json.load(open(os.path.join(ROOT_DIR, "datasets.mine.json"),"r"))
|
||||
for key, dataset in dataset_json.items():
|
||||
|
||||
if key not in datasets:
|
||||
dataset_object = _load_dataset(dataset)
|
||||
data_loader = _load_dataset(dataset)
|
||||
|
||||
if dataset_object is False:
|
||||
if data_loader is False:
|
||||
LOGGER.warning("Files are missing, we will not be able to serve the following dataset: '%s'" % (key))
|
||||
else:
|
||||
datasets[key] = dataset_object
|
||||
datasets[key] = data_loader
|
||||
else:
|
||||
LOGGER.warning("There is a conflicting dataset key in datasets.mine.json, skipping.")
|
||||
|
||||
return datasets
|
||||
|
||||
def is_valid_dataset(dataset_key):
|
||||
|
||||
dataset_json = json.load(open(os.path.join(ROOT_DIR, "datasets.json"), "r"))
|
||||
if os.path.exists(os.path.join(ROOT_DIR, "datasets.mine.json")):
|
||||
dataset_mine_json = json.load(open(os.path.join(ROOT_DIR, "datasets.mine.json"), "r"))
|
||||
|
|
|
|||
|
|
@ -85,5 +85,32 @@
|
|||
}
|
||||
],
|
||||
"validModels": ["naip_demo"]
|
||||
},
|
||||
"naip_all": {
|
||||
"metadata": {
|
||||
"displayName": "NAIP USA",
|
||||
"locationName": null
|
||||
},
|
||||
"dataLayer": {
|
||||
"type": "USA_LAYER",
|
||||
"padding": 20,
|
||||
"resolution": 1
|
||||
},
|
||||
"basemapLayers": [
|
||||
{
|
||||
"layerName": "NAIP 2017 Imagery",
|
||||
"initialZoom": 11,
|
||||
"url": "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
|
||||
"initialLocation": [38.477018, -75.402312],
|
||||
"args": {
|
||||
"attribution": "Tiles © Esri — Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community",
|
||||
"maxNativeZoom": 16,
|
||||
"maxZoom": 20,
|
||||
"minZoom": 4
|
||||
}
|
||||
}
|
||||
],
|
||||
"shapeLayers": null,
|
||||
"validModels": ["naip_demo"]
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Загрузка…
Ссылка в новой задаче