From 6ce6f12591c885e3c228eeab00b72fce7ebdc002 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 29 Nov 2023 19:15:08 +0000 Subject: [PATCH 01/58] Bump pypa/gh-action-pypi-publish from 1.8.10 to 1.8.11 Bumps [pypa/gh-action-pypi-publish](https://github.com/pypa/gh-action-pypi-publish) from 1.8.10 to 1.8.11. - [Release notes](https://github.com/pypa/gh-action-pypi-publish/releases) - [Commits](https://github.com/pypa/gh-action-pypi-publish/compare/v1.8.10...v1.8.11) --- updated-dependencies: - dependency-name: pypa/gh-action-pypi-publish dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] --- .github/workflows/distribute.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/distribute.yml b/.github/workflows/distribute.yml index 034c4ac0..bad4d8b1 100644 --- a/.github/workflows/distribute.yml +++ b/.github/workflows/distribute.yml @@ -29,7 +29,7 @@ jobs: python -m build - name: upload to PyPI.org - uses: pypa/gh-action-pypi-publish@v1.8.10 + uses: pypa/gh-action-pypi-publish@v1.8.11 with: user: __token__ password: ${{ secrets.TOOLS_PYPI_PAK }} From 4e0e439a49ca76c75f2a4b9ec4bd0efea8fa64b4 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 6 Dec 2023 19:50:25 +0000 Subject: [PATCH 02/58] Bump actions/setup-python from 4 to 5 Bumps [actions/setup-python](https://github.com/actions/setup-python) from 4 to 5. - [Release notes](https://github.com/actions/setup-python/releases) - [Commits](https://github.com/actions/setup-python/compare/v4...v5) --- updated-dependencies: - dependency-name: actions/setup-python dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/static-analysis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/static-analysis.yml b/.github/workflows/static-analysis.yml index 3c258e53..723ac27b 100644 --- a/.github/workflows/static-analysis.yml +++ b/.github/workflows/static-analysis.yml @@ -8,7 +8,7 @@ jobs: steps: - uses: actions/checkout@v4 - - uses: actions/setup-python@v4 + - uses: actions/setup-python@v5 with: python-version: "3.10" From ab4e51103558c2e98d487be57492cf60ceaf820d Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 21 Feb 2024 10:32:24 -0600 Subject: [PATCH 03/58] osm water mask creation --- .../watermasking/generate_osm_tiles.py | 227 ++++++++++++++++++ src/asf_tools/watermasking/utils.py | 19 ++ 2 files changed, 246 insertions(+) create mode 100644 src/asf_tools/watermasking/generate_osm_tiles.py create mode 100644 src/asf_tools/watermasking/utils.py diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py new file mode 100644 index 00000000..cde6644c --- /dev/null +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -0,0 +1,227 @@ +import argparse +import os +import time + +import geopandas as gpd +import numpy as np +from osgeo import gdal + +from hyp3_testing.watermasking.utils import lat_lon_to_tile_string + + +def process_pbf(planet_file: str, output_file: str): + """Process the global PBF file so that it only contains water features. + + Args: + planet_file: The path to the OSM Planet PBF file. + output_file: The desired path of the processed PBF file. + """ + + natural_file = 'planet_natural.pbf' + waterways_file = 'planet_waterways.pbf' + reservoirs_file = 'planet_reservoirs.pbf' + + natural_water_command = f'osmium tags-filter -o {natural_file} {planet_file} wr/natural=water' + waterways_command = f'osmium tags-filter -o {waterways_file} {planet_file} waterway="*"' + reservoirs_command = f'osmium tags-filter -o {reservoirs_file} {planet_file} landuse=reservoir' + merge_command = f'osmium merge {natural_file} {waterways_file} {reservoirs_file} -o {output_file}' + + os.system(natural_water_command) + os.system(waterways_command) + os.system(reservoirs_command) + os.system(merge_command) + + +def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_height_deg): + """Process and crop OSM ocean polygons into a tif tile. + + Args: + ocean_polygons_path: The path to the global ocean polygons file from OSM. + lat: The minimum latitude of the tile. + lon: The minimum longitude of the tile. + tile_width_deg: The width of the tile in degrees. + tile_height_deg: The height of the tile in degrees. + """ + + tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') + tile_tif = 'tiles/' + tile + '.tif' + + xmin, xmax, ymin, ymax = lon, lon+tile_width_deg, lat, lat+tile_height_deg + pixel_size_x = 0.00009009009 # 10m * 2 at the equator. + pixel_size_y = 0.00009009009 + + clipped_polygons_path = tile + '.shp' + command = f'ogr2ogr -clipsrc {xmin} {ymin} {xmax} {ymax} {clipped_polygons_path} {ocean_polygons_path}' + os.system(command) + + gdal.Rasterize( + tile_tif, + clipped_polygons_path, + xRes=pixel_size_x, + yRes=pixel_size_y, + burnValues=1, + outputBounds=[xmin, ymin, xmax, ymax], + outputType=gdal.GDT_Byte, + creationOptions=['COMPRESS=LZW'] + ) + + clipped_dbf = tile + '.dbf' + clipped_shx = tile + '.shx' + clipped_prj = tile + '.prj' + clipped_cpg = tile + '.cpg' + os.system(f'rm {clipped_polygons_path} {clipped_dbf} {clipped_shx} {clipped_prj} {clipped_cpg} > /dev/null 2>&1') + + +def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg): + """Rasterize a water tile from the processed global PBF file. + + Args: + water_file: The path to the processed global PBF file. + lat: The minimum latitude of the tile. + lon: The minimum longitude of the tile. + tile_width_deg: The desired width of the tile in degrees. + tile_height_deg: The desired height of the tile in degrees. + """ + + tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') + tile_pbf = tile+ '.osm.pbf' + tile_tif = tile + '.tif' + tile_shp = tile + '.shp' + tile_geojson = tile + '.geojson' + + # Extract tile from the main pbf, then convert it to a tif. + os.system(f'osmium extract -s smart -S tags=natural=water --bbox {lon},{lat},{lon+tile_width_deg},{lat+tile_height_deg} {water_file} -o {tile_pbf}') + os.system(f'osmium export --geometry-types="polygon" {tile_pbf} -o {tile_geojson}') + + # Islands and Islets can be members of the water features, so they must be removed. + water_gdf = gpd.read_file(tile_geojson, engine='pyogrio') + try: + water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'island'].index) + water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'islet'].index) + except: + pass + water_gdf.to_file(tile_shp, mode='w', engine='pyogrio') + water_gdf = None + + xmin, xmax, ymin, ymax = lon, lon+tile_width_deg, lat, lat+tile_height_deg + pixel_size_x = 0.00009009009 # 10m at the equator. + pixel_size_y = 0.00009009009 + + gdal.Rasterize( + tile_tif, + tile_shp, + xRes=pixel_size_x, + yRes=pixel_size_y, + burnValues=1, + outputBounds=[xmin, ymin, xmax, ymax], + outputType=gdal.GDT_Byte, + creationOptions=['COMPRESS=LZW'] + ) + + tile_dbf = tile + '.dbf' + tile_cpg = tile + '.cpg' + tile_prj = tile + '.prj' + tile_shx = tile + '.shx' + os.system(f'rm {tile_geojson} {tile_pbf} {tile_shp} {tile_dbf} {tile_prj} {tile_cpg} {tile_shx}') + + +def merge_tiles(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_to_cog: bool = False): + """Merge the interior water tiles and ocean water tiles. + + Args: + interior_tile_dir: The path to the directory containing the interior water tiles. + ocean_tile_dir: The path to the directory containing the ocean water tiles. + merged_tile_dir: The path to the directory containing the merged water tiles. + """ + index = 0 + num_tiles = len([f for f in os.listdir(internal_tile_dir) if f.endswith('tif')]) + for filename in os.listdir(internal_tile_dir): + if filename.endswith('.tif'): + try: + start_time = time.time() + + internal_tile = internal_tile_dir + filename + external_tile = ocean_tile_dir + filename + output_tile = finished_tile_dir + filename + command = f'gdal_calc.py -A {internal_tile} -B {external_tile} --format GTiff -co NUM_THREADS=all_cpus -co COMPRESS=LZW -co TILED=YES -co PREDICTOR=YES --outfile {output_tile} --calc "logical_or(A, B)"' + os.system(command) + + if translate_to_cog: + command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus' + + end_time = time.time() + total_time = end_time - start_time + + print(f'Elapsed Time: {total_time}(s)') + print(f'Completed {index} of {num_tiles}') + + index += 1 + except Exception as e: + print(f'Caught error while processing {filename}. Continuing anyways...\n{e}') + + +def setup_directories(): + dirs = ['interior_tiles/', 'ocean_tiles/', 'merged_tiles/', 'finished_tiles/'] + for dir in dirs: + try: + os.mkdir(dir) + except FileExistsError as e: + print(f'{dir} already exists. Skipping...') + + +def main(): + + parser = argparse.ArgumentParser( + prog='generate_osm_tiles.py', + description='Main script for creating a tiled watermask dataset from OSM data.' + ) + + parser.add_argument('--planet-file-path', help='The path to the global planet.pbf file.', default='planet.pbf') + parser.add_argument('--ocean-polygons-path', help='The path to the global OSM ocean polygons.', default='water_polygons.shp') + parser.add_argument('--lat-begin', help='The minimum latitude of the dataset.', default=-85) + parser.add_argument('--lat-end', help='The maximum latitude of the dataset.', default=85) + parser.add_argument('--lon-begin', help='The minimum longitude of the dataset.', default=-180) + parser.add_argument('--lon-end', help='The maximum longitude of the dataset.', default=180) + parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) + parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) + + args = parser.parse_args() + + # Setup Directories + interior_tile_dir = 'interior_tiles/' + ocean_tile_dir = 'ocean_tiles/' + finished_tile_dir = 'tiles/' + + dirs = [interior_tile_dir, ocean_tile_dir, finished_tile_dir] + for dir in dirs: + try: + os.mkdir(dir) + except FileExistsError as e: + print(f'{dir} already exists. Skipping...') + + # Process the PBF into one that contains only water. + processed_pbf_path = 'planet_processed.pbf' + process_pbf(args.planet_file_path, processed_pbf_path) + + # Extract tiles from the processed and ocean files. + lat_range = range(args.lat_begin, args.lat_end, args.tile_height) + lon_range = range(args.lon_begin, args.lon_end, args.tile_width) + for lat in lat_range: + for lon in lon_range: + tile_name = lat_lon_to_tile_string(lat, lon, is_worldcover=False) + try: + start_time = time.time() + extract_water(processed_pbf_path, lat, lon, args.tile_width, args.tile_height) + process_ocean_tiles(args.ocean_polygons_path, lat, lon, args.tile_width, args.tile_height) + end_time = time.time() + total_time = end_time - start_time #seconds + print(f'Finished initial creation of {tile_name} in {total_time}(s).') + except Exception as e: + print(f'Caught error while processing {tile_name}. Continuing anyways...\n{e}') + + # Merge the extracted processed/ocean tiles. + merge_tiles(interior_tile_dir, ocean_tile_dir) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py new file mode 100644 index 00000000..1935eaad --- /dev/null +++ b/src/asf_tools/watermasking/utils.py @@ -0,0 +1,19 @@ +import os +import time + +import numpy as np + +from osgeo import gdal + + +def lat_lon_to_tile_string(lat, lon, is_worldcover: bool = False, postfix: str ='.tif'): + prefixes = ['N', 'S', 'E', 'W'] if is_worldcover else ['n', 's', 'e', 'w'] + if lat >= 0: + lat_part = prefixes[0] + str(int(lat)).zfill(2) + else: + lat_part = prefixes[1] + str(int(np.abs(lat))).zfill(2) + if lon >= 0: + lon_part = prefixes[2] + str(int(lon)).zfill(3) + else: + lon_part = prefixes[3] + str(int(np.abs(lon))).zfill(3) + return lat_part + lon_part + postfix \ No newline at end of file From 8eb65f0158ea93e0ae5060fc7545bc9e995aa64f Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 21 Feb 2024 10:56:45 -0600 Subject: [PATCH 04/58] added function for removing temp files --- .../watermasking/generate_osm_tiles.py | 31 ++++++++++++------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index cde6644c..1f60d50e 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -6,7 +6,7 @@ import numpy as np from osgeo import gdal -from hyp3_testing.watermasking.utils import lat_lon_to_tile_string +from asf_tools.watermasking.utils import lat_lon_to_tile_string def process_pbf(planet_file: str, output_file: str): @@ -32,6 +32,19 @@ def process_pbf(planet_file: str, output_file: str): os.system(merge_command) +def remove_temp_files(temp_files: list): + """Remove each file in a list of files. + + Args: + temp_files: The list of temporary files to remove. + """ + for file in temp_files: + try: + os.remove(file) + except: + print('Error removing temporary file. Skipping it...') + + def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_height_deg): """Process and crop OSM ocean polygons into a tif tile. @@ -65,11 +78,8 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig creationOptions=['COMPRESS=LZW'] ) - clipped_dbf = tile + '.dbf' - clipped_shx = tile + '.shx' - clipped_prj = tile + '.prj' - clipped_cpg = tile + '.cpg' - os.system(f'rm {clipped_polygons_path} {clipped_dbf} {clipped_shx} {clipped_prj} {clipped_cpg} > /dev/null 2>&1') + temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx'] + remove_temp_files(temp_files) def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg): @@ -99,6 +109,7 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg): water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'island'].index) water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'islet'].index) except: + # When there are no islands to remove, an error will be occur, but we don't care about it. pass water_gdf.to_file(tile_shp, mode='w', engine='pyogrio') water_gdf = None @@ -118,11 +129,8 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg): creationOptions=['COMPRESS=LZW'] ) - tile_dbf = tile + '.dbf' - tile_cpg = tile + '.cpg' - tile_prj = tile + '.prj' - tile_shx = tile + '.shx' - os.system(f'rm {tile_geojson} {tile_pbf} {tile_shp} {tile_dbf} {tile_prj} {tile_cpg} {tile_shx}') + temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx'] + remove_temp_files(temp_files) def merge_tiles(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_to_cog: bool = False): @@ -161,6 +169,7 @@ def merge_tiles(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_ def setup_directories(): + """Setup the directories necessary for running the script.""" dirs = ['interior_tiles/', 'ocean_tiles/', 'merged_tiles/', 'finished_tiles/'] for dir in dirs: try: From e67e03dd9487e1ad86e30f3fcd4d7fd5d1af518c Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 21 Feb 2024 10:58:30 -0600 Subject: [PATCH 05/58] updated the dependencies --- environment.yml | 1 + pyproject.toml | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 7646b3e7..b26f47ef 100644 --- a/environment.yml +++ b/environment.yml @@ -23,6 +23,7 @@ dependencies: - boto3 - fiona - gdal>=3.7 + - geopandas - numpy - pysheds>=0.3 - rasterio diff --git a/pyproject.toml b/pyproject.toml index 4a5f0959..cf9f2da4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,14 +25,14 @@ dependencies = [ "astropy", "fiona", "gdal>=3.3", + "geopandas", "numpy", "pysheds>=0.3", "rasterio", "scikit-fuzzy", "scikit-image", "scipy", - "shapely", - "tqdm", + "shapely" ] dynamic = ["version"] From 662ec27be51cf24a795c338412fb585935ece9df Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 21 Feb 2024 12:14:25 -0600 Subject: [PATCH 06/58] worldcover code --- .../watermasking/generate_worldcover_tiles.py | 183 ++++++++++++++++++ src/asf_tools/watermasking/utils.py | 42 +++- 2 files changed, 224 insertions(+), 1 deletion(-) create mode 100644 src/asf_tools/watermasking/generate_worldcover_tiles.py diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py new file mode 100644 index 00000000..94bc8571 --- /dev/null +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -0,0 +1,183 @@ +import argparse +import os +import time + +import numpy as np + +from osgeo import gdal + +from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, merge_tiles + +TILE_DIR = 'worldcover_tiles_unfinished/' +CROPPED_TILE_DIR = 'worldcover_tiles/' +FILENAME_POSTFIX = '.tif' +WORLDCOVER_TILE_SIZE = 3 + + +def tile_preprocessing(tile_dir, min_lat, max_lat): + """The worldcover tiles have lots of unnecessary classes, these need to be removed first. + Note: make a back-up copy of this directory. + + Args: + tile_dir: The directory containing all of the worldcover tiles. + """ + + filenames = [f for f in os.listdir(tile_dir) if f.endswith('.tif')] + filter = lambda filename: (int(filename.split('_')[5][1:3]) >= min_lat) and (int(filename.split('_')[5][1:3]) <= max_lat) + filenames_filtered = [f for f in filenames if filter(f)] + + index = 0 + num_tiles = len(filenames_filtered) + for filename in filenames_filtered: + start_time = time.time() + + dst_filename = filename.split('_')[5] + '.tif' + + print(f'Processing: {filename} --- {dst_filename} -- {index} of {num_tiles}') + + src_ds = gdal.Open(filename) + src_band = src_ds.GetRasterBand(1) + src_arr = src_band.ReadAsArray() + + not_water = np.logical_and(src_arr != 80, src_arr != 0) + water_arr = np.ones(src_arr.shape) + water_arr[not_water] = 0 + + driver = gdal.GetDriverByName('GTiff') + + dst_ds = driver.Create(dst_filename, water_arr.shape[0], water_arr.shape[1], 1, gdal.GDT_Byte, options=['COMPRESS=LZW', 'TILED=YES']) + dst_ds.SetGeoTransform(src_ds.GetGeoTransform()) + dst_ds.SetProjection(src_ds.GetProjection()) + dst_band = dst_ds.GetRasterBand(1) + dst_band.WriteArray(water_arr) + dst_band.FlushCache() + + del dst_ds + del src_ds + + end_time = time.time() + total_time = end_time - start_time # seconds + + print(f'Processing {dst_filename} took {total_time} seconds.') + + index += 1 + + +def get_tiles(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): + """Get a list of the worldcover tile locations necessary to fully cover an OSM tile. + + Args: + osm_tile_coord: The lower left corner coordinate (lat, lon) of the desired OSM tile. + wc_deg: The width/height of the Worldcover tiles in degrees. + osm_deg: The width/height of the OSM tiles in degrees. + + Returns: + tiles: A list of the lower left corner coordinates of the Worldcover tiles that overlap the OSM tile. + """ + + osm_lat = osm_tile_coord[0] + osm_lon = osm_tile_coord[1] + + min_lat = osm_lat - (osm_lat % wc_deg) + max_lat = osm_lat + osm_deg + min_lon = osm_lon - (osm_lon % wc_deg) + max_lon = osm_lon + osm_deg + + lats = range(min_lat, max_lat, wc_deg) + lons = range(min_lon, max_lon, wc_deg) + + tiles = [] + for lat in lats: + for lon in lons: + tiles.append((lat, lon)) + + return tiles + + +def lat_lon_to_filenames(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): + """Get a list of the Worldcover tile filenames that are necessary to overlap an OSM tile. + + Args: + osm_tile: The lower left corner (lat, lon) of the desired OSM tile. + wc_deg: The width of the Worldcover tiles in degrees. + osm_deg: The width of the OSM tiles in degrees. + + Returns: + filenames: The list of Worldcover filenames. + """ + filenames = [] + tiles = get_tiles(osm_tile_coord, wc_deg, osm_deg) + for tile in tiles: + filenames.append(lat_lon_to_tile_string(tile[0], tile[1], is_worldcover=True)) + return filenames + + +def crop_tile(tile): + """Crop the merged tiles""" + try: + ref_image = TILE_DIR + tile + pixel_size = gdal.Warp('tmp_px_size.tif', ref_image, dstSRS='EPSG:4326').GetGeoTransform()[1] + shapefile_command = ' '.join(['gdaltindex', 'tmp.shp', ref_image]) + os.system(shapefile_command) + out_filename = CROPPED_TILE_DIR + tile + gdal.Warp( + out_filename, + tile, + cutlineDSName='tmp.shp', + cropToCutline=True, + xRes=pixel_size, + yRes=pixel_size, + targetAlignedPixels=True, + dstSRS='EPSG:4326', + format='COG' + ) + remove_temp_files(['tmp_px_size.tif', 'tmp.shp']) + except Exception as e: # When a tile fails to process, we don't necessarily want the program to stop. + print(f'Could not process {tile}. Continuing...') + index += 1 + + +def build_dataset(lat_range, lon_range, worldcover_degrees, osm_degrees): + for lat in lat_range: + for lon in lon_range: + start_time = time.time() + tile_filename = TILE_DIR + lat_lon_to_tile_string(lat, lon, is_worldcover=False) + worldcover_tiles = lat_lon_to_filenames(lat, lon, worldcover_degrees, osm_degrees) + print(f'Processing: {tile_filename} {worldcover_tiles}') + merge_tiles(worldcover_tiles, tile_filename) + end_time = time.time() + total_time = end_time - start_time + print(f'Time Elapsed: {total_time}s') + + +def main(): + parser = argparse.ArgumentParser( + prog='generate_worldcover_tiles.py', + description='Main script for creating a tiled watermask dataset from the ESA WorldCover dataset.' + ) + parser.add_argument('--worldcover-tiles-dir', help='The path to the directory containing the worldcover tifs.') + parser.add_argument('--lat-begin', help='The minimum latitude of the dataset.', default=-85, required=True) + parser.add_argument('--lat-end', help='The maximum latitude of the dataset.', default=85) + parser.add_argument('--lon-begin', help='The minimum longitude of the dataset.', default=-180) + parser.add_argument('--lon-end', help='The maximum longitude of the dataset.', default=180) + parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) + parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) + + args = parser.parse_args() + + for dir in [TILE_DIR, CROPPED_TILE_DIR]: + try: + os.mkdir(dir) + except FileExistsError as e: + print(f'{dir} already exists. Skipping...') + + tile_preprocessing(args.worldcover_tiles_dir, args.lat_begin, args.lat_end) + + lat_range = range(args.lat_begin, args.lat_end, args.tile_width) + lon_range = range(args.lon_begin, args.lon_end, args.tile_heigth) + + build_dataset(lat_range, lon_range, worldcover_degrees=WORLDCOVER_TILE_SIZE, osm_degrees=args.tile_width) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py index 1935eaad..15b56bd4 100644 --- a/src/asf_tools/watermasking/utils.py +++ b/src/asf_tools/watermasking/utils.py @@ -7,6 +7,17 @@ def lat_lon_to_tile_string(lat, lon, is_worldcover: bool = False, postfix: str ='.tif'): + """Get the name of the tile with lower left corner (lat, lon). + + Args: + lat: The minimum latitude of the tile. + lon: The minimum longitude of the tile. + is_worldcover: Wheter the tile is Worldcover or OSM. + postfix: A postfix to append to the tile name to make it a filename. + + Returns: + The name of the tile. + """ prefixes = ['N', 'S', 'E', 'W'] if is_worldcover else ['n', 's', 'e', 'w'] if lat >= 0: lat_part = prefixes[0] + str(int(lat)).zfill(2) @@ -16,4 +27,33 @@ def lat_lon_to_tile_string(lat, lon, is_worldcover: bool = False, postfix: str = lon_part = prefixes[2] + str(int(lon)).zfill(3) else: lon_part = prefixes[3] + str(int(np.abs(lon))).zfill(3) - return lat_part + lon_part + postfix \ No newline at end of file + return lat_part + lon_part + postfix + + +def merge_tiles(tiles, out_format, out_filename): + """Merge tiles via buildvrt and translate. + + Args: + tiles: The list of tiles to be merged. + out_format: The format of the output image. + out_filename: The name of the output COG. + """ + vrt = 'merged.vrt' + build_vrt_command = ' '.join(['gdalbuildvrt', vrt] + tiles) + translate_command = ' '.join(['gdal_translate', '-of', out_format, vrt, out_filename]) + os.system(build_vrt_command) + os.system(translate_command) + os.remove(vrt) + + +def remove_temp_files(temp_files: list): + """Remove each file in a list of files. + + Args: + temp_files: The list of temporary files to remove. + """ + for file in temp_files: + try: + os.remove(file) + except: + print('Error removing temporary file. Skipping it...') \ No newline at end of file From 4b5e8b5b97aa1b2ecfc1ddee025bf354d6168223 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Mon, 26 Feb 2024 09:00:01 -0600 Subject: [PATCH 07/58] added entrypoints for watermasking --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index cf9f2da4..4140cac9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,10 +41,12 @@ make_composite = "asf_tools.composite:main" water_map = "asf_tools.hydrosar.water_map:main" calculate_hand = "asf_tools.hydrosar.hand.calculate:main" flood_map = "asf_tools.hydrosar.flood_map:main" +water_masking = "asf_tools.watermasking.generate_osm_tiles:main" [project.entry-points.hyp3] water_map = "asf_tools.hydrosar.water_map:hyp3" flood_map = "asf_tools.hydrosar.flood_map:hyp3" +water_masking = "asf_tools.watermasking.generate_osm_tiles:main" [project.optional-dependencies] develop = [ From 6614f368e144d4cf62ae8d826b8cc1529734b1f2 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Mon, 26 Feb 2024 09:00:35 -0600 Subject: [PATCH 08/58] refactor and print statements --- .../watermasking/generate_osm_tiles.py | 31 +++++++++---------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index 1f60d50e..268568d0 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -8,6 +8,10 @@ from asf_tools.watermasking.utils import lat_lon_to_tile_string +INTERIOR_TILE_DIR = 'interior_tiles/' +OCEAN_TILE_DIR = 'ocean_tiles/' +FINISHED_TILE_DIR = 'tiles/' + def process_pbf(planet_file: str, output_file: str): """Process the global PBF file so that it only contains water features. @@ -170,7 +174,7 @@ def merge_tiles(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_ def setup_directories(): """Setup the directories necessary for running the script.""" - dirs = ['interior_tiles/', 'ocean_tiles/', 'merged_tiles/', 'finished_tiles/'] + dirs = [INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR] for dir in dirs: try: os.mkdir(dir) @@ -196,25 +200,17 @@ def main(): args = parser.parse_args() - # Setup Directories - interior_tile_dir = 'interior_tiles/' - ocean_tile_dir = 'ocean_tiles/' - finished_tile_dir = 'tiles/' - - dirs = [interior_tile_dir, ocean_tile_dir, finished_tile_dir] - for dir in dirs: - try: - os.mkdir(dir) - except FileExistsError as e: - print(f'{dir} already exists. Skipping...') + setup_directories() - # Process the PBF into one that contains only water. + print('Extracting water from planet file...') processed_pbf_path = 'planet_processed.pbf' process_pbf(args.planet_file_path, processed_pbf_path) - # Extract tiles from the processed and ocean files. + print('Processing tiles...') lat_range = range(args.lat_begin, args.lat_end, args.tile_height) lon_range = range(args.lon_begin, args.lon_end, args.tile_width) + num_tiles = len(lat_range) * len(lon_range) + index = 0 for lat in lat_range: for lon in lon_range: tile_name = lat_lon_to_tile_string(lat, lon, is_worldcover=False) @@ -224,12 +220,13 @@ def main(): process_ocean_tiles(args.ocean_polygons_path, lat, lon, args.tile_width, args.tile_height) end_time = time.time() total_time = end_time - start_time #seconds - print(f'Finished initial creation of {tile_name} in {total_time}(s).') + print(f'Finished initial creation of {tile_name} in {total_time}(s). {index} of {num_tiles}') + index += 1 except Exception as e: print(f'Caught error while processing {tile_name}. Continuing anyways...\n{e}') - # Merge the extracted processed/ocean tiles. - merge_tiles(interior_tile_dir, ocean_tile_dir) + print('Merging processed tiles...') + merge_tiles(INTERIOR_TILE_DIR, OCEAN_TILE_DIR) if __name__ == '__main__': From 9980a72b14d2accc161fa113b668449ddb962966 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Mon, 26 Feb 2024 10:35:17 -0600 Subject: [PATCH 09/58] added script for filling polar regions --- .../watermasking/fill_missing_tiles.py | 68 +++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 src/asf_tools/watermasking/fill_missing_tiles.py diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py new file mode 100644 index 00000000..ce46ee57 --- /dev/null +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -0,0 +1,68 @@ +import argparse +import os + +import numpy as np +from osgeo import gdal, osr + +from asf_tools.watermasking.utils import lat_lon_to_tile_string + + +def main(): + + parser = argparse.ArgumentParser( + prog='fill_missing_tiles.py', + description='Script for creating filled tifs in areas with missing tiles.' + ) + + parser.add_argument('--fill-value', help='The value to fill the data array with.', default=0) + parser.add_argument('--lat-begin', help='The minimum latitude of the dataset.', default=-85) + parser.add_argument('--lat-end', help='The maximum latitude of the dataset.', default=85) + parser.add_argument('--lon-begin', help='The minimum longitude of the dataset.', default=-180) + parser.add_argument('--lon-end', help='The maximum longitude of the dataset.', default=180) + parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) + parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) + + args = parser.parse_args() + + fill_value = int(args.fill_value) + lat_begin = int(args.lat_begin) + lat_end = int(args.lat_end) + lon_begin = int(args.lon_begin) + lon_end = int(args.lon_end) + tile_width = int(args.tile_width) + tile_height = int(args.tile_height) + + lat_range = range(lat_begin, lat_end, tile_width) + lon_range = range(lon_begin, lon_end, tile_width) + + for lat in lat_range: + for lon in lon_range: + + tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') + tile_tif = 'tiles/' + tile + '.tif' + + print(f'Processing: {tile}') + + xmin, xmax, ymin, ymax = lon, lon+tile_width, lat, lat+tile_height + pixel_size_x = 0.00009009009 # 10m * 2 at the equator. + pixel_size_y = 0.00009009009 + + # All images in the dataset should be this size. + data = np.empty((55500, 55500)) + data.fill(fill_value) + + driver = gdal.GetDriverByName('GTiff') + dst_ds = driver.Create(tile_tif, xsize=data.shape[0], ysize=data.shape[1], bands=1, eType=gdal.GDT_Byte) + dst_ds.SetGeoTransform( [ xmin, pixel_size_x, 0, ymax, 0, pixel_size_y ] ) + # wkt = f'POLYGON(({xmin} {ymin}, {xmax} {ymin}, {xmax} {ymax}, {xmin} {ymax}, {xmin} {ymin}))' + srs = osr.SpatialReference() + srs.ImportFromEPSG(4326) + dst_ds.SetProjection(srs) + dst_band = dst_ds.GetRasterBand(1) + dst_band.WriteArray(data) + dst_ds = None + del dst_ds + + +if __name__ == '__main__': + main() \ No newline at end of file From afe6f124ed88bf0ab54c7380ac883e108b77cdf2 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Mon, 26 Feb 2024 11:23:08 -0600 Subject: [PATCH 10/58] gdal expects lower left corner at these lats --- src/asf_tools/watermasking/fill_missing_tiles.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py index ce46ee57..050e9676 100644 --- a/src/asf_tools/watermasking/fill_missing_tiles.py +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -40,6 +40,7 @@ def main(): tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') tile_tif = 'tiles/' + tile + '.tif' + tile_cog = 'tiles/cogs/' + tile + '.tif' print(f'Processing: {tile}') @@ -53,16 +54,18 @@ def main(): driver = gdal.GetDriverByName('GTiff') dst_ds = driver.Create(tile_tif, xsize=data.shape[0], ysize=data.shape[1], bands=1, eType=gdal.GDT_Byte) - dst_ds.SetGeoTransform( [ xmin, pixel_size_x, 0, ymax, 0, pixel_size_y ] ) - # wkt = f'POLYGON(({xmin} {ymin}, {xmax} {ymin}, {xmax} {ymax}, {xmin} {ymax}, {xmin} {ymin}))' + dst_ds.SetGeoTransform( [ xmin, pixel_size_x, 0, ymin, 0, pixel_size_y ] ) srs = osr.SpatialReference() srs.ImportFromEPSG(4326) - dst_ds.SetProjection(srs) + dst_ds.SetProjection(srs.ExportToWkt()) dst_band = dst_ds.GetRasterBand(1) dst_band.WriteArray(data) dst_ds = None del dst_ds + command = f'gdal_translate -of COG -co NUM_THREADS=all_cpus {tile_tif} {tile_cog}' + os.system(command) + os.remove(tile_tif) if __name__ == '__main__': main() \ No newline at end of file From 023728ee189c222c93271db99a0e8fb3660465ad Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Mon, 26 Feb 2024 11:23:27 -0600 Subject: [PATCH 11/58] added entrypoint for filling missing tiles --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 4140cac9..a3aaac2a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,12 +41,12 @@ make_composite = "asf_tools.composite:main" water_map = "asf_tools.hydrosar.water_map:main" calculate_hand = "asf_tools.hydrosar.hand.calculate:main" flood_map = "asf_tools.hydrosar.flood_map:main" -water_masking = "asf_tools.watermasking.generate_osm_tiles:main" +generate_osm_tiles = "asf_tools.watermasking.generate_osm_tiles:main" +fill_missing_tiles = "asf_tools.watermasking.fill_missing_tiles:main" [project.entry-points.hyp3] water_map = "asf_tools.hydrosar.water_map:hyp3" flood_map = "asf_tools.hydrosar.flood_map:hyp3" -water_masking = "asf_tools.watermasking.generate_osm_tiles:main" [project.optional-dependencies] develop = [ From aeeff8da946dd10bc045d29e990f5a97627701cc Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Mon, 26 Feb 2024 11:55:01 -0600 Subject: [PATCH 12/58] refactoring --- .../watermasking/generate_worldcover_tiles.py | 26 ++++++++++++++----- 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 94bc8571..8ed26318 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -94,7 +94,7 @@ def get_tiles(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): return tiles -def lat_lon_to_filenames(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): +def lat_lon_to_filenames(worldcover_tile_dir, osm_tile_coord: tuple, wc_deg: int, osm_deg: int): """Get a list of the Worldcover tile filenames that are necessary to overlap an OSM tile. Args: @@ -108,12 +108,16 @@ def lat_lon_to_filenames(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): filenames = [] tiles = get_tiles(osm_tile_coord, wc_deg, osm_deg) for tile in tiles: - filenames.append(lat_lon_to_tile_string(tile[0], tile[1], is_worldcover=True)) + filenames.append(worldcover_tile_dir + lat_lon_to_tile_string(tile[0], tile[1], is_worldcover=True)) return filenames def crop_tile(tile): - """Crop the merged tiles""" + """Crop the merged tiles + + Args: + tile: The filename of the desired tile to crop. + """ try: ref_image = TILE_DIR + tile pixel_size = gdal.Warp('tmp_px_size.tif', ref_image, dstSRS='EPSG:4326').GetGeoTransform()[1] @@ -137,12 +141,20 @@ def crop_tile(tile): index += 1 -def build_dataset(lat_range, lon_range, worldcover_degrees, osm_degrees): +def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): + """ Main function for generating a dataset with worldcover tiles. + + Args: + worldcover_tile_dir: The directory containing the unprocessed worldcover tiles. + lat_range: The range of latitudes the dataset should cover. + lon_range: The range of longitudes the dataset should cover. + out_degrees: The width of the outputed dataset tiles in degrees. + """ for lat in lat_range: for lon in lon_range: start_time = time.time() tile_filename = TILE_DIR + lat_lon_to_tile_string(lat, lon, is_worldcover=False) - worldcover_tiles = lat_lon_to_filenames(lat, lon, worldcover_degrees, osm_degrees) + worldcover_tiles = lat_lon_to_filenames(worldcover_tile_dir, lat, lon, WORLDCOVER_TILE_SIZE, out_degrees) print(f'Processing: {tile_filename} {worldcover_tiles}') merge_tiles(worldcover_tiles, tile_filename) end_time = time.time() @@ -151,10 +163,12 @@ def build_dataset(lat_range, lon_range, worldcover_degrees, osm_degrees): def main(): + parser = argparse.ArgumentParser( prog='generate_worldcover_tiles.py', description='Main script for creating a tiled watermask dataset from the ESA WorldCover dataset.' ) + parser.add_argument('--worldcover-tiles-dir', help='The path to the directory containing the worldcover tifs.') parser.add_argument('--lat-begin', help='The minimum latitude of the dataset.', default=-85, required=True) parser.add_argument('--lat-end', help='The maximum latitude of the dataset.', default=85) @@ -176,7 +190,7 @@ def main(): lat_range = range(args.lat_begin, args.lat_end, args.tile_width) lon_range = range(args.lon_begin, args.lon_end, args.tile_heigth) - build_dataset(lat_range, lon_range, worldcover_degrees=WORLDCOVER_TILE_SIZE, osm_degrees=args.tile_width) + build_dataset(args.worldcover_tile_dir, lat_range, lon_range, worldcover_degrees=WORLDCOVER_TILE_SIZE, osm_degrees=args.tile_width) if __name__ == '__main__': From 64d5ed51428eb8b49ace46739ef8e81aa945c817 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Mon, 26 Feb 2024 11:55:25 -0600 Subject: [PATCH 13/58] refactoring --- .../watermasking/generate_osm_tiles.py | 42 ++++++++++++------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index 268568d0..b30802dc 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -49,7 +49,7 @@ def remove_temp_files(temp_files: list): print('Error removing temporary file. Skipping it...') -def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_height_deg): +def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_height_deg, output_dir): """Process and crop OSM ocean polygons into a tif tile. Args: @@ -61,7 +61,7 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig """ tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') - tile_tif = 'tiles/' + tile + '.tif' + tile_tif = output_dir + tile + '.tif' xmin, xmax, ymin, ymax = lon, lon+tile_width_deg, lat, lat+tile_height_deg pixel_size_x = 0.00009009009 # 10m * 2 at the equator. @@ -78,15 +78,14 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig yRes=pixel_size_y, burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], - outputType=gdal.GDT_Byte, - creationOptions=['COMPRESS=LZW'] + outputType=gdal.GDT_Byte ) temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx'] remove_temp_files(temp_files) -def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg): +def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interior_tile_dir): """Rasterize a water tile from the processed global PBF file. Args: @@ -99,7 +98,7 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg): tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') tile_pbf = tile+ '.osm.pbf' - tile_tif = tile + '.tif' + tile_tif = interior_tile_dir + tile + '.tif' tile_shp = tile + '.shp' tile_geojson = tile + '.geojson' @@ -130,10 +129,9 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg): burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], outputType=gdal.GDT_Byte, - creationOptions=['COMPRESS=LZW'] ) - temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx'] + temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx', tile_shp, tile_pbf, tile_geojson] remove_temp_files(temp_files) @@ -155,11 +153,18 @@ def merge_tiles(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_ internal_tile = internal_tile_dir + filename external_tile = ocean_tile_dir + filename output_tile = finished_tile_dir + filename - command = f'gdal_calc.py -A {internal_tile} -B {external_tile} --format GTiff -co NUM_THREADS=all_cpus -co COMPRESS=LZW -co TILED=YES -co PREDICTOR=YES --outfile {output_tile} --calc "logical_or(A, B)"' + command = f'gdal_calc.py -A {internal_tile} -B {external_tile} --format GTiff --outfile {output_tile} --calc "logical_or(A, B)"' os.system(command) if translate_to_cog: - command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus' + cogs_dir = finished_tile_dir + 'cogs/' + try: + os.mkdir(cogs_dir) + except FileExistsError as e: + pass + out_file = cogs_dir + filename + command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus {output_tile} {out_file}' + os.system(command) end_time = time.time() total_time = end_time - start_time @@ -200,6 +205,13 @@ def main(): args = parser.parse_args() + lat_begin = int(args.lat_begin) + lat_end = int(args.lat_end) + lon_begin = int(args.lon_begin) + lon_end = int(args.lon_end) + tile_width = int(args.tile_width) + tile_height = int(args.tile_height) + setup_directories() print('Extracting water from planet file...') @@ -207,8 +219,8 @@ def main(): process_pbf(args.planet_file_path, processed_pbf_path) print('Processing tiles...') - lat_range = range(args.lat_begin, args.lat_end, args.tile_height) - lon_range = range(args.lon_begin, args.lon_end, args.tile_width) + lat_range = range(lat_begin, lat_end, tile_height) + lon_range = range(lon_begin, lon_end, tile_width) num_tiles = len(lat_range) * len(lon_range) index = 0 for lat in lat_range: @@ -216,8 +228,8 @@ def main(): tile_name = lat_lon_to_tile_string(lat, lon, is_worldcover=False) try: start_time = time.time() - extract_water(processed_pbf_path, lat, lon, args.tile_width, args.tile_height) - process_ocean_tiles(args.ocean_polygons_path, lat, lon, args.tile_width, args.tile_height) + extract_water(processed_pbf_path, lat, lon, tile_width, tile_height, interior_tile_dir=INTERIOR_TILE_DIR) + process_ocean_tiles(args.ocean_polygons_path, lat, lon, tile_width, tile_height, output_dir=OCEAN_TILE_DIR) end_time = time.time() total_time = end_time - start_time #seconds print(f'Finished initial creation of {tile_name} in {total_time}(s). {index} of {num_tiles}') @@ -226,7 +238,7 @@ def main(): print(f'Caught error while processing {tile_name}. Continuing anyways...\n{e}') print('Merging processed tiles...') - merge_tiles(INTERIOR_TILE_DIR, OCEAN_TILE_DIR) + merge_tiles(INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR, translate_to_cog=True) if __name__ == '__main__': From ce58908b3047eeffd26449165121ba7c339b126d Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Tue, 27 Feb 2024 12:07:56 -0600 Subject: [PATCH 14/58] functionality to fill missing worldcover tiles --- .../watermasking/generate_worldcover_tiles.py | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 8ed26318..414928f2 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -141,6 +141,42 @@ def crop_tile(tile): index += 1 +def create_missing_tiles(tile_dir, lat_range, lon_range): + """Check for, and build, tiles that may be missing from WorldCover, such as over the ocean. + + Args: + lat_range: The range of latitudes to check. + lon_range: The range of longitudes to check. + Returns: + current_existing_tiles: The list of tiles that exist after the function has completed. + """ + current_existing_tiles = [f for f in os.listdir(tile_dir) if f.endswith(FILENAME_POSTFIX)] + for lon in lon_range: + for lat in lat_range: + tile = lat_lon_to_tile_string(lat, lon, is_worldcover=True) + print(f'Checking {tile}') + if tile not in current_existing_tiles: + print(f'Could not find {tile}') + x_size = 36000 + y_size = 36000 + x_res = 8.333333333333333055e-05 + y_res = -8.333333333333333055e-05 + ul_lon = lon + ul_lat = lat + WORLDCOVER_TILE_SIZE + geotransform = (ul_lon, x_res, 0, ul_lat, 0, y_res) + driver = gdal.GetDriverByName('GTiff') + ds = driver.Create(tile, xsize=x_size, ysize=y_size, bands=1, eType=gdal.GDT_Byte, options=['COMPRESS=LZW']) + ds.SetProjection('EPSG:4326') + ds.SetGeoTransform(geotransform) + band = ds.GetRasterBand(1) # Write ones, as tiles should only be missing over water. + band.WriteArray(np.ones((x_size, y_size))) + del ds + del band + current_existing_tiles.append(tile) + print(f'Added {tile}') + return current_existing_tiles + + def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): """ Main function for generating a dataset with worldcover tiles. @@ -157,6 +193,7 @@ def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): worldcover_tiles = lat_lon_to_filenames(worldcover_tile_dir, lat, lon, WORLDCOVER_TILE_SIZE, out_degrees) print(f'Processing: {tile_filename} {worldcover_tiles}') merge_tiles(worldcover_tiles, tile_filename) + create_missing_tiles() end_time = time.time() total_time = end_time - start_time print(f'Time Elapsed: {total_time}s') From b7ed7d3ef8412bee9c30b7a783fa95a9eff08b32 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 27 Feb 2024 19:32:10 +0000 Subject: [PATCH 15/58] Bump pypa/gh-action-pypi-publish from 1.8.11 to 1.8.12 Bumps [pypa/gh-action-pypi-publish](https://github.com/pypa/gh-action-pypi-publish) from 1.8.11 to 1.8.12. - [Release notes](https://github.com/pypa/gh-action-pypi-publish/releases) - [Commits](https://github.com/pypa/gh-action-pypi-publish/compare/v1.8.11...v1.8.12) --- updated-dependencies: - dependency-name: pypa/gh-action-pypi-publish dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] --- .github/workflows/distribute.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/distribute.yml b/.github/workflows/distribute.yml index bad4d8b1..1074e21e 100644 --- a/.github/workflows/distribute.yml +++ b/.github/workflows/distribute.yml @@ -29,7 +29,7 @@ jobs: python -m build - name: upload to PyPI.org - uses: pypa/gh-action-pypi-publish@v1.8.11 + uses: pypa/gh-action-pypi-publish@v1.8.12 with: user: __token__ password: ${{ secrets.TOOLS_PYPI_PAK }} From 936721815a966823c2ad6e4abe44fc16fe194381 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 12:06:55 -0600 Subject: [PATCH 16/58] refactor --- .../watermasking/generate_osm_tiles.py | 42 +++++-------------- 1 file changed, 11 insertions(+), 31 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index b30802dc..89779c7c 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -6,7 +6,7 @@ import numpy as np from osgeo import gdal -from asf_tools.watermasking.utils import lat_lon_to_tile_string +from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, setup_directories INTERIOR_TILE_DIR = 'interior_tiles/' OCEAN_TILE_DIR = 'ocean_tiles/' @@ -36,19 +36,6 @@ def process_pbf(planet_file: str, output_file: str): os.system(merge_command) -def remove_temp_files(temp_files: list): - """Remove each file in a list of files. - - Args: - temp_files: The list of temporary files to remove. - """ - for file in temp_files: - try: - os.remove(file) - except: - print('Error removing temporary file. Skipping it...') - - def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_height_deg, output_dir): """Process and crop OSM ocean polygons into a tif tile. @@ -78,7 +65,8 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig yRes=pixel_size_y, burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], - outputType=gdal.GDT_Byte + outputType=gdal.GDT_Byte, + creationOptions=['COMPRESS=LZW'] ) temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx'] @@ -129,15 +117,16 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], outputType=gdal.GDT_Byte, + creationOptions=['COMPRESS=LZW'] ) temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx', tile_shp, tile_pbf, tile_geojson] remove_temp_files(temp_files) -def merge_tiles(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_to_cog: bool = False): +def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_to_cog: bool = False): """Merge the interior water tiles and ocean water tiles. - + Args: interior_tile_dir: The path to the directory containing the interior water tiles. ocean_tile_dir: The path to the directory containing the ocean water tiles. @@ -165,6 +154,7 @@ def merge_tiles(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_ out_file = cogs_dir + filename command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus {output_tile} {out_file}' os.system(command) + os.remove(output_tile) end_time = time.time() total_time = end_time - start_time @@ -177,16 +167,6 @@ def merge_tiles(internal_tile_dir, ocean_tile_dir, finished_tile_dir, translate_ print(f'Caught error while processing {filename}. Continuing anyways...\n{e}') -def setup_directories(): - """Setup the directories necessary for running the script.""" - dirs = [INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR] - for dir in dirs: - try: - os.mkdir(dir) - except FileExistsError as e: - print(f'{dir} already exists. Skipping...') - - def main(): parser = argparse.ArgumentParser( @@ -212,7 +192,7 @@ def main(): tile_width = int(args.tile_width) tile_height = int(args.tile_height) - setup_directories() + setup_directories([INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR]) print('Extracting water from planet file...') processed_pbf_path = 'planet_processed.pbf' @@ -237,9 +217,9 @@ def main(): except Exception as e: print(f'Caught error while processing {tile_name}. Continuing anyways...\n{e}') - print('Merging processed tiles...') - merge_tiles(INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR, translate_to_cog=True) + print('Merging processed tiles...') + merge_interior_and_ocean(INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR, translate_to_cog=True) if __name__ == '__main__': - main() \ No newline at end of file + main() From 2314c808f71a14559bf084b8fb8b4e5e155460f9 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 12:07:04 -0600 Subject: [PATCH 17/58] refactor --- src/asf_tools/watermasking/fill_missing_tiles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py index 050e9676..1e3c2e9f 100644 --- a/src/asf_tools/watermasking/fill_missing_tiles.py +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -68,4 +68,4 @@ def main(): os.remove(tile_tif) if __name__ == '__main__': - main() \ No newline at end of file + main() From 5308f758a4be6529352d0c53467bd6e5761149f6 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 12:07:23 -0600 Subject: [PATCH 18/58] refactor --- .../watermasking/generate_worldcover_tiles.py | 107 ++++++++++-------- 1 file changed, 57 insertions(+), 50 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 414928f2..b422b646 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -6,9 +6,10 @@ from osgeo import gdal -from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, merge_tiles +from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, merge_tiles, setup_directories -TILE_DIR = 'worldcover_tiles_unfinished/' +PROCESSED_TILE_DIR = 'worldcover_tiles_preprocessed/' +TILE_DIR = 'worldcover_tiles_uncropped/' CROPPED_TILE_DIR = 'worldcover_tiles/' FILENAME_POSTFIX = '.tif' WORLDCOVER_TILE_SIZE = 3 @@ -29,9 +30,11 @@ def tile_preprocessing(tile_dir, min_lat, max_lat): index = 0 num_tiles = len(filenames_filtered) for filename in filenames_filtered: + start_time = time.time() - dst_filename = filename.split('_')[5] + '.tif' + filename = tile_dir + filename + dst_filename = PROCESSED_TILE_DIR + filename.split('_')[5] + '.tif' print(f'Processing: {filename} --- {dst_filename} -- {index} of {num_tiles}') @@ -56,13 +59,51 @@ def tile_preprocessing(tile_dir, min_lat, max_lat): del src_ds end_time = time.time() - total_time = end_time - start_time # seconds + total_time = end_time - start_time print(f'Processing {dst_filename} took {total_time} seconds.') index += 1 +def create_missing_tiles(tile_dir, lat_range, lon_range): + """Check for, and build, tiles that may be missing from WorldCover, such as over the ocean. + + Args: + lat_range: The range of latitudes to check. + lon_range: The range of longitudes to check. + Returns: + current_existing_tiles: The list of tiles that exist after the function has completed. + """ + current_existing_tiles = [f for f in os.listdir(tile_dir) if f.endswith(FILENAME_POSTFIX)] + for lon in lon_range: + for lat in lat_range: + tile = lat_lon_to_tile_string(lat, lon, is_worldcover=True) + print(f'Checking {tile}') + if tile not in current_existing_tiles: + print(f'Could not find {tile}') + + x_size, y_size = 36000, 36000 + x_res, y_res = 8.333333333333333055e-05, -8.333333333333333055e-05 + ul_lon = lon + ul_lat = lat + WORLDCOVER_TILE_SIZE + geotransform = (ul_lon, x_res, 0, ul_lat, 0, y_res) + + driver = gdal.GetDriverByName('GTiff') + ds = driver.Create(tile, xsize=x_size, ysize=y_size, bands=1, eType=gdal.GDT_Byte, options=['COMPRESS=LZW']) + ds.SetProjection('EPSG:4326') + ds.SetGeoTransform(geotransform) + band = ds.GetRasterBand(1) # Write ones, as tiles should only be missing over water. + band.WriteArray(np.ones((x_size, y_size))) + + del ds + del band + + current_existing_tiles.append(tile) + print(f'Added {tile}') + return current_existing_tiles + + def get_tiles(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): """Get a list of the worldcover tile locations necessary to fully cover an OSM tile. @@ -114,16 +155,18 @@ def lat_lon_to_filenames(worldcover_tile_dir, osm_tile_coord: tuple, wc_deg: int def crop_tile(tile): """Crop the merged tiles - + Args: tile: The filename of the desired tile to crop. """ try: ref_image = TILE_DIR + tile + out_filename = CROPPED_TILE_DIR + tile pixel_size = gdal.Warp('tmp_px_size.tif', ref_image, dstSRS='EPSG:4326').GetGeoTransform()[1] + shapefile_command = ' '.join(['gdaltindex', 'tmp.shp', ref_image]) os.system(shapefile_command) - out_filename = CROPPED_TILE_DIR + tile + gdal.Warp( out_filename, tile, @@ -137,46 +180,10 @@ def crop_tile(tile): ) remove_temp_files(['tmp_px_size.tif', 'tmp.shp']) except Exception as e: # When a tile fails to process, we don't necessarily want the program to stop. - print(f'Could not process {tile}. Continuing...') + print(f'Could not process {tile} due to {e}. Skipping...') index += 1 -def create_missing_tiles(tile_dir, lat_range, lon_range): - """Check for, and build, tiles that may be missing from WorldCover, such as over the ocean. - - Args: - lat_range: The range of latitudes to check. - lon_range: The range of longitudes to check. - Returns: - current_existing_tiles: The list of tiles that exist after the function has completed. - """ - current_existing_tiles = [f for f in os.listdir(tile_dir) if f.endswith(FILENAME_POSTFIX)] - for lon in lon_range: - for lat in lat_range: - tile = lat_lon_to_tile_string(lat, lon, is_worldcover=True) - print(f'Checking {tile}') - if tile not in current_existing_tiles: - print(f'Could not find {tile}') - x_size = 36000 - y_size = 36000 - x_res = 8.333333333333333055e-05 - y_res = -8.333333333333333055e-05 - ul_lon = lon - ul_lat = lat + WORLDCOVER_TILE_SIZE - geotransform = (ul_lon, x_res, 0, ul_lat, 0, y_res) - driver = gdal.GetDriverByName('GTiff') - ds = driver.Create(tile, xsize=x_size, ysize=y_size, bands=1, eType=gdal.GDT_Byte, options=['COMPRESS=LZW']) - ds.SetProjection('EPSG:4326') - ds.SetGeoTransform(geotransform) - band = ds.GetRasterBand(1) # Write ones, as tiles should only be missing over water. - band.WriteArray(np.ones((x_size, y_size))) - del ds - del band - current_existing_tiles.append(tile) - print(f'Added {tile}') - return current_existing_tiles - - def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): """ Main function for generating a dataset with worldcover tiles. @@ -193,7 +200,7 @@ def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): worldcover_tiles = lat_lon_to_filenames(worldcover_tile_dir, lat, lon, WORLDCOVER_TILE_SIZE, out_degrees) print(f'Processing: {tile_filename} {worldcover_tiles}') merge_tiles(worldcover_tiles, tile_filename) - create_missing_tiles() + crop_tile(tile_filename) end_time = time.time() total_time = end_time - start_time print(f'Time Elapsed: {total_time}s') @@ -216,14 +223,14 @@ def main(): args = parser.parse_args() - for dir in [TILE_DIR, CROPPED_TILE_DIR]: - try: - os.mkdir(dir) - except FileExistsError as e: - print(f'{dir} already exists. Skipping...') + setup_directories([PROCESSED_TILE_DIR, TILE_DIR, CROPPED_TILE_DIR]) + # Process the multi-class masks into water/not-water masks. tile_preprocessing(args.worldcover_tiles_dir, args.lat_begin, args.lat_end) + # Ocean only tiles are missing from WorldCover, so we need to create blank (water-only) ones. + create_missing_tiles(PROCESSED_TILE_DIR, lat_range, lon_range) + lat_range = range(args.lat_begin, args.lat_end, args.tile_width) lon_range = range(args.lon_begin, args.lon_end, args.tile_heigth) @@ -231,4 +238,4 @@ def main(): if __name__ == '__main__': - main() \ No newline at end of file + main() From 4543ee95c6141f56a39580e9b3b2262ab3b2e283 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 12:07:37 -0600 Subject: [PATCH 19/58] added setup_directories to utils --- src/asf_tools/watermasking/utils.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py index 15b56bd4..511323da 100644 --- a/src/asf_tools/watermasking/utils.py +++ b/src/asf_tools/watermasking/utils.py @@ -55,5 +55,15 @@ def remove_temp_files(temp_files: list): for file in temp_files: try: os.remove(file) - except: - print('Error removing temporary file. Skipping it...') \ No newline at end of file + except Exception as e: + print(f'Caught {e} while removing temporary file: {file}. Skipping it...') + + +def setup_directories(dirs: list[str]): + """Setup the directories necessary for running the script.""" + for dir in dirs: + try: + os.mkdir(dir) + except FileExistsError as e: + # Directories already exists. + pass \ No newline at end of file From 796496303785319db26d235585f0a95796effa3e Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 12:10:48 -0600 Subject: [PATCH 20/58] refactor --- src/asf_tools/watermasking/utils.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py index 511323da..0cce7a86 100644 --- a/src/asf_tools/watermasking/utils.py +++ b/src/asf_tools/watermasking/utils.py @@ -1,10 +1,7 @@ import os -import time import numpy as np -from osgeo import gdal - def lat_lon_to_tile_string(lat, lon, is_worldcover: bool = False, postfix: str ='.tif'): """Get the name of the tile with lower left corner (lat, lon). @@ -66,4 +63,4 @@ def setup_directories(dirs: list[str]): os.mkdir(dir) except FileExistsError as e: # Directories already exists. - pass \ No newline at end of file + pass From 29ea0cd233777de57b0ea4f22d44e4b19c094a0c Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 12:12:13 -0600 Subject: [PATCH 21/58] refactoring --- src/asf_tools/watermasking/generate_osm_tiles.py | 1 - src/asf_tools/watermasking/generate_worldcover_tiles.py | 3 +-- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index 89779c7c..4bb4272e 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -3,7 +3,6 @@ import time import geopandas as gpd -import numpy as np from osgeo import gdal from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, setup_directories diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index b422b646..e5c1c633 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -3,10 +3,9 @@ import time import numpy as np - from osgeo import gdal -from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, merge_tiles, setup_directories +from asf_tools.watermasking.utils import lat_lon_to_tile_string, merge_tiles, remove_temp_files, setup_directories PROCESSED_TILE_DIR = 'worldcover_tiles_preprocessed/' TILE_DIR = 'worldcover_tiles_uncropped/' From e0b8bfdc625d9445dc47ba03d9e513c9e2813a8f Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 13:11:23 -0600 Subject: [PATCH 22/58] making flake8 happy --- .../watermasking/fill_missing_tiles.py | 9 +-- .../watermasking/generate_osm_tiles.py | 61 +++++++++++++------ 2 files changed, 49 insertions(+), 21 deletions(-) diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py index 1e3c2e9f..0ba7c540 100644 --- a/src/asf_tools/watermasking/fill_missing_tiles.py +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -32,7 +32,7 @@ def main(): tile_width = int(args.tile_width) tile_height = int(args.tile_height) - lat_range = range(lat_begin, lat_end, tile_width) + lat_range = range(lat_begin, lat_end, tile_height) lon_range = range(lon_begin, lon_end, tile_width) for lat in lat_range: @@ -44,8 +44,8 @@ def main(): print(f'Processing: {tile}') - xmin, xmax, ymin, ymax = lon, lon+tile_width, lat, lat+tile_height - pixel_size_x = 0.00009009009 # 10m * 2 at the equator. + xmin, ymin = lon, lat + pixel_size_x = 0.00009009009 pixel_size_y = 0.00009009009 # All images in the dataset should be this size. @@ -54,7 +54,7 @@ def main(): driver = gdal.GetDriverByName('GTiff') dst_ds = driver.Create(tile_tif, xsize=data.shape[0], ysize=data.shape[1], bands=1, eType=gdal.GDT_Byte) - dst_ds.SetGeoTransform( [ xmin, pixel_size_x, 0, ymin, 0, pixel_size_y ] ) + dst_ds.SetGeoTransform([xmin, pixel_size_x, 0, ymin, 0, pixel_size_y]) srs = osr.SpatialReference() srs.ImportFromEPSG(4326) dst_ds.SetProjection(srs.ExportToWkt()) @@ -67,5 +67,6 @@ def main(): os.system(command) os.remove(tile_tif) + if __name__ == '__main__': main() diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index 4bb4272e..ef1c0f3b 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -17,7 +17,7 @@ def process_pbf(planet_file: str, output_file: str): Args: planet_file: The path to the OSM Planet PBF file. - output_file: The desired path of the processed PBF file. + output_file: The desired path of the processed PBF file. """ natural_file = 'planet_natural.pbf' @@ -37,7 +37,7 @@ def process_pbf(planet_file: str, output_file: str): def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_height_deg, output_dir): """Process and crop OSM ocean polygons into a tif tile. - + Args: ocean_polygons_path: The path to the global ocean polygons file from OSM. lat: The minimum latitude of the tile. @@ -50,7 +50,7 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig tile_tif = output_dir + tile + '.tif' xmin, xmax, ymin, ymax = lon, lon+tile_width_deg, lat, lat+tile_height_deg - pixel_size_x = 0.00009009009 # 10m * 2 at the equator. + pixel_size_x = 0.00009009009 pixel_size_y = 0.00009009009 clipped_polygons_path = tile + '.shp' @@ -61,9 +61,9 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig tile_tif, clipped_polygons_path, xRes=pixel_size_x, - yRes=pixel_size_y, + yRes=pixel_size_y, burnValues=1, - outputBounds=[xmin, ymin, xmax, ymax], + outputBounds=[xmin, ymin, xmax, ymax], outputType=gdal.GDT_Byte, creationOptions=['COMPRESS=LZW'] ) @@ -84,13 +84,14 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio """ tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') - tile_pbf = tile+ '.osm.pbf' + tile_pbf = tile + '.osm.pbf' tile_tif = interior_tile_dir + tile + '.tif' tile_shp = tile + '.shp' tile_geojson = tile + '.geojson' # Extract tile from the main pbf, then convert it to a tif. - os.system(f'osmium extract -s smart -S tags=natural=water --bbox {lon},{lat},{lon+tile_width_deg},{lat+tile_height_deg} {water_file} -o {tile_pbf}') + bbox = f'--bbox {lon},{lat},{lon+tile_width_deg},{lat+tile_height_deg}' + os.system(f'osmium extract -s smart -S tags=natural=water {bbox} {water_file} -o {tile_pbf}') os.system(f'osmium export --geometry-types="polygon" {tile_pbf} -o {tile_geojson}') # Islands and Islets can be members of the water features, so they must be removed. @@ -105,16 +106,16 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio water_gdf = None xmin, xmax, ymin, ymax = lon, lon+tile_width_deg, lat, lat+tile_height_deg - pixel_size_x = 0.00009009009 # 10m at the equator. + pixel_size_x = 0.00009009009 pixel_size_y = 0.00009009009 gdal.Rasterize( tile_tif, tile_shp, - xRes=pixel_size_x, - yRes=pixel_size_y, + xRes=pixel_size_x, + yRes=pixel_size_y, burnValues=1, - outputBounds=[xmin, ymin, xmax, ymax], + outputBounds=[xmin, ymin, xmax, ymax], outputType=gdal.GDT_Byte, creationOptions=['COMPRESS=LZW'] ) @@ -141,14 +142,26 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di internal_tile = internal_tile_dir + filename external_tile = ocean_tile_dir + filename output_tile = finished_tile_dir + filename - command = f'gdal_calc.py -A {internal_tile} -B {external_tile} --format GTiff --outfile {output_tile} --calc "logical_or(A, B)"' + command = ' '.join([ + 'gdal_calc.py', + '-A', + internal_tile, + '-B', + external_tile, + '--format', + 'GTiff', + '--outfile', + output_tile, + '--calc', + '"logical_or(A, B)"' + ]) os.system(command) if translate_to_cog: cogs_dir = finished_tile_dir + 'cogs/' try: os.mkdir(cogs_dir) - except FileExistsError as e: + except FileExistsError: pass out_file = cogs_dir + filename command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus {output_tile} {out_file}' @@ -173,8 +186,8 @@ def main(): description='Main script for creating a tiled watermask dataset from OSM data.' ) - parser.add_argument('--planet-file-path', help='The path to the global planet.pbf file.', default='planet.pbf') - parser.add_argument('--ocean-polygons-path', help='The path to the global OSM ocean polygons.', default='water_polygons.shp') + parser.add_argument('--planet-file-path', help='The path to the global planet.pbf file.') + parser.add_argument('--ocean-polygons-path', help='The path to the global OSM ocean polygons.') parser.add_argument('--lat-begin', help='The minimum latitude of the dataset.', default=-85) parser.add_argument('--lat-end', help='The maximum latitude of the dataset.', default=85) parser.add_argument('--lon-begin', help='The minimum longitude of the dataset.', default=-180) @@ -207,8 +220,22 @@ def main(): tile_name = lat_lon_to_tile_string(lat, lon, is_worldcover=False) try: start_time = time.time() - extract_water(processed_pbf_path, lat, lon, tile_width, tile_height, interior_tile_dir=INTERIOR_TILE_DIR) - process_ocean_tiles(args.ocean_polygons_path, lat, lon, tile_width, tile_height, output_dir=OCEAN_TILE_DIR) + extract_water( + processed_pbf_path, + lat, + lon, + tile_width, + tile_height, + interior_tile_dir=INTERIOR_TILE_DIR + ) + process_ocean_tiles( + args.ocean_polygons_path, + lat, + lon, + tile_width, + tile_height, + output_dir=OCEAN_TILE_DIR + ) end_time = time.time() total_time = end_time - start_time #seconds print(f'Finished initial creation of {tile_name} in {total_time}(s). {index} of {num_tiles}') From 22d646bf2025f0a15a4542dc849d1b66500e663f Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 13:21:48 -0600 Subject: [PATCH 23/58] making flake8 happy --- .../watermasking/generate_osm_tiles.py | 4 +- .../watermasking/generate_worldcover_tiles.py | 50 +++++++++++++------ 2 files changed, 37 insertions(+), 17 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index ef1c0f3b..e2b90b53 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -99,7 +99,7 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio try: water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'island'].index) water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'islet'].index) - except: + except Exception as e: # When there are no islands to remove, an error will be occur, but we don't care about it. pass water_gdf.to_file(tile_shp, mode='w', engine='pyogrio') @@ -237,7 +237,7 @@ def main(): output_dir=OCEAN_TILE_DIR ) end_time = time.time() - total_time = end_time - start_time #seconds + total_time = end_time - start_time print(f'Finished initial creation of {tile_name} in {total_time}(s). {index} of {num_tiles}') index += 1 except Exception as e: diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index e5c1c633..690e021c 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -8,7 +8,7 @@ from asf_tools.watermasking.utils import lat_lon_to_tile_string, merge_tiles, remove_temp_files, setup_directories PROCESSED_TILE_DIR = 'worldcover_tiles_preprocessed/' -TILE_DIR = 'worldcover_tiles_uncropped/' +TILE_DIR = 'worldcover_tiles_uncropped/' CROPPED_TILE_DIR = 'worldcover_tiles/' FILENAME_POSTFIX = '.tif' WORLDCOVER_TILE_SIZE = 3 @@ -23,8 +23,9 @@ def tile_preprocessing(tile_dir, min_lat, max_lat): """ filenames = [f for f in os.listdir(tile_dir) if f.endswith('.tif')] - filter = lambda filename: (int(filename.split('_')[5][1:3]) >= min_lat) and (int(filename.split('_')[5][1:3]) <= max_lat) - filenames_filtered = [f for f in filenames if filter(f)] + get_lat = lambda filename: int(filename.split('_')[5][1:3]) + filter = lambda filename: (get_lat(filename) >= min_lat) and (get_lat(filename) <= max_lat) + filenames_filtered = [f for f in filenames if filter(f)] index = 0 num_tiles = len(filenames_filtered) @@ -47,7 +48,14 @@ def tile_preprocessing(tile_dir, min_lat, max_lat): driver = gdal.GetDriverByName('GTiff') - dst_ds = driver.Create(dst_filename, water_arr.shape[0], water_arr.shape[1], 1, gdal.GDT_Byte, options=['COMPRESS=LZW', 'TILED=YES']) + dst_ds = driver.Create( + dst_filename, + water_arr.shape[0], + water_arr.shape[1], + 1, + gdal.GDT_Byte, + options=['COMPRESS=LZW', 'TILED=YES'] + ) dst_ds.SetGeoTransform(src_ds.GetGeoTransform()) dst_ds.SetProjection(src_ds.GetProjection()) dst_band = dst_ds.GetRasterBand(1) @@ -89,7 +97,14 @@ def create_missing_tiles(tile_dir, lat_range, lon_range): geotransform = (ul_lon, x_res, 0, ul_lat, 0, y_res) driver = gdal.GetDriverByName('GTiff') - ds = driver.Create(tile, xsize=x_size, ysize=y_size, bands=1, eType=gdal.GDT_Byte, options=['COMPRESS=LZW']) + ds = driver.Create( + tile, + xsize=x_size, + ysize=y_size, + bands=1, + eType=gdal.GDT_Byte, + options=['COMPRESS=LZW'] + ) ds.SetProjection('EPSG:4326') ds.SetGeoTransform(geotransform) band = ds.GetRasterBand(1) # Write ones, as tiles should only be missing over water. @@ -105,12 +120,12 @@ def create_missing_tiles(tile_dir, lat_range, lon_range): def get_tiles(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): """Get a list of the worldcover tile locations necessary to fully cover an OSM tile. - + Args: osm_tile_coord: The lower left corner coordinate (lat, lon) of the desired OSM tile. wc_deg: The width/height of the Worldcover tiles in degrees. osm_deg: The width/height of the OSM tiles in degrees. - + Returns: tiles: A list of the lower left corner coordinates of the Worldcover tiles that overlap the OSM tile. """ @@ -141,7 +156,7 @@ def lat_lon_to_filenames(worldcover_tile_dir, osm_tile_coord: tuple, wc_deg: int osm_tile: The lower left corner (lat, lon) of the desired OSM tile. wc_deg: The width of the Worldcover tiles in degrees. osm_deg: The width of the OSM tiles in degrees. - + Returns: filenames: The list of Worldcover filenames. """ @@ -180,7 +195,6 @@ def crop_tile(tile): remove_temp_files(['tmp_px_size.tif', 'tmp.shp']) except Exception as e: # When a tile fails to process, we don't necessarily want the program to stop. print(f'Could not process {tile} due to {e}. Skipping...') - index += 1 def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): @@ -190,14 +204,14 @@ def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): worldcover_tile_dir: The directory containing the unprocessed worldcover tiles. lat_range: The range of latitudes the dataset should cover. lon_range: The range of longitudes the dataset should cover. - out_degrees: The width of the outputed dataset tiles in degrees. + out_degrees: The width of the outputed dataset tiles in degrees. """ for lat in lat_range: for lon in lon_range: start_time = time.time() tile_filename = TILE_DIR + lat_lon_to_tile_string(lat, lon, is_worldcover=False) worldcover_tiles = lat_lon_to_filenames(worldcover_tile_dir, lat, lon, WORLDCOVER_TILE_SIZE, out_degrees) - print(f'Processing: {tile_filename} {worldcover_tiles}') + print(f'Processing: {tile_filename} {worldcover_tiles}') merge_tiles(worldcover_tiles, tile_filename) crop_tile(tile_filename) end_time = time.time() @@ -227,13 +241,19 @@ def main(): # Process the multi-class masks into water/not-water masks. tile_preprocessing(args.worldcover_tiles_dir, args.lat_begin, args.lat_end) - # Ocean only tiles are missing from WorldCover, so we need to create blank (water-only) ones. - create_missing_tiles(PROCESSED_TILE_DIR, lat_range, lon_range) - lat_range = range(args.lat_begin, args.lat_end, args.tile_width) lon_range = range(args.lon_begin, args.lon_end, args.tile_heigth) - build_dataset(args.worldcover_tile_dir, lat_range, lon_range, worldcover_degrees=WORLDCOVER_TILE_SIZE, osm_degrees=args.tile_width) + # Ocean only tiles are missing from WorldCover, so we need to create blank (water-only) ones. + create_missing_tiles(PROCESSED_TILE_DIR, lat_range, lon_range) + + build_dataset( + args.worldcover_tile_dir, + lat_range, + lon_range, + worldcover_degrees=WORLDCOVER_TILE_SIZE, + osm_degrees=args.tile_width + ) if __name__ == '__main__': From 2ed5362741cf436230a6f3215878aa58fb81fe32 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 13:25:54 -0600 Subject: [PATCH 24/58] making flake8 happy --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 690e021c..94d30546 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -23,8 +23,8 @@ def tile_preprocessing(tile_dir, min_lat, max_lat): """ filenames = [f for f in os.listdir(tile_dir) if f.endswith('.tif')] - get_lat = lambda filename: int(filename.split('_')[5][1:3]) - filter = lambda filename: (get_lat(filename) >= min_lat) and (get_lat(filename) <= max_lat) + def get_lat(filename): int(filename.split('_')[5][1:3]) + def filter(filename): (get_lat(filename) >= min_lat) and (get_lat(filename) <= max_lat) filenames_filtered = [f for f in filenames if filter(f)] index = 0 @@ -101,7 +101,7 @@ def create_missing_tiles(tile_dir, lat_range, lon_range): tile, xsize=x_size, ysize=y_size, - bands=1, + bands=1, eType=gdal.GDT_Byte, options=['COMPRESS=LZW'] ) From 45e55c9006e48655a13b9c091c645d1c96da1d9c Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 13:29:52 -0600 Subject: [PATCH 25/58] restructured lambda --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 94d30546..8dd4b4fa 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -23,9 +23,10 @@ def tile_preprocessing(tile_dir, min_lat, max_lat): """ filenames = [f for f in os.listdir(tile_dir) if f.endswith('.tif')] - def get_lat(filename): int(filename.split('_')[5][1:3]) - def filter(filename): (get_lat(filename) >= min_lat) and (get_lat(filename) <= max_lat) - filenames_filtered = [f for f in filenames if filter(f)] + def filename_filter(filename): + latitude = int(filename.split('_')[5][1:3]) + (latitude >= min_lat) and (latitude <= max_lat) + filenames_filtered = [f for f in filenames if filename_filter(f)] index = 0 num_tiles = len(filenames_filtered) From be30e4858b961520ea08326f58f58287fc3f0d4c Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 14:02:22 -0600 Subject: [PATCH 26/58] remove blind exceptions --- .../watermasking/generate_osm_tiles.py | 120 +++++++++--------- .../watermasking/generate_worldcover_tiles.py | 41 +++--- src/asf_tools/watermasking/utils.py | 14 +- 3 files changed, 83 insertions(+), 92 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index e2b90b53..7bdcdaba 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -99,8 +99,8 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio try: water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'island'].index) water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'islet'].index) - except Exception as e: - # When there are no islands to remove, an error will be occur, but we don't care about it. + except AttributeError: + # When there are no islands to remove, an AttributeError should throw, but we don't care about it. pass water_gdf.to_file(tile_shp, mode='w', engine='pyogrio') water_gdf = None @@ -136,47 +136,44 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di num_tiles = len([f for f in os.listdir(internal_tile_dir) if f.endswith('tif')]) for filename in os.listdir(internal_tile_dir): if filename.endswith('.tif'): - try: - start_time = time.time() - - internal_tile = internal_tile_dir + filename - external_tile = ocean_tile_dir + filename - output_tile = finished_tile_dir + filename - command = ' '.join([ - 'gdal_calc.py', - '-A', - internal_tile, - '-B', - external_tile, - '--format', - 'GTiff', - '--outfile', - output_tile, - '--calc', - '"logical_or(A, B)"' - ]) + start_time = time.time() + + internal_tile = internal_tile_dir + filename + external_tile = ocean_tile_dir + filename + output_tile = finished_tile_dir + filename + command = ' '.join([ + 'gdal_calc.py', + '-A', + internal_tile, + '-B', + external_tile, + '--format', + 'GTiff', + '--outfile', + output_tile, + '--calc', + '"logical_or(A, B)"' + ]) + os.system(command) + + if translate_to_cog: + cogs_dir = finished_tile_dir + 'cogs/' + try: + os.mkdir(cogs_dir) + except FileExistsError: + pass + out_file = cogs_dir + filename + command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus {output_tile} {out_file}' os.system(command) + os.remove(output_tile) - if translate_to_cog: - cogs_dir = finished_tile_dir + 'cogs/' - try: - os.mkdir(cogs_dir) - except FileExistsError: - pass - out_file = cogs_dir + filename - command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus {output_tile} {out_file}' - os.system(command) - os.remove(output_tile) + end_time = time.time() + total_time = end_time - start_time - end_time = time.time() - total_time = end_time - start_time + print(f'Elapsed Time: {total_time}(s)') + print(f'Completed {index} of {num_tiles}') - print(f'Elapsed Time: {total_time}(s)') - print(f'Completed {index} of {num_tiles}') - - index += 1 - except Exception as e: - print(f'Caught error while processing {filename}. Continuing anyways...\n{e}') + index += 1 def main(): @@ -218,30 +215,27 @@ def main(): for lat in lat_range: for lon in lon_range: tile_name = lat_lon_to_tile_string(lat, lon, is_worldcover=False) - try: - start_time = time.time() - extract_water( - processed_pbf_path, - lat, - lon, - tile_width, - tile_height, - interior_tile_dir=INTERIOR_TILE_DIR - ) - process_ocean_tiles( - args.ocean_polygons_path, - lat, - lon, - tile_width, - tile_height, - output_dir=OCEAN_TILE_DIR - ) - end_time = time.time() - total_time = end_time - start_time - print(f'Finished initial creation of {tile_name} in {total_time}(s). {index} of {num_tiles}') - index += 1 - except Exception as e: - print(f'Caught error while processing {tile_name}. Continuing anyways...\n{e}') + start_time = time.time() + extract_water( + processed_pbf_path, + lat, + lon, + tile_width, + tile_height, + interior_tile_dir=INTERIOR_TILE_DIR + ) + process_ocean_tiles( + args.ocean_polygons_path, + lat, + lon, + tile_width, + tile_height, + output_dir=OCEAN_TILE_DIR + ) + end_time = time.time() + total_time = end_time - start_time + print(f'Finished initial creation of {tile_name} in {total_time}(s). {index} of {num_tiles}') + index += 1 print('Merging processed tiles...') merge_interior_and_ocean(INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR, translate_to_cog=True) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 8dd4b4fa..f4f7019e 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -174,28 +174,25 @@ def crop_tile(tile): Args: tile: The filename of the desired tile to crop. """ - try: - ref_image = TILE_DIR + tile - out_filename = CROPPED_TILE_DIR + tile - pixel_size = gdal.Warp('tmp_px_size.tif', ref_image, dstSRS='EPSG:4326').GetGeoTransform()[1] - - shapefile_command = ' '.join(['gdaltindex', 'tmp.shp', ref_image]) - os.system(shapefile_command) - - gdal.Warp( - out_filename, - tile, - cutlineDSName='tmp.shp', - cropToCutline=True, - xRes=pixel_size, - yRes=pixel_size, - targetAlignedPixels=True, - dstSRS='EPSG:4326', - format='COG' - ) - remove_temp_files(['tmp_px_size.tif', 'tmp.shp']) - except Exception as e: # When a tile fails to process, we don't necessarily want the program to stop. - print(f'Could not process {tile} due to {e}. Skipping...') + ref_image = TILE_DIR + tile + out_filename = CROPPED_TILE_DIR + tile + pixel_size = gdal.Warp('tmp_px_size.tif', ref_image, dstSRS='EPSG:4326').GetGeoTransform()[1] + + shapefile_command = ' '.join(['gdaltindex', 'tmp.shp', ref_image]) + os.system(shapefile_command) + + gdal.Warp( + out_filename, + tile, + cutlineDSName='tmp.shp', + cropToCutline=True, + xRes=pixel_size, + yRes=pixel_size, + targetAlignedPixels=True, + dstSRS='EPSG:4326', + format='COG' + ) + remove_temp_files(['tmp_px_size.tif', 'tmp.shp']) def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py index 0cce7a86..9832e0aa 100644 --- a/src/asf_tools/watermasking/utils.py +++ b/src/asf_tools/watermasking/utils.py @@ -11,7 +11,7 @@ def lat_lon_to_tile_string(lat, lon, is_worldcover: bool = False, postfix: str = lon: The minimum longitude of the tile. is_worldcover: Wheter the tile is Worldcover or OSM. postfix: A postfix to append to the tile name to make it a filename. - + Returns: The name of the tile. """ @@ -45,22 +45,22 @@ def merge_tiles(tiles, out_format, out_filename): def remove_temp_files(temp_files: list): """Remove each file in a list of files. - + Args: temp_files: The list of temporary files to remove. """ for file in temp_files: try: os.remove(file) - except Exception as e: - print(f'Caught {e} while removing temporary file: {file}. Skipping it...') + except FileNotFoundError: + print(f'Temp file {file} was not found, skipping removal...') def setup_directories(dirs: list[str]): """Setup the directories necessary for running the script.""" - for dir in dirs: + for directory in dirs: try: - os.mkdir(dir) - except FileExistsError as e: + os.mkdir(directory) + except FileExistsError: # Directories already exists. pass From 49ce3301c951bf182e3ccdd00d3e2159d62bd542 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 14:15:13 -0600 Subject: [PATCH 27/58] basic instructions for the watermask scripts --- src/asf_tools/watermasking/README.MD | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 src/asf_tools/watermasking/README.MD diff --git a/src/asf_tools/watermasking/README.MD b/src/asf_tools/watermasking/README.MD new file mode 100644 index 00000000..037a8079 --- /dev/null +++ b/src/asf_tools/watermasking/README.MD @@ -0,0 +1,11 @@ +These scripts are for creating a global (or regional) water mask dataset based off of a OpenStreetMaps, and optionally augmented by ESA WorldCover. + +The basic steps are as follows: + Downloads: + 1. Download the "Latest Weekly Planet PBF File" file from here: "https://planet.openstreetmap.org/". + 2. Download the WGS84 water polygons shapefile from: "https://osmdata.openstreetmap.de/data/water-polygons.html". + 3. The files should be unzipped and you should have something like `planet.osm.pbf` or `planet.pbf` and `water_polygons.shp` (and the support files for `water_polygons.shp`). + Scripts: + 1. Run ```generate_osm_dataset --planet-file-path [path-to-planet.pbf] --ocean-polygons-path [path-to-water-polygons.shp] --lat-begin -85 --lat-end 85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` + 2. Run ```fill_missing_tiles --fill-value 0 --lat-begin -90 --lat-end -85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` + 3. Run ```fill_missing_tiles --fill-value 1 --lat-begin 85 --lat-end 90 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` \ No newline at end of file From 0acf55be5ff8bf7454ad32844f2521be1cab46fe Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 14:16:38 -0600 Subject: [PATCH 28/58] refactor --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index f4f7019e..2bbc1732 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -23,6 +23,7 @@ def tile_preprocessing(tile_dir, min_lat, max_lat): """ filenames = [f for f in os.listdir(tile_dir) if f.endswith('.tif')] + def filename_filter(filename): latitude = int(filename.split('_')[5][1:3]) (latitude >= min_lat) and (latitude <= max_lat) From fc904edd3cd87f831a479f359ec98d075d7801d7 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 14:17:22 -0600 Subject: [PATCH 29/58] refactor --- src/asf_tools/watermasking/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py index 9832e0aa..81ce4ed1 100644 --- a/src/asf_tools/watermasking/utils.py +++ b/src/asf_tools/watermasking/utils.py @@ -3,7 +3,7 @@ import numpy as np -def lat_lon_to_tile_string(lat, lon, is_worldcover: bool = False, postfix: str ='.tif'): +def lat_lon_to_tile_string(lat, lon, is_worldcover: bool = False, postfix: str = '.tif'): """Get the name of the tile with lower left corner (lat, lon). Args: From 17f1f5e5a41bd9449f6205d6b06d68b7e7a59654 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 14:17:30 -0600 Subject: [PATCH 30/58] better entrypoint names --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index a3aaac2a..56642ae9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,7 +41,8 @@ make_composite = "asf_tools.composite:main" water_map = "asf_tools.hydrosar.water_map:main" calculate_hand = "asf_tools.hydrosar.hand.calculate:main" flood_map = "asf_tools.hydrosar.flood_map:main" -generate_osm_tiles = "asf_tools.watermasking.generate_osm_tiles:main" +generate_osm_dataset = "asf_tools.watermasking.generate_osm_tiles:main" +generate_worldcover_dataset = "asf_tools.watermasking.generate_worldcover_tiles:main" fill_missing_tiles = "asf_tools.watermasking.fill_missing_tiles:main" [project.entry-points.hyp3] From 12c1e1617ac64eb5877f3343ec0f5af9e7a35644 Mon Sep 17 00:00:00 2001 From: Joseph H Kennedy Date: Wed, 28 Feb 2024 11:27:26 -0900 Subject: [PATCH 31/58] Don't verify PyPI distributions --- .github/workflows/distribute.yml | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/.github/workflows/distribute.yml b/.github/workflows/distribute.yml index bad4d8b1..798aedf0 100644 --- a/.github/workflows/distribute.yml +++ b/.github/workflows/distribute.yml @@ -33,22 +33,3 @@ jobs: with: user: __token__ password: ${{ secrets.TOOLS_PYPI_PAK }} - - verify-distribution: - runs-on: ubuntu-latest - needs: - - call-version-info-workflow - - distribute - defaults: - run: - shell: bash -l {0} - steps: - - uses: actions/checkout@v4 - - - uses: mamba-org/setup-micromamba@v1 - with: - environment-file: environment.yml - - - name: Ensure asf_tools v${{ needs.call-version-info-workflow.outputs.version }}} is pip installable - run: | - python -m pip install asf_tools==${{ needs.call-version-info-workflow.outputs.version_tag }} From 91e7bcc3a787475f06992d60859c4279c439f755 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 14:44:37 -0600 Subject: [PATCH 32/58] added filter by longitude --- .../watermasking/generate_worldcover_tiles.py | 23 +++++++++++++------ 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 2bbc1732..25440c2b 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -14,7 +14,7 @@ WORLDCOVER_TILE_SIZE = 3 -def tile_preprocessing(tile_dir, min_lat, max_lat): +def tile_preprocessing(tile_dir, min_lat, max_lat, min_lon, max_lon): """The worldcover tiles have lots of unnecessary classes, these need to be removed first. Note: make a back-up copy of this directory. @@ -26,7 +26,10 @@ def tile_preprocessing(tile_dir, min_lat, max_lat): def filename_filter(filename): latitude = int(filename.split('_')[5][1:3]) - (latitude >= min_lat) and (latitude <= max_lat) + longitude = int(filename.split('_')[5][4:7]) + in_lat_range = (latitude >= min_lat) and (latitude <= max_lat) + in_lon_range = (longitude >= min_lon) and (longitude <= max_lon) + return in_lat_range and in_lon_range filenames_filtered = [f for f in filenames if filename_filter(f)] index = 0 @@ -235,13 +238,19 @@ def main(): args = parser.parse_args() + lat_begin = int(args.lat_begin) + lat_end = int(args.lat_end) + lon_begin = int(args.lon_begin) + lon_end = int(args.lon_end) + tile_width = int(args.tile_width) + tile_height = int(args.tile_height) + lat_range = range(lat_begin, lat_end, tile_width) + lon_range = range(lon_begin, lon_end, tile_height) + setup_directories([PROCESSED_TILE_DIR, TILE_DIR, CROPPED_TILE_DIR]) # Process the multi-class masks into water/not-water masks. - tile_preprocessing(args.worldcover_tiles_dir, args.lat_begin, args.lat_end) - - lat_range = range(args.lat_begin, args.lat_end, args.tile_width) - lon_range = range(args.lon_begin, args.lon_end, args.tile_heigth) + tile_preprocessing(args.worldcover_tiles_dir, lat_begin, lat_end, lon_begin, lon_end) # Ocean only tiles are missing from WorldCover, so we need to create blank (water-only) ones. create_missing_tiles(PROCESSED_TILE_DIR, lat_range, lon_range) @@ -251,7 +260,7 @@ def main(): lat_range, lon_range, worldcover_degrees=WORLDCOVER_TILE_SIZE, - osm_degrees=args.tile_width + osm_degrees=tile_width ) From 2768eb243704a6368dfc6cc27b6becda6c380fac Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 16:09:26 -0600 Subject: [PATCH 33/58] fixes for wc generation --- .../watermasking/generate_worldcover_tiles.py | 76 +++++++++---------- src/asf_tools/watermasking/utils.py | 16 +++- 2 files changed, 51 insertions(+), 41 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 25440c2b..d40a9120 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -39,7 +39,7 @@ def filename_filter(filename): start_time = time.time() filename = tile_dir + filename - dst_filename = PROCESSED_TILE_DIR + filename.split('_')[5] + '.tif' + dst_filename = PROCESSED_TILE_DIR + filename.split('_')[6] + '.tif' print(f'Processing: {filename} --- {dst_filename} -- {index} of {num_tiles}') @@ -123,13 +123,13 @@ def create_missing_tiles(tile_dir, lat_range, lon_range): return current_existing_tiles -def get_tiles(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): +def get_tiles(osm_tile_coord: tuple, wc_tile_width: int, tile_width: int): """Get a list of the worldcover tile locations necessary to fully cover an OSM tile. Args: osm_tile_coord: The lower left corner coordinate (lat, lon) of the desired OSM tile. - wc_deg: The width/height of the Worldcover tiles in degrees. - osm_deg: The width/height of the OSM tiles in degrees. + wc_tile_width: The width/height of the Worldcover tiles in degrees. + tile_width: The width/height of the OSM tiles in degrees. Returns: tiles: A list of the lower left corner coordinates of the Worldcover tiles that overlap the OSM tile. @@ -138,13 +138,13 @@ def get_tiles(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): osm_lat = osm_tile_coord[0] osm_lon = osm_tile_coord[1] - min_lat = osm_lat - (osm_lat % wc_deg) - max_lat = osm_lat + osm_deg - min_lon = osm_lon - (osm_lon % wc_deg) - max_lon = osm_lon + osm_deg + min_lat = osm_lat - (osm_lat % wc_tile_width) + max_lat = osm_lat + tile_width + min_lon = osm_lon - (osm_lon % wc_tile_width) + max_lon = osm_lon + tile_width - lats = range(min_lat, max_lat, wc_deg) - lons = range(min_lon, max_lon, wc_deg) + lats = range(min_lat, max_lat, wc_tile_width) + lons = range(min_lon, max_lon, wc_tile_width) tiles = [] for lat in lats: @@ -154,52 +154,49 @@ def get_tiles(osm_tile_coord: tuple, wc_deg: int, osm_deg: int): return tiles -def lat_lon_to_filenames(worldcover_tile_dir, osm_tile_coord: tuple, wc_deg: int, osm_deg: int): +def lat_lon_to_filenames(worldcover_tile_dir, osm_tile_coord: tuple, wc_tile_width: int, tile_width: int): """Get a list of the Worldcover tile filenames that are necessary to overlap an OSM tile. Args: osm_tile: The lower left corner (lat, lon) of the desired OSM tile. - wc_deg: The width of the Worldcover tiles in degrees. - osm_deg: The width of the OSM tiles in degrees. + wc_tile_width: The width of the Worldcover tiles in degrees. + tile_width: The width of the OSM tiles in degrees. Returns: filenames: The list of Worldcover filenames. """ filenames = [] - tiles = get_tiles(osm_tile_coord, wc_deg, osm_deg) + tiles = get_tiles(osm_tile_coord, wc_tile_width, tile_width) for tile in tiles: filenames.append(worldcover_tile_dir + lat_lon_to_tile_string(tile[0], tile[1], is_worldcover=True)) return filenames -def crop_tile(tile): +def crop_tile(tile, lat, lon, tile_width, tile_height): """Crop the merged tiles Args: tile: The filename of the desired tile to crop. """ - ref_image = TILE_DIR + tile + in_filename = TILE_DIR + tile out_filename = CROPPED_TILE_DIR + tile - pixel_size = gdal.Warp('tmp_px_size.tif', ref_image, dstSRS='EPSG:4326').GetGeoTransform()[1] + pixel_size_x, pixel_size_y = 0.00009009009, -0.00009009009 - shapefile_command = ' '.join(['gdaltindex', 'tmp.shp', ref_image]) - os.system(shapefile_command) - - gdal.Warp( + src_ds = gdal.Open(in_filename) + gdal.Translate( out_filename, - tile, - cutlineDSName='tmp.shp', - cropToCutline=True, - xRes=pixel_size, - yRes=pixel_size, - targetAlignedPixels=True, - dstSRS='EPSG:4326', - format='COG' + src_ds, + projWin=[lon, lat+tile_height, lon+tile_width, lat], + xRes=pixel_size_x, + yRes=pixel_size_y, + outputSRS='EPSG:4326', + format='COG', + creationOptions=['NUM_THREADS=all_cpus'] ) remove_temp_files(['tmp_px_size.tif', 'tmp.shp']) -def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): +def build_dataset(worldcover_tile_dir, lat_range, lon_range, tile_width, tile_height): """ Main function for generating a dataset with worldcover tiles. Args: @@ -211,11 +208,12 @@ def build_dataset(worldcover_tile_dir, lat_range, lon_range, out_degrees): for lat in lat_range: for lon in lon_range: start_time = time.time() - tile_filename = TILE_DIR + lat_lon_to_tile_string(lat, lon, is_worldcover=False) - worldcover_tiles = lat_lon_to_filenames(worldcover_tile_dir, lat, lon, WORLDCOVER_TILE_SIZE, out_degrees) + tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False) + tile_filename = TILE_DIR + tile + worldcover_tiles = lat_lon_to_filenames(worldcover_tile_dir, (lat, lon), WORLDCOVER_TILE_SIZE, tile_width) print(f'Processing: {tile_filename} {worldcover_tiles}') - merge_tiles(worldcover_tiles, tile_filename) - crop_tile(tile_filename) + merge_tiles(worldcover_tiles, tile_filename, 'GTiff', compress=True) + crop_tile(tile, lat, lon, tile_width, tile_height) end_time = time.time() total_time = end_time - start_time print(f'Time Elapsed: {total_time}s') @@ -252,15 +250,17 @@ def main(): # Process the multi-class masks into water/not-water masks. tile_preprocessing(args.worldcover_tiles_dir, lat_begin, lat_end, lon_begin, lon_end) + wc_lat_range = range(lat_begin, lat_end, WORLDCOVER_TILE_SIZE) + wc_lon_range = range(lon_begin, lon_end, WORLDCOVER_TILE_SIZE) # Ocean only tiles are missing from WorldCover, so we need to create blank (water-only) ones. - create_missing_tiles(PROCESSED_TILE_DIR, lat_range, lon_range) + create_missing_tiles(PROCESSED_TILE_DIR, wc_lat_range, wc_lon_range) build_dataset( - args.worldcover_tile_dir, + PROCESSED_TILE_DIR, lat_range, lon_range, - worldcover_degrees=WORLDCOVER_TILE_SIZE, - osm_degrees=tile_width + tile_width=tile_width, + tile_height=tile_height ) diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py index 81ce4ed1..9b45482d 100644 --- a/src/asf_tools/watermasking/utils.py +++ b/src/asf_tools/watermasking/utils.py @@ -27,7 +27,7 @@ def lat_lon_to_tile_string(lat, lon, is_worldcover: bool = False, postfix: str = return lat_part + lon_part + postfix -def merge_tiles(tiles, out_format, out_filename): +def merge_tiles(tiles, out_filename, out_format, compress=False): """Merge tiles via buildvrt and translate. Args: @@ -37,10 +37,20 @@ def merge_tiles(tiles, out_format, out_filename): """ vrt = 'merged.vrt' build_vrt_command = ' '.join(['gdalbuildvrt', vrt] + tiles) - translate_command = ' '.join(['gdal_translate', '-of', out_format, vrt, out_filename]) + if not compress: + translate_command = ' '.join(['gdal_translate', '-of', out_format, vrt, out_filename]) + else: + translate_command = ' '.join([ + 'gdal_translate', + '-of', out_format, + '-co', 'COMPRESS=LZW', + '-co', 'NUM_THREADS=all_cpus', + vrt, + out_filename + ]) os.system(build_vrt_command) os.system(translate_command) - os.remove(vrt) + remove_temp_files([vrt]) def remove_temp_files(temp_files: list): From 18cbed2df1a0f9bba513d14852723a3c631e2d6c Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 16:11:04 -0600 Subject: [PATCH 34/58] added dependencies --- environment.yml | 2 ++ pyproject.toml | 1 + 2 files changed, 3 insertions(+) diff --git a/environment.yml b/environment.yml index b26f47ef..1f49abd9 100644 --- a/environment.yml +++ b/environment.yml @@ -25,6 +25,8 @@ dependencies: - gdal>=3.7 - geopandas - numpy + - osmium-tool + - pyogrio - pysheds>=0.3 - rasterio - scikit-fuzzy diff --git a/pyproject.toml b/pyproject.toml index 56642ae9..893453cf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,7 @@ dependencies = [ "gdal>=3.3", "geopandas", "numpy", + "pyogrio", "pysheds>=0.3", "rasterio", "scikit-fuzzy", From 6ac972ae9af840380baeeb4f5537c94c736b8495 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 16:13:47 -0600 Subject: [PATCH 35/58] remove whitespace --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index d40a9120..5a731592 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -183,7 +183,7 @@ def crop_tile(tile, lat, lon, tile_width, tile_height): pixel_size_x, pixel_size_y = 0.00009009009, -0.00009009009 src_ds = gdal.Open(in_filename) - gdal.Translate( + gdal.Translate( out_filename, src_ds, projWin=[lon, lat+tile_height, lon+tile_width, lat], From f91b2fa7789e0a99f283562420cb63ae5c694976 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 16:14:01 -0600 Subject: [PATCH 36/58] remove =None before del --- src/asf_tools/watermasking/fill_missing_tiles.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py index 0ba7c540..72d6e6a0 100644 --- a/src/asf_tools/watermasking/fill_missing_tiles.py +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -60,7 +60,6 @@ def main(): dst_ds.SetProjection(srs.ExportToWkt()) dst_band = dst_ds.GetRasterBand(1) dst_band.WriteArray(data) - dst_ds = None del dst_ds command = f'gdal_translate -of COG -co NUM_THREADS=all_cpus {tile_tif} {tile_cog}' From 95c291f4cb9e62cd67b7a661e2525baa98fbb71e Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Wed, 28 Feb 2024 16:17:21 -0600 Subject: [PATCH 37/58] updated changelog --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 281f3d5d..dd1d0757 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/) and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.7.0] + +## Added +* Scripts and entrypoints for generating our global watermasking dataset added to `watermasking`. + ## [0.6.0] From 3801c665dbaac544c475afdbf40f1d139fe2fa30 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 12:08:37 -0600 Subject: [PATCH 38/58] replaced system with subprocess --- src/asf_tools/watermasking/utils.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py index 9b45482d..7806be0d 100644 --- a/src/asf_tools/watermasking/utils.py +++ b/src/asf_tools/watermasking/utils.py @@ -1,4 +1,5 @@ import os +import subprocess import numpy as np @@ -36,20 +37,20 @@ def merge_tiles(tiles, out_filename, out_format, compress=False): out_filename: The name of the output COG. """ vrt = 'merged.vrt' - build_vrt_command = ' '.join(['gdalbuildvrt', vrt] + tiles) + build_vrt_command = ['gdalbuildvrt', vrt] + tiles if not compress: - translate_command = ' '.join(['gdal_translate', '-of', out_format, vrt, out_filename]) + translate_command = ['gdal_translate', '-of', out_format, vrt, out_filename] else: - translate_command = ' '.join([ + translate_command = [ 'gdal_translate', '-of', out_format, '-co', 'COMPRESS=LZW', '-co', 'NUM_THREADS=all_cpus', vrt, out_filename - ]) - os.system(build_vrt_command) - os.system(translate_command) + ] + subprocess.run(build_vrt_command) + subprocess.run(translate_command) remove_temp_files([vrt]) From 4b3c9bdfb4584c3a90b35e1738e304414729868b Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 12:09:05 -0600 Subject: [PATCH 39/58] refactoring --- .../watermasking/generate_worldcover_tiles.py | 29 +++++++++++-------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 5a731592..8f5eb371 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -5,13 +5,17 @@ import numpy as np from osgeo import gdal +gdal.UseExceptions() + from asf_tools.watermasking.utils import lat_lon_to_tile_string, merge_tiles, remove_temp_files, setup_directories -PROCESSED_TILE_DIR = 'worldcover_tiles_preprocessed/' -TILE_DIR = 'worldcover_tiles_uncropped/' + +PREPROCESSED_TILE_DIR = 'worldcover_tiles_preprocessed/' +UNCROPPED_TILE_DIR = 'worldcover_tiles_uncropped/' CROPPED_TILE_DIR = 'worldcover_tiles/' FILENAME_POSTFIX = '.tif' WORLDCOVER_TILE_SIZE = 3 +GDAL_OPTIONS = ['COMPRESS=LZW', 'TILED=YES', 'NUM_THREADS=all_cpus'] def tile_preprocessing(tile_dir, min_lat, max_lat, min_lon, max_lon): @@ -39,7 +43,7 @@ def filename_filter(filename): start_time = time.time() filename = tile_dir + filename - dst_filename = PROCESSED_TILE_DIR + filename.split('_')[6] + '.tif' + dst_filename = PREPROCESSED_TILE_DIR + filename.split('_')[6] + '.tif' print(f'Processing: {filename} --- {dst_filename} -- {index} of {num_tiles}') @@ -52,14 +56,13 @@ def filename_filter(filename): water_arr[not_water] = 0 driver = gdal.GetDriverByName('GTiff') - dst_ds = driver.Create( dst_filename, water_arr.shape[0], water_arr.shape[1], 1, gdal.GDT_Byte, - options=['COMPRESS=LZW', 'TILED=YES'] + options=GDAL_OPTIONS ) dst_ds.SetGeoTransform(src_ds.GetGeoTransform()) dst_ds.SetProjection(src_ds.GetProjection()) @@ -95,6 +98,7 @@ def create_missing_tiles(tile_dir, lat_range, lon_range): if tile not in current_existing_tiles: print(f'Could not find {tile}') + filename = PREPROCESSED_TILE_DIR + tile x_size, y_size = 36000, 36000 x_res, y_res = 8.333333333333333055e-05, -8.333333333333333055e-05 ul_lon = lon @@ -103,12 +107,12 @@ def create_missing_tiles(tile_dir, lat_range, lon_range): driver = gdal.GetDriverByName('GTiff') ds = driver.Create( - tile, + filename, xsize=x_size, ysize=y_size, bands=1, eType=gdal.GDT_Byte, - options=['COMPRESS=LZW'] + options=GDAL_OPTIONS ) ds.SetProjection('EPSG:4326') ds.SetGeoTransform(geotransform) @@ -178,7 +182,7 @@ def crop_tile(tile, lat, lon, tile_width, tile_height): Args: tile: The filename of the desired tile to crop. """ - in_filename = TILE_DIR + tile + in_filename = UNCROPPED_TILE_DIR + tile out_filename = CROPPED_TILE_DIR + tile pixel_size_x, pixel_size_y = 0.00009009009, -0.00009009009 @@ -209,7 +213,7 @@ def build_dataset(worldcover_tile_dir, lat_range, lon_range, tile_width, tile_he for lon in lon_range: start_time = time.time() tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False) - tile_filename = TILE_DIR + tile + tile_filename = UNCROPPED_TILE_DIR + tile worldcover_tiles = lat_lon_to_filenames(worldcover_tile_dir, (lat, lon), WORLDCOVER_TILE_SIZE, tile_width) print(f'Processing: {tile_filename} {worldcover_tiles}') merge_tiles(worldcover_tiles, tile_filename, 'GTiff', compress=True) @@ -245,18 +249,19 @@ def main(): lat_range = range(lat_begin, lat_end, tile_width) lon_range = range(lon_begin, lon_end, tile_height) - setup_directories([PROCESSED_TILE_DIR, TILE_DIR, CROPPED_TILE_DIR]) + setup_directories([PREPROCESSED_TILE_DIR, UNCROPPED_TILE_DIR, CROPPED_TILE_DIR]) # Process the multi-class masks into water/not-water masks. tile_preprocessing(args.worldcover_tiles_dir, lat_begin, lat_end, lon_begin, lon_end) wc_lat_range = range(lat_begin, lat_end, WORLDCOVER_TILE_SIZE) wc_lon_range = range(lon_begin, lon_end, WORLDCOVER_TILE_SIZE) + # Ocean only tiles are missing from WorldCover, so we need to create blank (water-only) ones. - create_missing_tiles(PROCESSED_TILE_DIR, wc_lat_range, wc_lon_range) + create_missing_tiles(PREPROCESSED_TILE_DIR, wc_lat_range, wc_lon_range) build_dataset( - PROCESSED_TILE_DIR, + PREPROCESSED_TILE_DIR, lat_range, lon_range, tile_width=tile_width, From 71bae1d7409262f40c546b72debfd2f2430d48cf Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 12:14:48 -0600 Subject: [PATCH 40/58] switch to subprocess --- .../watermasking/generate_osm_tiles.py | 39 +++++++++++-------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index 7bdcdaba..d3e6e695 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -1,15 +1,20 @@ import argparse import os +import subprocess import time import geopandas as gpd from osgeo import gdal +gdal.UseExceptions() + from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, setup_directories + INTERIOR_TILE_DIR = 'interior_tiles/' OCEAN_TILE_DIR = 'ocean_tiles/' FINISHED_TILE_DIR = 'tiles/' +GDAL_OPTIONS = ['COMPRESS=LZW', 'TILED=YES', 'NUM_THREADS=all_cpus'] def process_pbf(planet_file: str, output_file: str): @@ -24,15 +29,15 @@ def process_pbf(planet_file: str, output_file: str): waterways_file = 'planet_waterways.pbf' reservoirs_file = 'planet_reservoirs.pbf' - natural_water_command = f'osmium tags-filter -o {natural_file} {planet_file} wr/natural=water' - waterways_command = f'osmium tags-filter -o {waterways_file} {planet_file} waterway="*"' - reservoirs_command = f'osmium tags-filter -o {reservoirs_file} {planet_file} landuse=reservoir' - merge_command = f'osmium merge {natural_file} {waterways_file} {reservoirs_file} -o {output_file}' + natural_water_command = f'osmium tags-filter -o {natural_file} {planet_file} wr/natural=water'.split(' ') + waterways_command = f'osmium tags-filter -o {waterways_file} {planet_file} waterway="*"'.split(' ') + reservoirs_command = f'osmium tags-filter -o {reservoirs_file} {planet_file} landuse=reservoir'.split(' ') + merge_command = f'osmium merge {natural_file} {waterways_file} {reservoirs_file} -o {output_file}'.split(' ') - os.system(natural_water_command) - os.system(waterways_command) - os.system(reservoirs_command) - os.system(merge_command) + subprocess.run(natural_water_command) + subprocess.run(waterways_command) + subprocess.run(reservoirs_command) + subprocess.run(merge_command) def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_height_deg, output_dir): @@ -54,8 +59,8 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig pixel_size_y = 0.00009009009 clipped_polygons_path = tile + '.shp' - command = f'ogr2ogr -clipsrc {xmin} {ymin} {xmax} {ymax} {clipped_polygons_path} {ocean_polygons_path}' - os.system(command) + command = f'ogr2ogr -clipsrc {xmin} {ymin} {xmax} {ymax} {clipped_polygons_path} {ocean_polygons_path}'.split(' ') + subprocess.run(command) gdal.Rasterize( tile_tif, @@ -65,7 +70,7 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], outputType=gdal.GDT_Byte, - creationOptions=['COMPRESS=LZW'] + creationOptions=GDAL_OPTIONS ) temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx'] @@ -91,8 +96,10 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio # Extract tile from the main pbf, then convert it to a tif. bbox = f'--bbox {lon},{lat},{lon+tile_width_deg},{lat+tile_height_deg}' - os.system(f'osmium extract -s smart -S tags=natural=water {bbox} {water_file} -o {tile_pbf}') - os.system(f'osmium export --geometry-types="polygon" {tile_pbf} -o {tile_geojson}') + extract_command = f'osmium extract -s smart -S tags=natural=water {bbox} {water_file} -o {tile_pbf}'.split(' ') + export_command = f'osmium export --geometry-types="polygon" {tile_pbf} -o {tile_geojson}'.split(' ') + subprocess.run(extract_command) + subprocess.run(export_command) # Islands and Islets can be members of the water features, so they must be removed. water_gdf = gpd.read_file(tile_geojson, engine='pyogrio') @@ -117,7 +124,7 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], outputType=gdal.GDT_Byte, - creationOptions=['COMPRESS=LZW'] + creationOptions=GDAL_OPTIONS ) temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx', tile_shp, tile_pbf, tile_geojson] @@ -154,7 +161,7 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di '--calc', '"logical_or(A, B)"' ]) - os.system(command) + subprocess.run(command) if translate_to_cog: cogs_dir = finished_tile_dir + 'cogs/' @@ -164,7 +171,7 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di pass out_file = cogs_dir + filename command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus {output_tile} {out_file}' - os.system(command) + subprocess.run(command) os.remove(output_tile) end_time = time.time() From be67b7084ec0e8652ca41b96af245796db5732f5 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 12:16:36 -0600 Subject: [PATCH 41/58] added gdal.UseExceptions --- src/asf_tools/watermasking/fill_missing_tiles.py | 3 +++ src/asf_tools/watermasking/generate_osm_tiles.py | 4 ++-- src/asf_tools/watermasking/generate_worldcover_tiles.py | 4 ++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py index 72d6e6a0..656cd4e2 100644 --- a/src/asf_tools/watermasking/fill_missing_tiles.py +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -7,6 +7,9 @@ from asf_tools.watermasking.utils import lat_lon_to_tile_string +gdal.UseExceptions() + + def main(): parser = argparse.ArgumentParser( diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index d3e6e695..bfc5dc04 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -6,11 +6,11 @@ import geopandas as gpd from osgeo import gdal -gdal.UseExceptions() - from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, setup_directories +gdal.UseExceptions() + INTERIOR_TILE_DIR = 'interior_tiles/' OCEAN_TILE_DIR = 'ocean_tiles/' FINISHED_TILE_DIR = 'tiles/' diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 8f5eb371..12f87324 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -5,11 +5,11 @@ import numpy as np from osgeo import gdal -gdal.UseExceptions() - from asf_tools.watermasking.utils import lat_lon_to_tile_string, merge_tiles, remove_temp_files, setup_directories +gdal.UseExceptions() + PREPROCESSED_TILE_DIR = 'worldcover_tiles_preprocessed/' UNCROPPED_TILE_DIR = 'worldcover_tiles_uncropped/' CROPPED_TILE_DIR = 'worldcover_tiles/' From 129cd0367200c824611df5fc706f5c5c790cbac8 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 12:45:36 -0600 Subject: [PATCH 42/58] instructions for worldcover --- src/asf_tools/watermasking/README.MD | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/asf_tools/watermasking/README.MD b/src/asf_tools/watermasking/README.MD index 037a8079..0edebdd4 100644 --- a/src/asf_tools/watermasking/README.MD +++ b/src/asf_tools/watermasking/README.MD @@ -1,11 +1,15 @@ These scripts are for creating a global (or regional) water mask dataset based off of a OpenStreetMaps, and optionally augmented by ESA WorldCover. -The basic steps are as follows: - Downloads: - 1. Download the "Latest Weekly Planet PBF File" file from here: "https://planet.openstreetmap.org/". - 2. Download the WGS84 water polygons shapefile from: "https://osmdata.openstreetmap.de/data/water-polygons.html". - 3. The files should be unzipped and you should have something like `planet.osm.pbf` or `planet.pbf` and `water_polygons.shp` (and the support files for `water_polygons.shp`). - Scripts: - 1. Run ```generate_osm_dataset --planet-file-path [path-to-planet.pbf] --ocean-polygons-path [path-to-water-polygons.shp] --lat-begin -85 --lat-end 85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` - 2. Run ```fill_missing_tiles --fill-value 0 --lat-begin -90 --lat-end -85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` - 3. Run ```fill_missing_tiles --fill-value 1 --lat-begin 85 --lat-end 90 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` \ No newline at end of file +For the OSM water mask dataset, follow these steps to replicate our dataset: + 1. Download the "Latest Weekly Planet PBF File" file from here: "https://planet.openstreetmap.org/". + 2. Download the WGS84 water polygons shapefile from: "https://osmdata.openstreetmap.de/data/water-polygons.html". + 3. The files should be unzipped and you should have something like `planet.osm.pbf` or `planet.pbf` and `water_polygons.shp` (and the support files for `water_polygons.shp`). + 4. Run ```generate_osm_dataset --planet-file-path [path-to-planet.pbf] --ocean-polygons-path [path-to-water-polygons.shp] --lat-begin -85 --lat-end 85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` + 5. Run ```fill_missing_tiles --fill-value 0 --lat-begin -90 --lat-end -85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` + 6. Run ```fill_missing_tiles --fill-value 1 --lat-begin 85 --lat-end 90 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` + +For the WorldCover water mask dataset, follow these steps: + 1. Download the portions of the dataset for the areas you would like to cover from here: "https://worldcover2020.esa.int/downloader" + 2. Extract the contents into a folder. Note, if you download multiple portions of the dataset, extract them all into the same folder. + 3. Run ```generate_worldcover_dataset --worldcover-tiles-dir [path-to-worldcover-data] --lat-begin 55 --lat-end 80 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` + Note that we only use WorldCover data over Alaska, Canada, and Russia for our dataset. \ No newline at end of file From 8c735c082e15e7c774a694592f6e39b5053a9ff5 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 12:46:14 -0600 Subject: [PATCH 43/58] typo --- src/asf_tools/watermasking/README.MD | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/asf_tools/watermasking/README.MD b/src/asf_tools/watermasking/README.MD index 0edebdd4..537b2d60 100644 --- a/src/asf_tools/watermasking/README.MD +++ b/src/asf_tools/watermasking/README.MD @@ -1,4 +1,4 @@ -These scripts are for creating a global (or regional) water mask dataset based off of a OpenStreetMaps, and optionally augmented by ESA WorldCover. +These scripts are for creating a global (or regional) water mask dataset based off of OpenStreetMaps, and optionally augmented by ESA WorldCover. For the OSM water mask dataset, follow these steps to replicate our dataset: 1. Download the "Latest Weekly Planet PBF File" file from here: "https://planet.openstreetmap.org/". From 63b4f6b7f3ce660452125064ffe48cc7290768fc Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 13:21:20 -0600 Subject: [PATCH 44/58] removed quotes around polygon --- src/asf_tools/watermasking/generate_osm_tiles.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index bfc5dc04..cf6a2f9d 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -97,7 +97,7 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio # Extract tile from the main pbf, then convert it to a tif. bbox = f'--bbox {lon},{lat},{lon+tile_width_deg},{lat+tile_height_deg}' extract_command = f'osmium extract -s smart -S tags=natural=water {bbox} {water_file} -o {tile_pbf}'.split(' ') - export_command = f'osmium export --geometry-types="polygon" {tile_pbf} -o {tile_geojson}'.split(' ') + export_command = f'osmium export --geometry-types=polygon {tile_pbf} -o {tile_geojson}'.split(' ') subprocess.run(extract_command) subprocess.run(export_command) @@ -210,9 +210,9 @@ def main(): setup_directories([INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR]) - print('Extracting water from planet file...') + # print('Extracting water from planet file...') processed_pbf_path = 'planet_processed.pbf' - process_pbf(args.planet_file_path, processed_pbf_path) + # process_pbf(args.planet_file_path, processed_pbf_path) print('Processing tiles...') lat_range = range(lat_begin, lat_end, tile_height) From 47ead4d8ee1f4cd93a44f9796e1aba8250389b0d Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 13:52:25 -0600 Subject: [PATCH 45/58] fix line length --- src/asf_tools/watermasking/generate_osm_tiles.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index cf6a2f9d..c6d54ce7 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -73,7 +73,7 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig creationOptions=GDAL_OPTIONS ) - temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx'] + temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx', tile + '.shp'] remove_temp_files(temp_files) @@ -148,7 +148,7 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di internal_tile = internal_tile_dir + filename external_tile = ocean_tile_dir + filename output_tile = finished_tile_dir + filename - command = ' '.join([ + command = [ 'gdal_calc.py', '-A', internal_tile, @@ -160,7 +160,7 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di output_tile, '--calc', '"logical_or(A, B)"' - ]) + ] subprocess.run(command) if translate_to_cog: @@ -170,7 +170,8 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di except FileExistsError: pass out_file = cogs_dir + filename - command = f'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus {output_tile} {out_file}' + translate_string = 'gdal_translate -ot Byte -of COG -co NUM_THREADS=all_cpus' + command = f'{translate_string} {output_tile} {out_file}'.split(' ') subprocess.run(command) os.remove(output_tile) @@ -210,9 +211,9 @@ def main(): setup_directories([INTERIOR_TILE_DIR, OCEAN_TILE_DIR, FINISHED_TILE_DIR]) - # print('Extracting water from planet file...') + print('Extracting water from planet file...') processed_pbf_path = 'planet_processed.pbf' - # process_pbf(args.planet_file_path, processed_pbf_path) + process_pbf(args.planet_file_path, processed_pbf_path) print('Processing tiles...') lat_range = range(lat_begin, lat_end, tile_height) From d762ba97cfb4ab4aee34c086f06a85c31a286cd4 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 13:54:46 -0600 Subject: [PATCH 46/58] replaced system with subprocess --- src/asf_tools/watermasking/fill_missing_tiles.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py index 656cd4e2..6db3d2f8 100644 --- a/src/asf_tools/watermasking/fill_missing_tiles.py +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -1,5 +1,6 @@ import argparse import os +import subprocess import numpy as np from osgeo import gdal, osr @@ -65,8 +66,8 @@ def main(): dst_band.WriteArray(data) del dst_ds - command = f'gdal_translate -of COG -co NUM_THREADS=all_cpus {tile_tif} {tile_cog}' - os.system(command) + command = f'gdal_translate -of COG -co NUM_THREADS=all_cpus {tile_tif} {tile_cog}'.split(' ') + subprocess.run(command) os.remove(tile_tif) From ed5f543b5b0ee5d77baf45d5ca40bec44c5a69a1 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 14:53:49 -0600 Subject: [PATCH 47/58] clarifying comments --- src/asf_tools/watermasking/fill_missing_tiles.py | 8 ++++---- src/asf_tools/watermasking/generate_osm_tiles.py | 12 ++++++------ .../watermasking/generate_worldcover_tiles.py | 8 ++++---- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py index 6db3d2f8..ea9b64a2 100644 --- a/src/asf_tools/watermasking/fill_missing_tiles.py +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -19,10 +19,10 @@ def main(): ) parser.add_argument('--fill-value', help='The value to fill the data array with.', default=0) - parser.add_argument('--lat-begin', help='The minimum latitude of the dataset.', default=-85) - parser.add_argument('--lat-end', help='The maximum latitude of the dataset.', default=85) - parser.add_argument('--lon-begin', help='The minimum longitude of the dataset.', default=-180) - parser.add_argument('--lon-end', help='The maximum longitude of the dataset.', default=180) + parser.add_argument('--lat-begin', help='The minimum latitude of the dataset in EPSG:4326.', default=-85) + parser.add_argument('--lat-end', help='The maximum latitude of the dataset in EPSG:4326.', default=85) + parser.add_argument('--lon-begin', help='The minimum longitude of the dataset in EPSG:4326.', default=-180) + parser.add_argument('--lon-end', help='The maximum longitude of the dataset in EPSG:4326.', default=180) parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index c6d54ce7..d9d391b6 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -14,7 +14,7 @@ INTERIOR_TILE_DIR = 'interior_tiles/' OCEAN_TILE_DIR = 'ocean_tiles/' FINISHED_TILE_DIR = 'tiles/' -GDAL_OPTIONS = ['COMPRESS=LZW', 'TILED=YES', 'NUM_THREADS=all_cpus'] +GDAL_OPTIONS = ['COMPRESS=LZW', 'NUM_THREADS=all_cpus'] def process_pbf(planet_file: str, output_file: str): @@ -64,7 +64,7 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig gdal.Rasterize( tile_tif, - clipped_polygons_path, + ocean_polygons_path, xRes=pixel_size_x, yRes=pixel_size_y, burnValues=1, @@ -193,10 +193,10 @@ def main(): parser.add_argument('--planet-file-path', help='The path to the global planet.pbf file.') parser.add_argument('--ocean-polygons-path', help='The path to the global OSM ocean polygons.') - parser.add_argument('--lat-begin', help='The minimum latitude of the dataset.', default=-85) - parser.add_argument('--lat-end', help='The maximum latitude of the dataset.', default=85) - parser.add_argument('--lon-begin', help='The minimum longitude of the dataset.', default=-180) - parser.add_argument('--lon-end', help='The maximum longitude of the dataset.', default=180) + parser.add_argument('--lat-begin', help='The minimum latitude of the dataset in EPSG:4326.', default=-85) + parser.add_argument('--lat-end', help='The maximum latitude of the dataset in EPSG:4326.', default=85) + parser.add_argument('--lon-begin', help='The minimum longitude of the dataset in EPSG:4326.', default=-180) + parser.add_argument('--lon-end', help='The maximum longitude of the dataset in EPSG:4326.', default=180) parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 12f87324..3d81b1f2 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -231,10 +231,10 @@ def main(): ) parser.add_argument('--worldcover-tiles-dir', help='The path to the directory containing the worldcover tifs.') - parser.add_argument('--lat-begin', help='The minimum latitude of the dataset.', default=-85, required=True) - parser.add_argument('--lat-end', help='The maximum latitude of the dataset.', default=85) - parser.add_argument('--lon-begin', help='The minimum longitude of the dataset.', default=-180) - parser.add_argument('--lon-end', help='The maximum longitude of the dataset.', default=180) + parser.add_argument('--lat-begin', help='The minimum latitude of the dataset in EPSG:4326.', default=-85, required=True) + parser.add_argument('--lat-end', help='The maximum latitude of the dataset in EPSG:4326.', default=85) + parser.add_argument('--lon-begin', help='The minimum longitude of the dataset in EPSG:4326.', default=-180) + parser.add_argument('--lon-end', help='The maximum longitude of the dataset in EPSG:4326.', default=180) parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) From 7a6be8b121cd2fe4b29ed3d202513b162e61bbbf Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 15:07:46 -0600 Subject: [PATCH 48/58] expanded filter --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 3d81b1f2..e289ec82 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -31,8 +31,8 @@ def tile_preprocessing(tile_dir, min_lat, max_lat, min_lon, max_lon): def filename_filter(filename): latitude = int(filename.split('_')[5][1:3]) longitude = int(filename.split('_')[5][4:7]) - in_lat_range = (latitude >= min_lat) and (latitude <= max_lat) - in_lon_range = (longitude >= min_lon) and (longitude <= max_lon) + in_lat_range = (latitude >= min_lat - WORLDCOVER_TILE_SIZE) and (latitude <= max_lat + WORLDCOVER_TILE_SIZE) + in_lon_range = (longitude >= min_lon - WORLDCOVER_TILE_SIZE) and (longitude <= max_lon + WORLDCOVER_TILE_SIZE) return in_lat_range and in_lon_range filenames_filtered = [f for f in filenames if filename_filter(f)] From f47f7af195f95e4a365c3f2da5d4d828bf07f621 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 15:19:50 -0600 Subject: [PATCH 49/58] flake8 --- .../watermasking/generate_worldcover_tiles.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index e289ec82..92f1d5e2 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -28,11 +28,15 @@ def tile_preprocessing(tile_dir, min_lat, max_lat, min_lon, max_lon): filenames = [f for f in os.listdir(tile_dir) if f.endswith('.tif')] - def filename_filter(filename): + def filename_filter(filename): # We only want to preprocess the files that are necessary for the roi. latitude = int(filename.split('_')[5][1:3]) longitude = int(filename.split('_')[5][4:7]) - in_lat_range = (latitude >= min_lat - WORLDCOVER_TILE_SIZE) and (latitude <= max_lat + WORLDCOVER_TILE_SIZE) - in_lon_range = (longitude >= min_lon - WORLDCOVER_TILE_SIZE) and (longitude <= max_lon + WORLDCOVER_TILE_SIZE) + mnlat = min_lat - (min_lat % WORLDCOVER_TILE_SIZE) + mnlon = min_lon - (min_lon % WORLDCOVER_TILE_SIZE) + mxlat = max_lat + (max_lat % WORLDCOVER_TILE_SIZE) + mxlon = max_lon + (max_lon % WORLDCOVER_TILE_SIZE) + in_lat_range = (latitude >= mnlat) and (latitude <= mxlat) + in_lon_range = (longitude >= mnlon) and (longitude <= mxlon) return in_lat_range and in_lon_range filenames_filtered = [f for f in filenames if filename_filter(f)] @@ -231,7 +235,12 @@ def main(): ) parser.add_argument('--worldcover-tiles-dir', help='The path to the directory containing the worldcover tifs.') - parser.add_argument('--lat-begin', help='The minimum latitude of the dataset in EPSG:4326.', default=-85, required=True) + parser.add_argument( + '--lat-begin', + help='The minimum latitude of the dataset in EPSG:4326.', + default=-85, + required=True + ) parser.add_argument('--lat-end', help='The maximum latitude of the dataset in EPSG:4326.', default=85) parser.add_argument('--lon-begin', help='The minimum longitude of the dataset in EPSG:4326.', default=-180) parser.add_argument('--lon-end', help='The maximum longitude of the dataset in EPSG:4326.', default=180) From 478b7f8ead9a26239b312f39ae5dba1561be751b Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 15:45:55 -0600 Subject: [PATCH 50/58] better filters --- .../watermasking/generate_worldcover_tiles.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 92f1d5e2..0ba403d3 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -27,8 +27,7 @@ def tile_preprocessing(tile_dir, min_lat, max_lat, min_lon, max_lon): """ filenames = [f for f in os.listdir(tile_dir) if f.endswith('.tif')] - - def filename_filter(filename): # We only want to preprocess the files that are necessary for the roi. + def filename_filter(filename): latitude = int(filename.split('_')[5][1:3]) longitude = int(filename.split('_')[5][4:7]) mnlat = min_lat - (min_lat % WORLDCOVER_TILE_SIZE) @@ -263,8 +262,16 @@ def main(): # Process the multi-class masks into water/not-water masks. tile_preprocessing(args.worldcover_tiles_dir, lat_begin, lat_end, lon_begin, lon_end) - wc_lat_range = range(lat_begin, lat_end, WORLDCOVER_TILE_SIZE) - wc_lon_range = range(lon_begin, lon_end, WORLDCOVER_TILE_SIZE) + wc_lat_range = range( + lat_begin - (lat_begin % WORLDCOVER_TILE_SIZE), + lat_end + (lat_end % WORLDCOVER_TILE_SIZE), + WORLDCOVER_TILE_SIZE + ) + wc_lon_range = range( + lon_begin - (lon_begin % WORLDCOVER_TILE_SIZE), + lon_end + (lon_end % WORLDCOVER_TILE_SIZE), + WORLDCOVER_TILE_SIZE + ) # Ocean only tiles are missing from WorldCover, so we need to create blank (water-only) ones. create_missing_tiles(PREPROCESSED_TILE_DIR, wc_lat_range, wc_lon_range) From 9b2e13c1f822ca3301574ed9dc3277dbe5676d31 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 15:48:20 -0600 Subject: [PATCH 51/58] uncommented line --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 0ba403d3..ff4a90d0 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -30,6 +30,7 @@ def tile_preprocessing(tile_dir, min_lat, max_lat, min_lon, max_lon): def filename_filter(filename): latitude = int(filename.split('_')[5][1:3]) longitude = int(filename.split('_')[5][4:7]) + if filename.split('_')[5][3] == 'W': longitude = -longitude mnlat = min_lat - (min_lat % WORLDCOVER_TILE_SIZE) mnlon = min_lon - (min_lon % WORLDCOVER_TILE_SIZE) mxlat = max_lat + (max_lat % WORLDCOVER_TILE_SIZE) From 871267653d0779de360c067a32a19f3bcec67c95 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 15:52:07 -0600 Subject: [PATCH 52/58] removed whitespace --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index ff4a90d0..25223927 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -269,8 +269,8 @@ def main(): WORLDCOVER_TILE_SIZE ) wc_lon_range = range( - lon_begin - (lon_begin % WORLDCOVER_TILE_SIZE), - lon_end + (lon_end % WORLDCOVER_TILE_SIZE), + lon_begin - (lon_begin % WORLDCOVER_TILE_SIZE), + lon_end + (lon_end % WORLDCOVER_TILE_SIZE), WORLDCOVER_TILE_SIZE ) From 576545013414387ddba4a0cf12239e0625329fb9 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 16:08:27 -0600 Subject: [PATCH 53/58] final touches --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 25223927..10619e62 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -46,8 +46,9 @@ def filename_filter(filename): start_time = time.time() + tile_name = filename.split('_')[5] filename = tile_dir + filename - dst_filename = PREPROCESSED_TILE_DIR + filename.split('_')[6] + '.tif' + dst_filename = PREPROCESSED_TILE_DIR + tile_name + '.tif' print(f'Processing: {filename} --- {dst_filename} -- {index} of {num_tiles}') From fa54a43e48d881e4c41eb9e127a1d794bf52ea43 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Thu, 29 Feb 2024 16:11:08 -0600 Subject: [PATCH 54/58] flake8 --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 10619e62..88d3fa36 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -27,10 +27,12 @@ def tile_preprocessing(tile_dir, min_lat, max_lat, min_lon, max_lon): """ filenames = [f for f in os.listdir(tile_dir) if f.endswith('.tif')] + def filename_filter(filename): latitude = int(filename.split('_')[5][1:3]) longitude = int(filename.split('_')[5][4:7]) - if filename.split('_')[5][3] == 'W': longitude = -longitude + if filename.split('_')[5][3] == 'W': + longitude = -longitude mnlat = min_lat - (min_lat % WORLDCOVER_TILE_SIZE) mnlon = min_lon - (min_lon % WORLDCOVER_TILE_SIZE) mxlat = max_lat + (max_lat % WORLDCOVER_TILE_SIZE) From 9981861c027b9153b7f64a60b8a203fec1adbdf1 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 1 Mar 2024 14:36:09 +0000 Subject: [PATCH 55/58] fix pathing issue --- src/asf_tools/watermasking/generate_worldcover_tiles.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 88d3fa36..508f8883 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -1,6 +1,7 @@ import argparse import os import time +from pathlib import Path import numpy as np from osgeo import gdal @@ -49,7 +50,7 @@ def filename_filter(filename): start_time = time.time() tile_name = filename.split('_')[5] - filename = tile_dir + filename + filename = str(Path(tile_dir) / filename) dst_filename = PREPROCESSED_TILE_DIR + tile_name + '.tif' print(f'Processing: {filename} --- {dst_filename} -- {index} of {num_tiles}') From 2352660d4cf4ac643a9847150a367267e0d660d5 Mon Sep 17 00:00:00 2001 From: Forrest Williams <31411324+forrestfwilliams@users.noreply.github.com> Date: Fri, 1 Mar 2024 09:06:51 -0600 Subject: [PATCH 56/58] Update README.MD --- src/asf_tools/watermasking/README.MD | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/asf_tools/watermasking/README.MD b/src/asf_tools/watermasking/README.MD index 537b2d60..8ef49a3c 100644 --- a/src/asf_tools/watermasking/README.MD +++ b/src/asf_tools/watermasking/README.MD @@ -1,6 +1,7 @@ These scripts are for creating a global (or regional) water mask dataset based off of OpenStreetMaps, and optionally augmented by ESA WorldCover. For the OSM water mask dataset, follow these steps to replicate our dataset: + 1. Download the "Latest Weekly Planet PBF File" file from here: "https://planet.openstreetmap.org/". 2. Download the WGS84 water polygons shapefile from: "https://osmdata.openstreetmap.de/data/water-polygons.html". 3. The files should be unzipped and you should have something like `planet.osm.pbf` or `planet.pbf` and `water_polygons.shp` (and the support files for `water_polygons.shp`). @@ -9,7 +10,8 @@ For the OSM water mask dataset, follow these steps to replicate our dataset: 6. Run ```fill_missing_tiles --fill-value 1 --lat-begin 85 --lat-end 90 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` For the WorldCover water mask dataset, follow these steps: + 1. Download the portions of the dataset for the areas you would like to cover from here: "https://worldcover2020.esa.int/downloader" 2. Extract the contents into a folder. Note, if you download multiple portions of the dataset, extract them all into the same folder. 3. Run ```generate_worldcover_dataset --worldcover-tiles-dir [path-to-worldcover-data] --lat-begin 55 --lat-end 80 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` - Note that we only use WorldCover data over Alaska, Canada, and Russia for our dataset. \ No newline at end of file + Note that we only use WorldCover data over Alaska, Canada, and Russia for our dataset. From 3022127cc9e827d8235489cdc09daf9c7f8be71c Mon Sep 17 00:00:00 2001 From: Forrest Williams <31411324+forrestfwilliams@users.noreply.github.com> Date: Fri, 1 Mar 2024 09:07:31 -0600 Subject: [PATCH 57/58] Update README.MD --- src/asf_tools/watermasking/README.MD | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/asf_tools/watermasking/README.MD b/src/asf_tools/watermasking/README.MD index 8ef49a3c..dde40451 100644 --- a/src/asf_tools/watermasking/README.MD +++ b/src/asf_tools/watermasking/README.MD @@ -2,16 +2,17 @@ These scripts are for creating a global (or regional) water mask dataset based o For the OSM water mask dataset, follow these steps to replicate our dataset: - 1. Download the "Latest Weekly Planet PBF File" file from here: "https://planet.openstreetmap.org/". - 2. Download the WGS84 water polygons shapefile from: "https://osmdata.openstreetmap.de/data/water-polygons.html". - 3. The files should be unzipped and you should have something like `planet.osm.pbf` or `planet.pbf` and `water_polygons.shp` (and the support files for `water_polygons.shp`). - 4. Run ```generate_osm_dataset --planet-file-path [path-to-planet.pbf] --ocean-polygons-path [path-to-water-polygons.shp] --lat-begin -85 --lat-end 85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` - 5. Run ```fill_missing_tiles --fill-value 0 --lat-begin -90 --lat-end -85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` - 6. Run ```fill_missing_tiles --fill-value 1 --lat-begin 85 --lat-end 90 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` +1. Download the "Latest Weekly Planet PBF File" file from here: "https://planet.openstreetmap.org/". +2. Download the WGS84 water polygons shapefile from: "https://osmdata.openstreetmap.de/data/water-polygons.html". +3. The files should be unzipped and you should have something like `planet.osm.pbf` or `planet.pbf` and `water_polygons.shp` (and the support files for `water_polygons.shp`). +4. Run ```generate_osm_dataset --planet-file-path [path-to-planet.pbf] --ocean-polygons-path [path-to-water-polygons.shp] --lat-begin -85 --lat-end 85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` +5. Run ```fill_missing_tiles --fill-value 0 --lat-begin -90 --lat-end -85 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` +6. Run ```fill_missing_tiles --fill-value 1 --lat-begin 85 --lat-end 90 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` For the WorldCover water mask dataset, follow these steps: - 1. Download the portions of the dataset for the areas you would like to cover from here: "https://worldcover2020.esa.int/downloader" - 2. Extract the contents into a folder. Note, if you download multiple portions of the dataset, extract them all into the same folder. - 3. Run ```generate_worldcover_dataset --worldcover-tiles-dir [path-to-worldcover-data] --lat-begin 55 --lat-end 80 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` - Note that we only use WorldCover data over Alaska, Canada, and Russia for our dataset. +1. Download the portions of the dataset for the areas you would like to cover from here: "https://worldcover2020.esa.int/downloader" +2. Extract the contents into a folder. Note, if you download multiple portions of the dataset, extract them all into the same folder. +3. Run ```generate_worldcover_dataset --worldcover-tiles-dir [path-to-worldcover-data] --lat-begin 55 --lat-end 80 --lon-begin -180 --lon-end 180 --tile-width 5 --tile-height 5``` + +Note that we only use WorldCover data over Alaska, Canada, and Russia for our dataset. From 12d959bec872677f062faf59884a841fd43673a5 Mon Sep 17 00:00:00 2001 From: Andrew Player Date: Fri, 1 Mar 2024 10:20:53 -0600 Subject: [PATCH 58/58] attriberror to keyerror and typo fix --- src/asf_tools/watermasking/generate_osm_tiles.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index d9d391b6..4b934900 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -64,7 +64,7 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig gdal.Rasterize( tile_tif, - ocean_polygons_path, + clipped_polygons_path, xRes=pixel_size_x, yRes=pixel_size_y, burnValues=1, @@ -106,7 +106,7 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio try: water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'island'].index) water_gdf = water_gdf.drop(water_gdf[water_gdf['place'] == 'islet'].index) - except AttributeError: + except KeyError: # When there are no islands to remove, an AttributeError should throw, but we don't care about it. pass water_gdf.to_file(tile_shp, mode='w', engine='pyogrio') @@ -140,7 +140,7 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di merged_tile_dir: The path to the directory containing the merged water tiles. """ index = 0 - num_tiles = len([f for f in os.listdir(internal_tile_dir) if f.endswith('tif')]) + num_tiles = len([f for f in os.listdir(internal_tile_dir) if f.endswith('tif')]) - 1 for filename in os.listdir(internal_tile_dir): if filename.endswith('.tif'): start_time = time.time() @@ -218,7 +218,7 @@ def main(): print('Processing tiles...') lat_range = range(lat_begin, lat_end, tile_height) lon_range = range(lon_begin, lon_end, tile_width) - num_tiles = len(lat_range) * len(lon_range) + num_tiles = len(lat_range) * len(lon_range) - 1 index = 0 for lat in lat_range: for lon in lon_range: