diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 732992f..299b0a1 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -1,7 +1,11 @@ +import os from abc import abstractmethod -from typing import Union, Tuple, List +from typing import Union, Tuple +from uuid import uuid4 +from osgeo import gdal import ee +import boto3 from dask.diagnostics import ProgressBar from ee import ImageCollection from geocube.api.core import make_geocube @@ -52,6 +56,35 @@ def groupby(self, zones, layer=None): """ return LayerGroupBy(self.aggregate, zones, layer, self.masks) + def write(self, bbox, output_path, tile_degrees=None): + """ + Write the layer to a path. Does not apply masks. + + :param bbox: (min x, min y, max x, max y) + :param output_path: local or s3 path to output to + :param tile_degrees: optional param to tile the results into multiple files with a VRT. + Degrees to tile by. `output_path` should be a folder path to store the tiles. + :return: + """ + + if tile_degrees is not None: + tiles = create_fishnet_grid(*bbox, tile_degrees) + + if not os.path.exists(output_path): + os.makedirs(output_path) + + file_names = [] + for tile in tiles["geometry"]: + data = self.aggregate.get_data(tile.bounds) + + file_name = f"{output_path}/{uuid4()}.tif" + file_names.append(file_name) + + write_layer(file_name, data) + else: + data = self.aggregate.get_data(bbox) + write_layer(output_path, data) + class LayerGroupBy: def __init__(self, aggregate, zones, layer=None, masks=[]): @@ -268,3 +301,21 @@ def get_image_collection( return data + +def write_layer(path, data): + if isinstance(data, xr.DataArray): + # for rasters, need to write to locally first then copy to cloud storage + if path.startswith("s3://"): + tmp_path = f"{uuid4()}.tif" + data.rio.to_raster(raster_path=tmp_path, driver="COG") + + s3 = boto3.client('s3') + s3.upload_file(tmp_path, path.split('/')[2], '/'.join(path.split('/')[3:])) + + os.remove(tmp_path) + else: + data.rio.to_raster(raster_path=path, driver="COG") + elif isinstance(data, gpd.GeoDataFrame): + data.to_file(path, driver="GeoJSON") + else: + raise NotImplementedError("Can only write DataArray or GeoDataFrame") diff --git a/environment.yml b/environment.yml index 5848dbc..436dbbc 100644 --- a/environment.yml +++ b/environment.yml @@ -17,6 +17,7 @@ dependencies: - dask[complete]=2023.11.0 - matplotlib=3.8.2 - jupyterlab=4.0.10 + - s3fs=2024.5.0 - geemap=0.32.0 - pip=23.3.1 - pip: diff --git a/notebooks/tutorial/get layers.ipynb b/notebooks/tutorial/get layers.ipynb index 883d107..e67ce8a 100644 --- a/notebooks/tutorial/get layers.ipynb +++ b/notebooks/tutorial/get layers.ipynb @@ -86,7 +86,24 @@ }, { "cell_type": "markdown", - "id": "6280fc2f", + "id": "5c506fca", + "metadata": {}, + "source": [] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "286103d4", + "metadata": {}, + "outputs": [], + "source": [ + "import boto3" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "d4f01f2f-5164-4a05-9997-f70b2abe6b37", "metadata": {}, "source": [ "In most cases you shouldn't need to set these. GEE will open your browser to authenticate. " @@ -240,7 +257,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 5, "id": "5eb883fd-c533-47b7-8735-65393afca89d", "metadata": {}, "outputs": [], @@ -250,7 +267,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 6, "id": "0703cfe0-201d-4cae-ab43-7a4d0e5082cf", "metadata": {}, "outputs": [ @@ -258,8 +275,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "Extracting tree cover layer:\n", - "[########################################] | 100% Completed | 19.58 s\n" + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 9.24 ss\n" ] }, { @@ -628,40 +645,40 @@ " stroke: currentColor;\n", " fill: currentColor;\n", "}\n", - "
<xarray.DataArray 'ttc' (y: 2878, x: 3721)>\n",
-       "array([[ nan,  nan,  nan, ..., 100., 100.,  90.],\n",
-       "       [ nan,  nan,  nan, ...,  90.,  90.,  90.],\n",
+       "
<xarray.DataArray 'ttc' (y: 2878, x: 3722)>\n",
+       "array([[ nan,  nan,  nan, ..., 100.,  90.,  90.],\n",
        "       [ nan,  nan,  nan, ...,  90.,  90.,  90.],\n",
+       "       [ nan,  nan,  nan, ...,  90.,  90.,  80.],\n",
        "       ...,\n",
        "       [ nan,  nan,  nan, ...,  nan,  nan,  nan],\n",
        "       [ nan,  nan,  nan, ...,  nan,  nan,  nan],\n",
        "       [ nan,  nan,  nan, ...,  nan,  nan,  nan]])\n",
        "Coordinates:\n",
-       "    time     int32 0\n",
+       "    time     int64 0\n",
        "  * x        (x) float32 -38.65 -38.65 -38.65 -38.65 ... -38.3 -38.3 -38.3 -38.3\n",
        "  * y        (y) float32 -12.76 -12.76 -12.76 -12.76 ... -13.02 -13.02 -13.02\n",
        "Attributes:\n",
        "    id:             ttc\n",
        "    data_type:      {'type': 'PixelType', 'precision': 'double', 'min': 0, 'm...\n",
        "    crs:            EPSG:4326\n",
-       "    crs_transform:  [1, 0, 0, 0, 1, 0]
  • id :
    ttc
    data_type :
    {'type': 'PixelType', 'precision': 'double', 'min': 0, 'max': 255}
    crs :
    EPSG:4326
    crs_transform :
    [1, 0, 0, 0, 1, 0]
    • " ], "text/plain": [ - "\n", - "array([[ nan, nan, nan, ..., 100., 100., 90.],\n", - " [ nan, nan, nan, ..., 90., 90., 90.],\n", + "\n", + "array([[ nan, nan, nan, ..., 100., 90., 90.],\n", " [ nan, nan, nan, ..., 90., 90., 90.],\n", + " [ nan, nan, nan, ..., 90., 90., 80.],\n", " ...,\n", " [ nan, nan, nan, ..., nan, nan, nan],\n", " [ nan, nan, nan, ..., nan, nan, nan],\n", " [ nan, nan, nan, ..., nan, nan, nan]])\n", "Coordinates:\n", - " time int32 0\n", + " time int64 0\n", " * x (x) float32 -38.65 -38.65 -38.65 -38.65 ... -38.3 -38.3 -38.3 -38.3\n", " * y (y) float32 -12.76 -12.76 -12.76 -12.76 ... -13.02 -13.02 -13.02\n", "Attributes:\n", @@ -692,7 +709,7 @@ " crs_transform: [1, 0, 0, 0, 1, 0]" ] }, - "execution_count": 8, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } @@ -703,6 +720,75 @@ "city_TreeCover" ] }, + { + "cell_type": "code", + "execution_count": 7, + "id": "8c2477dc", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 8.91 ss\n" + ] + } + ], + "source": [ + "city_TreeCover = TreeCover().write(city_gdf.total_bounds, \"output.tif\")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "cb58b256", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.63 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.52 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.53 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.38 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.35 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 2.47 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.53 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.45 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.46 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.54 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.69 sms\n", + "Extracting layer tree cover from Google Earth Engine:\n", + "[########################################] | 100% Completed | 1.59 sms\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/usr/local/anaconda3/envs/cities-cif/lib/python3.10/site-packages/osgeo/gdal.py:287: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.\n", + " warnings.warn(\n", + "ERROR 4: Failed to open output to write.\n" + ] + } + ], + "source": [ + "city_TreeCover = TreeCover().write(city_gdf.total_bounds, \"output\", tile_degrees=0.1)" + ] + }, { "cell_type": "code", "execution_count": 9, diff --git a/setup.py b/setup.py index bff48fe..c2d24d8 100644 --- a/setup.py +++ b/setup.py @@ -20,5 +20,6 @@ "utm", "osmnx", "geopandas", + "s3fs", ], )