diff --git a/notebook/crop_image.ipynb b/notebook/crop_image.ipynb index 92ce7e7..c676722 100644 --- a/notebook/crop_image.ipynb +++ b/notebook/crop_image.ipynb @@ -23,20 +23,30 @@ "execution_count": 1, "id": "fbcc9b0b", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/Users/krismannino/Code/CADS/IMPACT/uavsar-wildfire-main-jpl/uavsar-wildfire/python\n" + ] + } + ], "source": [ "import sys\n", "from pathlib import Path\n", "\n", "# Add the path to the utils folder to sys.path\n", "utils_path = Path('../python').resolve()\n", + "print(utils_path)\n", "sys.path.append(str(utils_path))\n", "\n", "from pathlib import Path\n", "from rasterio.crs import CRS\n", "from crop_utils import (crop_image_by_coordinates, \n", " crop_image_by_geojson_shp, \n", - " reproject_geotiff)" + " reproject_geotiff)\n", + "from edit_path_utils import (edit_paths)" ] }, { @@ -45,15 +55,17 @@ "metadata": {}, "source": [ "---\n", - "## Crop image by radius from center\n", - "Returns a cropped version of the input image based intersection of the specificed circle\n", + "## Crop image by box from center\n", + "Returns a cropped version of the input image based on rectangle centered on coordinate\n", "\n", "**Parameters**:\n", "- `path_to_images` (list): a list containing the paths to the images cropping from \n", "- `output_names` (list): a list containing the output names for the cropped images [**.tif** files]\n", "- `center_longitude` (float): center longitude of the fire\n", "- `center_latitude` (float): center latitude of the fire\n", - "- `radius_in_km` (float): radius from the coordinate in km to crop" + "- `width_km` (float): full width of box in km to crop centered on coordinate \n", + "- `height_km` (float): full height of box in km to crop centered on coordinate\n", + "" ] }, { @@ -63,22 +75,29 @@ "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "[PosixPath('/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_HHHH_rtc.grd'),\n", - " PosixPath('/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_HVHV_rtc.grd'),\n", - " PosixPath('/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_VVVV_rtc.grd')]" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_L090_CX_01/SanAnd_08525_18076_003_181011_HVHV_rtc.grd\n", + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_L090_CX_01/SanAnd_08525_21041_007_210526_HVHV_rtc.grd\n" + ] } ], "source": [ - "data_dir = Path('/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01')\n", + "data_dir = Path('/Volumes/BlackT7/bobcat_flight_paths')\n", "tifs = sorted(list(data_dir.rglob('./*_rtc.grd')))\n", - "tifs" + "\n", + "# adjust for wanted and unwanted tifs\n", + "include = [\n", + " \"181011\",\n", + " \"210526\"\n", + " ]\n", + "exclude = [\n", + " \"26526\"\n", + " ]\n", + "\n", + "tifs = edit_paths(include,exclude,tifs)\n", + "print(*tifs, sep=\"\\n\")" ] }, { @@ -91,17 +110,23 @@ "name": "stdout", "output_type": "stream", "text": [ - "['/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_HHHH_rtc_cropped_5km.tif', '/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_HVHV_rtc_cropped_5km.tif', '/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_VVVV_rtc_cropped_5km.tif']\n" + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_HVHV_rtc_box_crop_45km_40km.tif\n", + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_HVHV_rtc_box_crop_45km_40km.tif\n" ] } ], "source": [ "path_to_images = tifs\n", - "center_longitude = -112.050134\n", - "center_latitude = 38.478988\n", - "radius_in_km = 5\n", - "output_names = [str(data_dir) + '/' + file.stem + '_cropped_' + str(radius_in_km) + 'km.tif' for file in path_to_images]\n", - "print(output_names)##" + "\n", + "# Enter center longitude and latitude for Area of Interest\n", + "center_longitude = -117.92976\n", + "center_latitude = 34.33253\n", + "# radius_in_km = 25 #radius functionality in python/crop_utlis\n", + "width_km = 45 ## total width and height\n", + "height_km =40\n", + "\n", + "output_names = [str(data_dir) + '/' + file.stem + '_box_crop_' + str(width_km) +\"km_\"+ str(height_km) + 'km.tif' for file in path_to_images]\n", + "print(*output_names, sep=\"\\n\")" ] }, { @@ -117,15 +142,16 @@ "cell_type": "code", "execution_count": 4, "id": "15972a4d-6427-45b8-ad51-91cb533d7a7a", - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_HHHH_rtc_cropped_5km.tif is outputted!\n", - "/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_HVHV_rtc_cropped_5km.tif is outputted!\n", - "/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_VVVV_rtc_cropped_5km.tif is outputted!\n" + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_HVHV_rtc_box_crop_45km_40km.tif is outputted!\n", + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_HVHV_rtc_box_crop_45km_40km.tif is outputted!\n" ] } ], @@ -135,7 +161,8 @@ " output_names[i], \n", " center_longitude, \n", " center_latitude, \n", - " radius_in_km)" + " width_km,\n", + " height_km)" ] }, { @@ -145,20 +172,24 @@ "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "[PosixPath('/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_L090_CX_01.inc')]" - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_L090_CX_01/SanAnd_08525_18076_003_181011_L090_CX_01.inc\n", + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_L090_CX_01/SanAnd_08525_21041_007_210526_L090_CX_01.inc\n" + ] } ], "source": [ "##data_dir = Path('../data/latuna')\n", "incs = sorted(list(data_dir.rglob('./*.inc')))\n", - "incs" + "# include = [\n", + "# ]\n", + "# exclude = [\n", + "# ]\n", + "\n", + "incs = edit_paths(include,exclude,incs)\n", + "print(*incs, sep=\"\\n\")" ] }, { @@ -171,17 +202,15 @@ "name": "stdout", "output_type": "stream", "text": [ - "['/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_L090_CX_01_inc_5km.tif']\n" + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_L090_CX_01_inc_box_crop_45km_40km.tif\n", + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_L090_CX_01_inc_box_crop_45km_40km.tif\n" ] } ], "source": [ "path_to_images = incs\n", - "##center_longitude = -118.2674\n", - "##center_latitude = 34.22957\n", - "##radius_in_km = 10\n", - "output_names = [str(data_dir) + '/' + file.stem + '_inc_' + str(radius_in_km) + 'km.tif' for file in path_to_images]\n", - "print(output_names)##" + "output_names = [str(data_dir) + '/' + file.stem + '_inc_box_crop_' + str(width_km) +\"km_\"+ str(height_km) + 'km.tif' for file in path_to_images]\n", + "print(*output_names, sep=\"\\n\")" ] }, { @@ -197,13 +226,16 @@ "cell_type": "code", "execution_count": 7, "id": "167de085", - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "/mnt/karenan/fasmee/fishlake/fishla_19601_23024_001_231017_L090_CX_01/fishla_19601_23024_001_231017_L090_CX_01_inc_5km.tif is outputted!\n" + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_L090_CX_01_inc_box_crop_45km_40km.tif is outputted!\n", + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_L090_CX_01_inc_box_crop_45km_40km.tif is outputted!\n" ] } ], @@ -213,7 +245,8 @@ " output_names[i], \n", " center_longitude, \n", " center_latitude, \n", - " radius_in_km)" + " width_km,\n", + " height_km)" ] }, { @@ -235,96 +268,114 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 8, "id": "12698266", "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "[PosixPath('/mnt/karenan/fasmee/fasmee_fall_2023_burn_perimeter/fasmee_fall_2023_burn_perimeter.shp')]" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "/Volumes/BlackT7/bobcat_perimeter_wen/._bobcat_perimeter_wen.shp\n", + "/Volumes/BlackT7/bobcat_perimeter_wen/bobcat_perimeter_wen.shp\n" + ] } ], "source": [ - "shp_dir = Path('/mnt/karenan/fasmee/fasmee_fall_2023_burn_perimeter')\n", + "# OPTIONAL\n", + "# Regressive search directory here \n", + "folder_path = '/Volumes/BlackT7/bobcat_perimeter_wen' # path to shapefiles\n", + "\n", + "shp_dir = Path(folder_path)\n", "shps = sorted(list(shp_dir.glob('./*.shp')))\n", - "shps" + "\n", + "# OPTIONAL\n", + "# adjust for wanted and unwanted shps\n", + "# include = [ \n", + "# ]\n", + "# exclude = [\n", + "# ]\n", + "# shps = edit_paths(include,exclude,shps)\n", + "\n", + "print(*shps, sep=\"\\n\")" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 9, "id": "21463bff-be6f-4545-8498-5591405ab11f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "'/mnt/karenan/fasmee/fasmee_fall_2023_burn_perimeter/fasmee_fall_2023_burn_perimeter.shp'" + "'/Volumes/BlueT7/bobcat_burn_perim/mtbs/2020/ca3424811795920200906/ca3424811795920200906_20200703_20210706_burn_bndy.shp'" ] }, - "execution_count": 12, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "##shps = '../data/bobcat/uavsar_perimeter/bobcat_perimeter_bilinear_inc_south.geojson'\n", - "shps = '/mnt/karenan/fasmee/fasmee_fall_2023_burn_perimeter/fasmee_fall_2023_burn_perimeter.shp'\n", - "shps" + "# Direct enter (NOT OPTIONAL)\n", + "# from path or list above\n", + "shp = '/Volumes/BlueT7/bobcat_burn_perim/mtbs/2020/ca3424811795920200906/ca3424811795920200906_20200703_20210706_burn_bndy.shp'\n", + "shp" ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 10, "id": "5a1260af", "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "[PosixPath('/mnt/karenan/fasmee/fishlake/classify_fishla_01601_19601_230727_231017/fishlake_weighted_inc_merge_hv_0.tif'),\n", - " PosixPath('/mnt/karenan/fasmee/fishlake/classify_fishla_01601_19601_230727_231017/fishlake_weighted_inc_merge_hv_1.tif')]" - ] - }, - "execution_count": 13, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_L090_CX_01/SanAnd_08525_18076_003_181011_HVHV_rtc.grd\n", + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_L090_CX_01/SanAnd_08525_21041_007_210526_HVHV_rtc.grd\n" + ] } ], "source": [ - "data_dir = Path('/mnt/karenan/fasmee/fishlake/classify_fishla_01601_19601_230727_231017')\n", - "tifs = sorted(list(data_dir.glob('./*weighted*.tif')))\n", - "tifs" + "data_dir = Path('/Volumes/BlackT7/bobcat_flight_paths')\n", + "# tifs = sorted(list(data_dir.glob('./*weighted*.tif')))\n", + "# tifs = sorted(list(data_dir.rglob('./*_rtc.grd')))\n", + "\n", + "# # adjust for wanted and unwanted tifs\n", + "# include = [\n", + "# ]\n", + "# exclude = [\n", + "# ]\n", + "\n", + "tifs = edit_paths(include,exclude,tifs)\n", + "print(*tifs, sep=\"\\n\")" ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 11, "id": "5d7e63f9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "['fishlake_weighted_inc_merge_hv_0_perimeter_intersection_uavsar.tif',\n", - " 'fishlake_weighted_inc_merge_hv_1_perimeter_intersection_uavsar.tif']" + "['/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_HVHV_rtc_perimeter_crop.tif',\n", + " '/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_HVHV_rtc_perimeter_crop.tif']" ] }, - "execution_count": 14, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "path_to_polygon_file = shps\n", + "path_to_polygon_file = shp\n", "path_to_images = tifs\n", - "output_names = [file.stem + '_perimeter_intersection_uavsar.tif' for file in tifs]\n", + "output_names = [str(data_dir) + '/' + file.stem + '_perimeter_crop.tif' for file in tifs]\n", "output_names" ] }, @@ -341,16 +392,18 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 12, "id": "910f5202", - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "fishlake_weighted_inc_merge_hv_0_perimeter_intersection_uavsar.tif is outputted!\n", - "fishlake_weighted_inc_merge_hv_1_perimeter_intersection_uavsar.tif is outputted!\n" + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_18076_003_181011_HVHV_rtc_perimeter_crop.tif is outputted.\n", + "/Volumes/BlackT7/bobcat_flight_paths/SanAnd_08525_21041_007_210526_HVHV_rtc_perimeter_crop.tif is outputted.\n" ] } ], @@ -380,28 +433,48 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "id": "6c4d9328", - "metadata": {}, - "outputs": [], + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/Volumes/BlueT7/test_bob/SanAnd_08525_18076_003_181011_HVHV_rtc_cropped_25km.tif\n" + ] + } + ], "source": [ - "data_dir = Path('../data/MTBS_sbs')\n", + "data_dir = Path('/Volumes/BlueT7/test_bob/')\n", "tifs = sorted(list(data_dir.glob('./*.tif')))\n", - "tifs" + "\n", + "# adjust for wanted and unwanted tifs\n", + "include = [\n", + " \"181011\"\n", + " ]\n", + "exclude = [\n", + " \"inc\"\n", + " ]\n", + "\n", + "tifs = edit_paths(include,exclude,tifs)\n", + "print(*tifs, sep=\"\\n\")" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "id": "2be64663", "metadata": { "scrolled": true }, "outputs": [], "source": [ - "target_crs = CRS.from_epsg(4326) \n", + "target_crs = CRS.from_epsg(26911) \n", "input_raster_path = tifs[0]\n", - "output_raster_path = 'mtbs_CA_2009_crs4326.tif'" + "output_raster_path = 'SanAnd_08525_18076_003_181011_HVHV_rtc_cropped_25km_crs4326.tif'" ] }, { @@ -415,10 +488,18 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "id": "f92a5dd1", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reprojected raster saved to: SanAnd_08525_18076_003_181011_HVHV_rtc_cropped_25km_crs4326.tif\n" + ] + } + ], "source": [ "reproject_geotiff(target_crs, input_raster_path, output_raster_path)" ] @@ -426,9 +507,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "classification", "language": "python", - "name": "python3" + "name": "classification" }, "language_info": { "codemirror_mode": { @@ -440,7 +521,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.2" + "version": "3.12.7" } }, "nbformat": 4, diff --git a/python/crop_utils.py b/python/crop_utils.py index 0510a74..2dfdc65 100644 --- a/python/crop_utils.py +++ b/python/crop_utils.py @@ -1,9 +1,9 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ -@author: Wen Tao Lin +@author: Wen Tao Lin, Kris Mannino """ -from shapely.geometry import Point +from shapely.geometry import Point, box from shapely.geometry import mapping from rasterio.mask import mask from rasterio.crs import CRS @@ -12,6 +12,22 @@ import numpy as np import geopandas as gpd +def generate_rectangle_area(longitude: float, + latitude: float, + width_km: float, + height_km: float) -> box: + width_deg = width_km / 111.13 + height_deg = height_km / (111.13 * np.cos(np.radians(latitude))) + + # Calculate bounds + min_lon = longitude - width_deg/2 + max_lon = longitude + width_deg/2 + min_lat = latitude - height_deg/2 + max_lat = latitude + height_deg/2 + + return box(min_lon, min_lat, max_lon, max_lat) + +# Radius functionality disabled def generate_area(longitude: float, latitude: float, radius_in_km: float) -> Point: @@ -23,7 +39,10 @@ def crop_image_by_coordinates(path_to_image: str, output_name: str, center_longitude: float, center_latitude: float, - radius_in_km: float) -> None: + # radius_in_km: float + width_km: float, + height_km: float + ) -> None: """ Crops a raster image by finding the intersection of the original image and the circle generated by the user inputted radius and center coordinates. @@ -48,8 +67,12 @@ def crop_image_by_coordinates(path_to_image: str, based intersection of the specificed circle. """ + # Generate the rectangular polygon + aoi = generate_rectangle_area(center_longitude, center_latitude, width_km, height_km) + + # Radius functionality disabled # Generate the circular polygon - aoi = generate_area(center_longitude, center_latitude, radius_in_km) + # aoi = generate_area(center_longitude, center_latitude, radius_in_km) shapes = [mapping(aoi)] # Crop the raster using the polygon diff --git a/python/edit_path_utils.py b/python/edit_path_utils.py new file mode 100644 index 0000000..9167168 --- /dev/null +++ b/python/edit_path_utils.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +@author: Kris Mannino +""" +from pathlib import Path +def edit_paths(include: list[str], + exclude: list[str], + paths: list[Path]) -> list[Path]: + """ + Crops a raster image based on the shape of the first polygon in the + inputted GeoJSON or shapefile. + + Parameters + ---------- + include : list[str] + List of strings that must exist in Path. + exclude : list[str] + List of strings that if exists in Path, Path is removed. + paths: list[Path] + List of initial Paths. + + Returns + ------- + list[Path] + Returns adjust Path list with inclusions and exclusions. + + Notes + ------- + + """ + if include: + filtered_paths = [ + p for p in paths + if any(term in str(p) for term in include) + ] + else: + filtered_paths = paths + + if exclude: + filtered_paths = [ + p for p in filtered_paths + if not any(term in str(p) for term in exclude) + ] + + return filtered_paths