Skip to content

Added MZ example notebooks #62

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Added MZ example notebooks #62

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
249 changes: 249 additions & 0 deletions notebooks/example-workflows/part1_radar_image_processing.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,249 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "8e1316f1-3245-48e3-894d-73c906430829",
"metadata": {},
"source": [
"## Image Processing with Radar Reflectivity\n",
"This program accepts radar scans plotted onto a cartesian grid. The program then thresholds the image to perform image processing and analysis on the largest echo of reflectivity within the image, and outputs a CSV file with the parameters that were calculated for each radar file. </br>"
]
},
{
"cell_type": "markdown",
"id": "20d8b66a-458c-4b3f-9653-c65191fc4262",
"metadata": {},
"source": [
"#### Author\n",
"<b>Maggie Zoerner</b></br>\n",
"Argonne National Laboratory, <i>Environmental Science Division</i></br>\n",
"Student Undergraduate Laboratory Intern </br>"
]
},
{
"cell_type": "markdown",
"id": "f0f532b2-9a88-4fc5-a852-dfd37f3e5225",
"metadata": {},
"source": [
"#### Import the libraries."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "47afec8d-2dfa-4503-ad82-cd9212719502",
"metadata": {},
"outputs": [],
"source": [
"from netCDF4 import Dataset\n",
"from scipy import ndimage\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import cv2\n",
"import csv\n",
"from glob import glob"
]
},
{
"cell_type": "markdown",
"id": "f32221d3-af43-4bd6-9a18-23cc4585f5c3",
"metadata": {},
"source": [
"#### Import the files."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "716d89f4-6f90-48a8-8d6b-b95bbbc7d69d",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"files = glob('season_grids/KILX2020*')"
]
},
{
"cell_type": "markdown",
"id": "12f7b207-62e5-4eeb-8ee8-d625d1d7e2ab",
"metadata": {},
"source": [
"#### In the below step: </br>\n",
"\n",
"- Go through the data set.\n",
"- Remove reflectivities smaller than 35 dBz to remove non-convective precip.\n",
"- Clear smaller echoes to remove noise.\n",
"- Calculate the raw, central, normalized central, and invariant moments.\n",
"- Perform shape computations based on the moments.\n",
"- Find the largest external contour of the echo.\n",
"- Perform shape computations based on the contours.\n",
"- Export computations to CSV file."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "63ae9194-7cb8-4b9d-9022-ceacc2280bd6",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# LOOP THROUGH THE DATASET\n",
"\n",
"header_added = False\n",
" \n",
"for file in files:\n",
" ncfile = Dataset(file, \"r\")\n",
" ref = ncfile.variables['reflectivity'][:]\n",
" ref = ref.squeeze()\n",
" ref_thre = np.where(ref < 35, 0, ref)\n",
" labeled_echo = ndimage.label(ref_thre)[0]\n",
"\n",
" # Remove smaller echoes \n",
" \n",
" def clear_small_echoes(label_image, min_size):\n",
" \"\"\" Takes in binary image and clears objects less than min_size. \"\"\"\n",
" flat_image = pd.Series(label_image.flatten())\n",
" flat_image = flat_image[flat_image > 0]\n",
" size_table = flat_image.value_counts(sort=False)\n",
" small_objects = size_table.keys()[size_table < min_size]\n",
"\n",
" for obj in small_objects:\n",
" label_image[label_image == obj] = 0\n",
" label_image = ndimage.label(label_image)\n",
" return label_image[0]\n",
" \n",
" max_echo = clear_small_echoes(labeled_echo, 175)\n",
" \n",
" # Convert to float\n",
"\n",
" max_echo = max_echo.astype('float32')\n",
" \n",
" # Calculate geometric moments and moment invariants\n",
" \n",
" moments =cv2.moments(max_echo)\n",
"\n",
" Hu_moments = cv2.HuMoments(moments)\n",
" \n",
" # The 0th moment gives the area\n",
" \n",
" M = moments\n",
" area = M['m00']\n",
" \n",
" # Ensure that the area is not equal to zero (as this will mess up future calculations)\n",
" \n",
" try:\n",
" \n",
" # Calculate the center of mass\n",
" Cx = int(M['m10'] / M['m00'])\n",
" Cy = int(M['m01'] / M['m00'])\n",
" \n",
" # The 3rd normalized central moment gives the skewness (the deviation of the respective projection from symmetry)\n",
" \n",
" skewness_x = (M['nu30'])\n",
" skewness_y = (M['nu03'])\n",
"\n",
" # Find the deviation from circlular shape\n",
" \n",
" dev_circ = abs(M['nu11'])\n",
" \n",
" # Convert to a format accepted by OpenCV contour analysis\n",
" \n",
" max_echo_u8 = max_echo.astype(np.uint8)\n",
" \n",
" # Find and draw the contours in the image\n",
" \n",
" u8_copy = max_echo_u8.copy()\n",
" contours, hierarchy = cv2.findContours(u8_copy, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n",
" \n",
" cnt = contours[0]\n",
" contour_image = cv2.drawContours(u8_copy, cnt, -1, (0,255,0),3)\n",
"\n",
" # Find the perimeter of the shape\n",
"\n",
" perimeter = cv2.arcLength(cnt,True)\n",
" \n",
" # Find the elongation (aspect ratio) of the shape \n",
" \n",
" x,y,w,h = cv2.boundingRect(cnt)\n",
" \n",
" elongation = float(w)/h\n",
"\n",
" # Find the equivalent diameter: diameter of a circle whose diameter is the same as the contour area\n",
" \n",
" area = cv2.contourArea(cnt)\n",
" equi_diameter = np.sqrt(4*area/np.pi)\n",
" \n",
" # Find the roundness \n",
" \n",
" roundness = (4*np.pi*(area))/((perimeter)**2)\n",
" \n",
" # Find the compactness \n",
" \n",
" compactness = ((perimeter)**2)/(4*np.pi*(area))\n",
" \n",
" # Find the eccentricity \n",
" \n",
" (x,y),(MA,ma),angle = cv2.fitEllipse(cnt)\n",
" ecc = ma/MA\n",
" \n",
"\n",
" # WRITING TO A CSV FILE\n",
" \n",
" # Define header names\n",
" \n",
" headers = ['File Name', 'Area', 'Perimeter', 'Center of Mass (X)', 'Center of Mass (Y)', 'Skewness (X)', 'Skewness (Y)', \n",
" 'Elongation','Equivalent Diameter','Compactness', 'Roundness', 'Eccentricity', 'Deviation from Circle']\n",
" \n",
" # Define data to include\n",
" \n",
" rows = [(file), (area), (perimeter), (Cx), (Cy), (skewness_x), (skewness_y), (elongation), (equi_diameter), \n",
" (compactness), (roundness), (ecc), (dev_circ)]\n",
" \n",
" # Name the CSV file\n",
" \n",
" filename = \"parameter_computations.csv\"\n",
" \n",
" # Ensure header does not duplicate in loop and populate file\n",
" \n",
" with open(filename, 'a') as csvfile:\n",
" csvwriter = csv.writer(csvfile)\n",
" if not header_added:\n",
" csvwriter.writerow(headers)\n",
" header_added = True\n",
" csvwriter.writerow(rows)\n",
" \n",
" except ZeroDivisionError:\n",
" pass"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.12"
},
"toc-autonumbering": true,
"toc-showcode": true,
"toc-showmarkdowntxt": false,
"toc-showtags": false
},
"nbformat": 4,
"nbformat_minor": 5
}
Loading