Skip to content

Commit 2034c17

Browse files
JanushiShastri2001JanushiShastri2001
committed
Added topogrphic correction code for NEON (NIWO) and Drone data
1 parent f48c52d commit 2034c17

7 files changed

+4667
-0
lines changed

‎.gitignore

Whitespace-only changes.

‎Calulating_Sun_angles.ipynb

+277
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,277 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": 2,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import rasterio\n",
10+
"import numpy as np\n",
11+
"from pyproj import Proj, transform\n",
12+
"from rasterio.transform import from_origin\n",
13+
"import ephem\n",
14+
"import datetime\n",
15+
"import math"
16+
]
17+
},
18+
{
19+
"cell_type": "markdown",
20+
"metadata": {},
21+
"source": [
22+
"### Function to get the lat long from the drone image\n",
23+
"### Load drone data path in the function"
24+
]
25+
},
26+
{
27+
"cell_type": "code",
28+
"execution_count": 3,
29+
"metadata": {},
30+
"outputs": [
31+
{
32+
"name": "stderr",
33+
"output_type": "stream",
34+
"text": [
35+
"C:\\Users\\JANUSHI SHASTRI\\AppData\\Local\\Temp\\ipykernel_48828\\3174482380.py:13: DeprecationWarning: Right multiplication will be prohibited in version 3.0\n",
36+
" rc2en = lambda r, c: (c, r) * T1\n",
37+
"C:\\Users\\JANUSHI SHASTRI\\AppData\\Local\\Temp\\ipykernel_48828\\3174482380.py:17: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1\n",
38+
" longs, lats = transform(p1, p2, eastings, northings)\n"
39+
]
40+
},
41+
{
42+
"name": "stdout",
43+
"output_type": "stream",
44+
"text": [
45+
"Longitudes: [[-105.55816094 -105.55816011 -105.55815927 ... -105.55651632\n",
46+
" -105.55651549 -105.55651466]\n",
47+
" [-105.55816094 -105.55816011 -105.55815928 ... -105.55651633\n",
48+
" -105.5565155 -105.55651466]\n",
49+
" [-105.55816095 -105.55816012 -105.55815928 ... -105.55651633\n",
50+
" -105.5565155 -105.55651467]\n",
51+
" ...\n",
52+
" [-105.55817339 -105.55817256 -105.55817173 ... -105.55652874\n",
53+
" -105.55652791 -105.55652707]\n",
54+
" [-105.5581734 -105.55817256 -105.55817173 ... -105.55652875\n",
55+
" -105.55652791 -105.55652708]\n",
56+
" [-105.5581734 -105.55817257 -105.55817174 ... -105.55652875\n",
57+
" -105.55652792 -105.55652709]]\n",
58+
"Latitudes: [[40.03941242 40.03941242 40.03941243 ... 40.03942033 40.03942033\n",
59+
" 40.03942034]\n",
60+
" [40.03941306 40.03941306 40.03941307 ... 40.03942097 40.03942097\n",
61+
" 40.03942098]\n",
62+
" [40.0394137 40.0394137 40.03941371 ... 40.03942161 40.03942161\n",
63+
" 40.03942162]\n",
64+
" ...\n",
65+
" [40.04093971 40.04093971 40.04093971 ... 40.04094762 40.04094762\n",
66+
" 40.04094763]\n",
67+
" [40.04094035 40.04094035 40.04094035 ... 40.04094826 40.04094826\n",
68+
" 40.04094827]\n",
69+
" [40.04094099 40.04094099 40.04094099 ... 40.0409489 40.0409489\n",
70+
" 40.04094891]]\n"
71+
]
72+
}
73+
],
74+
"source": [
75+
"def get_coordinates(file_path):\n",
76+
" # Load geotiff\n",
77+
" with rasterio.open(file_path) as r:\n",
78+
" T0 = r.transform # upper-left pixel corner affine transform\n",
79+
" p1 = Proj(r.crs)\n",
80+
" A = r.read() # pixel values\n",
81+
"\n",
82+
" # All rows and columns\n",
83+
" cols, rows = np.meshgrid(np.arange(A.shape[2]), np.arange(A.shape[1]))\n",
84+
"\n",
85+
" # Get pixel coordinates from the rows and columns\n",
86+
" T1 = T0 * rasterio.Affine.scale(1, -1)\n",
87+
" rc2en = lambda r, c: (c, r) * T1 \n",
88+
" eastings, northings = np.vectorize(rc2en, otypes=[float, float])(rows, cols)\n",
89+
"\n",
90+
" p2 = Proj(proj='latlong', datum='WGS84')\n",
91+
" longs, lats = transform(p1, p2, eastings, northings)\n",
92+
"\n",
93+
" return longs, lats\n",
94+
"\n",
95+
"file_path = r'C:\\Users\\JANUSHI SHASTRI\\Desktop\\CIRES_PROJECT\\Drone Data\\drive-download-20230711T205208Z-001\\niwot_6_23_2022_ortho_cropped.tif'\n",
96+
"longitudes, latitudes = get_coordinates(file_path)\n",
97+
"\n",
98+
"print('Longitudes:', longitudes)\n",
99+
"print('Latitudes:', latitudes)"
100+
]
101+
},
102+
{
103+
"cell_type": "markdown",
104+
"metadata": {},
105+
"source": [
106+
"### Method 1:"
107+
]
108+
},
109+
{
110+
"cell_type": "code",
111+
"execution_count": 8,
112+
"metadata": {},
113+
"outputs": [
114+
{
115+
"name": "stdout",
116+
"output_type": "stream",
117+
"text": [
118+
"Sun azimuth in degrees: 84.64035306276917, Sun zenith in degrees: 58.434557373447454\n"
119+
]
120+
}
121+
],
122+
"source": [
123+
"# Calculate the central latitude and longitude\n",
124+
"center_long = np.mean(longitudes)\n",
125+
"center_lat = np.mean(latitudes)\n",
126+
"\n",
127+
"# Convert the date and time to UTC. The time given is 2:34 PM, which is 14:34 in 24-hour format\n",
128+
"date_time_str = '2022-06-23 14:34:00'\n",
129+
"\n",
130+
"# Calculation of sun azimuth and zenith angle\n",
131+
"# Ephemeris computations\n",
132+
"observer = ephem.Observer()\n",
133+
"observer.lat = str(center_lat)\n",
134+
"observer.lon = str(center_long)\n",
135+
"observer.date = ephem.Date(date_time_str)\n",
136+
"\n",
137+
"sun = ephem.Sun()\n",
138+
"sun.compute(observer)\n",
139+
"\n",
140+
"azimuth = sun.az\n",
141+
"zenith = ephem.degrees('90') - sun.alt\n",
142+
"azimuth_deg = math.degrees(azimuth)\n",
143+
"zenith_deg = math.degrees(zenith)\n",
144+
"\n",
145+
"print(f'Sun azimuth in degrees: {azimuth_deg}, Sun zenith in degrees: {zenith_deg}')\n",
146+
"# print(f'Sun azimuth: {azimuth}, Sun zenith: {zenith}')"
147+
]
148+
},
149+
{
150+
"cell_type": "code",
151+
"execution_count": null,
152+
"metadata": {},
153+
"outputs": [],
154+
"source": []
155+
},
156+
{
157+
"cell_type": "markdown",
158+
"metadata": {},
159+
"source": [
160+
"### Method 2"
161+
]
162+
},
163+
{
164+
"cell_type": "markdown",
165+
"metadata": {},
166+
"source": [
167+
"### Function to calculate sun angles using all pixels of image and lat long"
168+
]
169+
},
170+
{
171+
"cell_type": "code",
172+
"execution_count": 4,
173+
"metadata": {},
174+
"outputs": [
175+
{
176+
"name": "stderr",
177+
"output_type": "stream",
178+
"text": [
179+
"C:\\Users\\JANUSHI SHASTRI\\AppData\\Local\\Temp\\ipykernel_46172\\550025637.py:13: DeprecationWarning: Right multiplication will be prohibited in version 3.0\n",
180+
" rc2en = lambda r, c: (c, r) * T1\n",
181+
"C:\\Users\\JANUSHI SHASTRI\\AppData\\Local\\Temp\\ipykernel_46172\\550025637.py:17: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1\n",
182+
" longs, lats = transform(p1, p2, eastings, northings)\n"
183+
]
184+
},
185+
{
186+
"name": "stdout",
187+
"output_type": "stream",
188+
"text": [
189+
"Sun azimuth in degrees: 84.64035306276917, Sun zenith in degrees: 58.434557373447454\n"
190+
]
191+
}
192+
],
193+
"source": [
194+
"def get_sun_angles(file_path, date_time_str):\n",
195+
" # Load geotiff\n",
196+
" with rasterio.open(file_path) as r:\n",
197+
" T0 = r.transform # upper-left pixel corner affine transform\n",
198+
" p1 = Proj(r.crs)\n",
199+
" A = r.read() # pixel values\n",
200+
"\n",
201+
" # All rows and columns\n",
202+
" cols, rows = np.meshgrid(np.arange(A.shape[2]), np.arange(A.shape[1]))\n",
203+
"\n",
204+
" # Get pixel coordinates from the rows and columns\n",
205+
" T1 = T0 * rasterio.Affine.scale(1, -1)\n",
206+
" rc2en = lambda r, c: (c, r) * T1 \n",
207+
" eastings, northings = np.vectorize(rc2en, otypes=[float, float])(rows, cols)\n",
208+
"\n",
209+
" p2 = Proj(proj='latlong',datum='WGS84')\n",
210+
" longs, lats = transform(p1, p2, eastings, northings)\n",
211+
"\n",
212+
" # Calculation of sun azimuth and zenith angle\n",
213+
" # We'll do this for the center of the image\n",
214+
" center_lat = np.mean(lats)\n",
215+
" center_long = np.mean(longs)\n",
216+
"\n",
217+
" # Convert the date and time to UTC.\n",
218+
" date_time_obj = datetime.datetime.strptime(date_time_str, '%Y-%m-%d %H:%M:%S')\n",
219+
"\n",
220+
" # Ephemeris computations\n",
221+
" observer = ephem.Observer()\n",
222+
" observer.lat = str(center_lat)\n",
223+
" observer.lon = str(center_long)\n",
224+
" observer.date = ephem.Date(date_time_obj)\n",
225+
"\n",
226+
" sun = ephem.Sun()\n",
227+
" sun.compute(observer)\n",
228+
"\n",
229+
" azimuth = sun.az\n",
230+
" zenith = ephem.degrees('90') - sun.alt\n",
231+
" # Convert radians to degrees\n",
232+
" azimuth_deg = math.degrees(azimuth)\n",
233+
" zenith_deg = math.degrees(zenith)\n",
234+
"\n",
235+
" print(f'Sun azimuth in degrees: {azimuth_deg}, Sun zenith in degrees: {zenith_deg}')\n",
236+
"\n",
237+
" return azimuth, zenith\n",
238+
"\n",
239+
"file_path = r'C:\\Users\\JANUSHI SHASTRI\\Desktop\\CIRES_PROJECT\\Drone Data\\drive-download-20230711T205208Z-001\\niwot_6_23_2022_ortho_cropped.tif'\n",
240+
"date_time_str = '2022-06-23 14:34:00'\n",
241+
"\n",
242+
"azimuth, zenith = get_sun_angles(file_path, date_time_str)\n",
243+
"\n",
244+
"# print(f'Sun azimuth: {azimuth}, Sun zenith: {zenith}')"
245+
]
246+
},
247+
{
248+
"cell_type": "code",
249+
"execution_count": null,
250+
"metadata": {},
251+
"outputs": [],
252+
"source": []
253+
}
254+
],
255+
"metadata": {
256+
"kernelspec": {
257+
"display_name": "cires-demo",
258+
"language": "python",
259+
"name": "python3"
260+
},
261+
"language_info": {
262+
"codemirror_mode": {
263+
"name": "ipython",
264+
"version": 3
265+
},
266+
"file_extension": ".py",
267+
"mimetype": "text/x-python",
268+
"name": "python",
269+
"nbconvert_exporter": "python",
270+
"pygments_lexer": "ipython3",
271+
"version": "3.10.8"
272+
},
273+
"orig_nbformat": 4
274+
},
275+
"nbformat": 4,
276+
"nbformat_minor": 2
277+
}

‎Slope_Aspect_Drone_Data.ipynb

+153
Large diffs are not rendered by default.

0 commit comments

Comments
 (0)