| blob | e1dd6cdceb157d4f8c0b22de609527c6bae64fdd |
1 {
2 "cells": [
3 {
4 "cell_type": "markdown",
5 "id": "f03213ea-7260-4f78-9b3b-99cc2ac18358",
6 "metadata": {},
7 "source": [
8 "# Geotiff Interpolation Tutorial\n",
9 "\n",
10 "The purpose of this notebook is to demonstrate interpolation from geotiff files to lon/lat pairs. In this case, the geotiff files are bands from HRRR grib2 files collected using `wrfxpy` methods on Alderaan. The lon/lat pairs are collected from RAWS stations with `Mesopy`, and assumed to be in the WGS84 standard coordinate system AKA EPSG 4326.\n",
11 "\n",
12 "The tiff files are saved with naming a convention that stores the UTC time info as well as the associated band. See:\n",
13 "\n",
14 "https://www.nco.ncep.noaa.gov/pmb/products/hrrr/hrrr.t00z.wrfprsf00.grib2.shtml"
15 ]
16 },
17 {
18 "cell_type": "markdown",
19 "id": "9f0e2346-5d88-4b9f-9b07-751238700af9",
20 "metadata": {},
21 "source": [
22 "# Setup "
23 ]
24 },
25 {
26 "cell_type": "code",
27 "execution_count": null,
28 "id": "110e36b1-61ea-4578-afc1-f774b7f61ee2",
29 "metadata": {},
30 "outputs": [],
31 "source": [
32 "import os\n",
33 "import os.path as osp\n",
34 "import numpy as np\n",
35 "from osgeo import gdal, osr\n",
36 "import matplotlib.pyplot as plt\n",
37 "import pandas as pd\n",
38 "import sys\n",
39 "sys.path.append('..')\n",
40 "from utils import retrieve_url"
41 ]
42 },
43 {
44 "cell_type": "code",
45 "execution_count": null,
46 "id": "3aa52e7f-fbb2-44ff-8308-a27ff93f8ef3",
47 "metadata": {},
48 "outputs": [],
49 "source": [
50 "import importlib\n",
51 "import utils\n",
52 "importlib.reload(utils)\n",
53 "from utils import retrieve_url"
54 ]
55 },
56 {
57 "cell_type": "code",
58 "execution_count": null,
59 "id": "c30fdba7-5a53-4bb9-8723-ae9253ae8036",
60 "metadata": {},
61 "outputs": [],
62 "source": [
63 "tfile = \"hrrr.t00z.wrfprsf00.585.tif\"\n",
64 "retrieve_url(\n",
65 " url = f\"https://demo.openwfm.org/web/data/fmda/tif/20240101/{tfile}\", \n",
66 " dest_path = f\"{tfile}\")"
67 ]
68 },
69 {
70 "cell_type": "code",
71 "execution_count": null,
72 "id": "5eaf3f83-a862-488d-a860-b2301a35ead3",
73 "metadata": {
74 "scrolled": true
75 },
76 "outputs": [],
77 "source": [
78 "# Get RAWS Station lat/lon\n",
79 "sts = pd.read_csv(\"C:/Users/jhirs/Documents/Projects/openwfm/notebooks/fmda/data/raws_stations_CO.csv\")\n",
80 "sts"
81 ]
82 },
83 {
84 "cell_type": "markdown",
85 "id": "236d8863-309c-4cbc-8b34-bd5cddda420b",
86 "metadata": {},
87 "source": [
88 "## Read Raster File"
89 ]
90 },
91 {
92 "cell_type": "code",
93 "execution_count": null,
94 "id": "8135266a-3638-4332-b79e-58fbca955f01",
95 "metadata": {},
96 "outputs": [],
97 "source": [
98 "# Extract data from tif file\n",
99 "ds = gdal.Open(tfile)\n",
100 "width = ds.RasterXSize\n",
101 "height = ds.RasterYSize\n",
102 "gt = ds.GetGeoTransform()\n",
103 "gp = ds.GetProjection()\n",
104 "# data = np.array(ds.ReadAsArray())"
105 ]
106 },
107 {
108 "cell_type": "code",
109 "execution_count": null,
110 "id": "b1195375-c6e8-4d57-b296-a5ed3f467a39",
111 "metadata": {},
112 "outputs": [],
113 "source": [
114 "print(width)\n",
115 "print(height)\n",
116 "print(gt)\n",
117 "print(gp)"
118 ]
119 },
120 {
121 "cell_type": "code",
122 "execution_count": null,
123 "id": "c10688f1-9649-44b6-a32f-9ae794548da0",
124 "metadata": {},
125 "outputs": [],
126 "source": [
127 "print(gt)"
128 ]
129 },
130 {
131 "cell_type": "code",
132 "execution_count": null,
133 "id": "1bdb1170-6d6d-4dba-84a0-2dbe8ecac11c",
134 "metadata": {},
135 "outputs": [],
136 "source": [
137 "print('Raster count: ' + str(ds.RasterCount))"
138 ]
139 },
140 {
141 "cell_type": "code",
142 "execution_count": null,
143 "id": "7c9ac9cd-a4ec-4f95-923f-ba6ec1f51d30",
144 "metadata": {},
145 "outputs": [],
146 "source": [
147 "band = ds.GetRasterBand(1)\n",
148 "data = band.ReadAsArray()"
149 ]
150 },
151 {
152 "cell_type": "markdown",
153 "id": "74b529ac-9ad4-407c-bb92-31828f4b43d9",
154 "metadata": {},
155 "source": [
156 "# Plot Raster File\n",
157 "\n",
158 "Using `imshow`, add a point at 100,100 just to demonstrate image indexing.\n",
159 "\n",
160 "source: https://www.geeksforgeeks.org/visualizing-tiff-file-using-matplotlib-and-gdal-using-python/"
161 ]
162 },
163 {
164 "cell_type": "code",
165 "execution_count": null,
166 "id": "0dcb3b37-d7d9-41ac-871b-bcc064cb2207",
167 "metadata": {
168 "scrolled": true
169 },
170 "outputs": [],
171 "source": [
172 "# Plot\n",
173 "plt.imshow(data)\n",
174 "plt.plot(100, 100, marker='o', color='blue', markersize=6)\n",
175 "plt.annotate(\"(100,100)\", (100,100), color='blue')"
176 ]
177 },
178 {
179 "cell_type": "markdown",
180 "id": "87ba2725-1512-4218-9bf7-e5e3082d4bb8",
181 "metadata": {},
182 "source": [
183 "# Nearest Neighbor Lat/Lon\n",
184 "\n",
185 "Source (nearest neighbor method): https://stackoverflow.com/questions/69034965/given-a-geotiff-file-how-does-one-find-the-single-pixel-closest-to-a-given-lati"
186 ]
187 },
188 {
189 "cell_type": "code",
190 "execution_count": null,
191 "id": "d332d1d5-4e2a-454c-bbf7-a8d5401296e7",
192 "metadata": {},
193 "outputs": [],
194 "source": [
195 "point_srs = osr.SpatialReference()\n",
196 "point_srs.ImportFromEPSG(4326) # hardcode for lon/lat\n",
197 "\n",
198 "# GDAL>=3: make sure it's x/y\n",
199 "# see https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn\n",
200 "point_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) \n",
201 "\n",
202 "file_srs = osr.SpatialReference()\n",
203 "file_srs.ImportFromWkt(gp)"
204 ]
205 },
206 {
207 "cell_type": "code",
208 "execution_count": null,
209 "id": "3399032e-2832-473d-a0f9-701a107303e1",
210 "metadata": {},
211 "outputs": [],
212 "source": [
213 "ct = osr.CoordinateTransformation(point_srs, file_srs)\n",
214 "\n",
215 "point_x = -105.002133 # lon\n",
216 "point_y = 39.746153 # lat\n",
217 "mapx, mapy, z = ct.TransformPoint(point_x, point_y) # output: coordinate pair (m)"
218 ]
219 },
220 {
221 "cell_type": "code",
222 "execution_count": null,
223 "id": "352bef20-7e53-4b80-8669-fea2b3c0767f",
224 "metadata": {},
225 "outputs": [],
226 "source": [
227 "mapx"
228 ]
229 },
230 {
231 "cell_type": "code",
232 "execution_count": null,
233 "id": "a4203c63-a0bb-421f-b476-b27ea9f6abbe",
234 "metadata": {},
235 "outputs": [],
236 "source": [
237 "mapy"
238 ]
239 },
240 {
241 "cell_type": "code",
242 "execution_count": null,
243 "id": "f93b0e4a-d632-4787-8edc-e736eb881e03",
244 "metadata": {},
245 "outputs": [],
246 "source": [
247 "z"
248 ]
249 },
250 {
251 "cell_type": "code",
252 "execution_count": null,
253 "id": "e7d90e9a-031f-49c9-9229-1da02125376f",
254 "metadata": {},
255 "outputs": [],
256 "source": [
257 "gt_inv = gdal.InvGeoTransform(gt)\n",
258 "pixel_x, pixel_y = gdal.ApplyGeoTransform(gt_inv, mapx, mapy)"
259 ]
260 },
261 {
262 "cell_type": "markdown",
263 "id": "b188fe33-3d41-437e-8eec-6b050a79bfb3",
264 "metadata": {},
265 "source": [
266 "We plot the image with the pixel annotated. The lon/lat pair is from a RAWS station near Colorado Springs, which matches the image below."
267 ]
268 },
269 {
270 "cell_type": "code",
271 "execution_count": null,
272 "id": "7bbfb4ef-11f9-4f01-a8cc-01ddd972a000",
273 "metadata": {},
274 "outputs": [],
275 "source": [
276 "plt.imshow(data)\n",
277 "\n",
278 "# Plot pixel translation with annotations\n",
279 "plt.plot(pixel_x,pixel_y, marker='o', color='blue', markersize=6)\n",
280 "plt.annotate(f\"Pixels: ({round(pixel_x, 2)}, {round(pixel_y, 2)})\", xy=(pixel_x, pixel_y),\n",
281 " xytext=(pixel_x-100, pixel_y-100), fontsize=8, color='lightblue')\n",
282 "plt.annotate(f\"Lon/Lat: ({round(point_x, 2)}, {round(point_y, 2)})\", xy=(pixel_x, pixel_y),\n",
283 " xytext=(pixel_x-100, pixel_y-50), fontsize=8, color='lightblue')"
284 ]
285 },
286 {
287 "cell_type": "markdown",
288 "id": "19bc96c4-5b62-4484-bcec-94ce9bec27a5",
289 "metadata": {},
290 "source": [
291 "After this point, the tutorial goes on to describe a method for nearest neighbor. This is just one form of interpolation, so various methods will be explored below.\n",
292 "\n",
293 "We will plot a zoomed in version of the pixels to demonstrate this."
294 ]
295 },
296 {
297 "cell_type": "code",
298 "execution_count": null,
299 "id": "242d8106-79c7-4948-a15d-6f4e4ddf0cec",
300 "metadata": {},
301 "outputs": [],
302 "source": [
303 "plt.imshow(data)\n",
304 "\n",
305 "# Plot pixel translation with annotations\n",
306 "plt.plot(pixel_x,pixel_y, marker='o', color='lightblue', markersize=6,\n",
307 " label=f\"({round(pixel_x, 2)}, {round(pixel_y, 2)})\")\n",
308 "\n",
309 "# Zoom in, set limits from plotted pixel\n",
310 "offset = 5\n",
311 "plt.xlim(pixel_x-offset, pixel_x+offset) # Note different +/- bc origin is topright\n",
312 "plt.ylim(pixel_y+offset,pixel_y-offset)\n",
313 "\n",
314 "# Plot 4 points bracketing target pixel\n",
315 "x1, y1=np.floor(pixel_x), np.floor(pixel_y)\n",
316 "x2, y2=np.floor(pixel_x), np.ceil(pixel_y)\n",
317 "x3, y3=np.ceil(pixel_x), np.floor(pixel_y)\n",
318 "x4, y4=np.ceil(pixel_x), np.ceil(pixel_y)\n",
319 "\n",
320 "plt.plot(x1,y1, marker='o', color='red', markersize=4, label='Bracketing Points')\n",
321 "plt.plot(x2,y2, marker='o', color='red', markersize=4)\n",
322 "plt.plot(x3,y3, marker='o', color='red', markersize=4)\n",
323 "plt.plot(x4,y4, marker='o', color='red', markersize=4)\n",
324 "\n",
325 "plt.legend()"
326 ]
327 },
328 {
329 "cell_type": "markdown",
330 "id": "a3362432-1469-47e2-aa77-d5eb4aede959",
331 "metadata": {},
332 "source": [
333 "# Interpolation Methods"
334 ]
335 },
336 {
337 "cell_type": "markdown",
338 "id": "c0b16dbc-d206-414f-aa55-9c91482f3489",
339 "metadata": {},
340 "source": [
341 "## Nearest Neighbor\n",
342 "\n",
343 "The tutorial linked above simply rounds the pixel x and y coordinates to the nearest pixels and takes the value from that grid location. This is mathematically equivalent to an L1 nearest neighbor, or manhattan distance minimization."
344 ]
345 },
346 {
347 "cell_type": "code",
348 "execution_count": null,
349 "id": "7dc6b796-33fc-461b-be5c-5b175d8a1d1f",
350 "metadata": {},
351 "outputs": [],
352 "source": [
353 "# round to pixel\n",
354 "x_l1 = round(pixel_x)\n",
355 "y_l1 = round(pixel_y)"
356 ]
357 },
358 {
359 "cell_type": "code",
360 "execution_count": null,
361 "id": "e4bf3262-752d-46f5-a4f8-9f90a68a6b4c",
362 "metadata": {},
363 "outputs": [],
364 "source": [
365 "plt.imshow(data)\n",
366 "\n",
367 "# Plot pixel translation with annotations\n",
368 "plt.plot(pixel_x,pixel_y, marker='o', color='blue', markersize=6,\n",
369 " label=f\"({round(pixel_x, 2)}, {round(pixel_y, 2)})\")\n",
370 "\n",
371 "# Zoom in, set limits from plotted pixel\n",
372 "offset = 5\n",
373 "plt.xlim(pixel_x-offset, pixel_x+offset) # Note different +/- bc origin is topright\n",
374 "plt.ylim(pixel_y+offset,pixel_y-offset)\n",
375 "\n",
376 "# Plot 4 points bracketing target pixel\n",
377 "x1, y1=np.floor(pixel_x), np.floor(pixel_y)\n",
378 "x2, y2=np.floor(pixel_x), np.ceil(pixel_y)\n",
379 "x3, y3=np.ceil(pixel_x), np.floor(pixel_y)\n",
380 "x4, y4=np.ceil(pixel_x), np.ceil(pixel_y)\n",
381 "\n",
382 "plt.plot(x1,y1, marker='o', color='red', markersize=4, label='Bracketing Points')\n",
383 "plt.plot(x2,y2, marker='o', color='red', markersize=4)\n",
384 "plt.plot(x3,y3, marker='o', color='red', markersize=4)\n",
385 "plt.plot(x4,y4, marker='o', color='red', markersize=4)\n",
386 "\n",
387 "# Plot interpolated pixel\n",
388 "plt.plot(x_l1,y_l1, marker='o', color='purple', markersize=5,\n",
389 " label='Interpolated')\n",
390 "interp_val = data[y_l1, x_l1]\n",
391 "plt.annotate(\"Interp. Value=\"+str(round(interp_val, 5)), xy=(x_l1, y_l1), xytext=(x_l1-1, y_l1-2), color='purple')\n",
392 "\n",
393 "plt.title(\"L1 Nearest Neighbor\")\n",
394 "plt.legend()"
395 ]
396 },
397 {
398 "cell_type": "code",
399 "execution_count": null,
400 "id": "f75be7c9-1957-4daf-8ec8-368aab6beb90",
401 "metadata": {},
402 "outputs": [],
403 "source": []
404 },
405 {
406 "cell_type": "markdown",
407 "id": "ec7f3ce4-96b1-49b8-b1c5-46f31bad88c6",
408 "metadata": {},
409 "source": [
410 "## Nearest Neighbor (Euclidean)\n",
411 "\n",
412 "In `wrfxpy`, the function `find_closest_grid_point` is defined to find the L2 nearest neighbor, which finds the minimum sum of squared distance (Euclidean norm).\n",
413 "\n",
414 "https://github.com/openwfm/wrfxpy/blob/master/src/utils.py#L529\n",
415 "\n",
416 "NOTE: very slow implementation, but I wanted to reproduce the method clearly. In this case, the interpolated value is the same as L1."
417 ]
418 },
419 {
420 "cell_type": "code",
421 "execution_count": null,
422 "id": "09d65190-a4d0-4e8a-b6cf-494b636ce37a",
423 "metadata": {},
424 "outputs": [],
425 "source": [
426 "x = np.arange(0, band.XSize)\n",
427 "y = np.arange(0, band.YSize)\n",
428 "pixels = [(xx, yy) for xx in x for yy in y]\n",
429 "d = np.zeros(len(pixels))\n",
430 "for i in range(0, len(pixels)):\n",
431 " p = pixels[i]\n",
432 " d[i] = (pixel_x - p[0])**2 + (pixel_y - p[1])**2\n",
433 "\n",
434 "nearest = pixels[np.argmin(d)]"
435 ]
436 },
437 {
438 "cell_type": "code",
439 "execution_count": null,
440 "id": "c8165a6e-f63b-4b9e-841d-27570829b3ee",
441 "metadata": {},
442 "outputs": [],
443 "source": [
444 "plt.imshow(data)\n",
445 "\n",
446 "# Plot pixel translation with annotations\n",
447 "plt.plot(pixel_x,pixel_y, marker='o', color='blue', markersize=6,\n",
448 " label=f\"({round(pixel_x, 2)}, {round(pixel_y, 2)})\")\n",
449 "\n",
450 "# Zoom in, set limits from plotted pixel\n",
451 "offset = 5\n",
452 "plt.xlim(pixel_x-offset, pixel_x+offset) # Note different +/- bc origin is topright\n",
453 "plt.ylim(pixel_y+offset,pixel_y-offset)\n",
454 "\n",
455 "# Plot 4 points bracketing target pixel\n",
456 "x1, y1=np.floor(pixel_x), np.floor(pixel_y)\n",
457 "x2, y2=np.floor(pixel_x), np.ceil(pixel_y)\n",
458 "x3, y3=np.ceil(pixel_x), np.floor(pixel_y)\n",
459 "x4, y4=np.ceil(pixel_x), np.ceil(pixel_y)\n",
460 "\n",
461 "plt.plot(x1,y1, marker='o', color='red', markersize=4, label='Bracketing Pixels')\n",
462 "plt.plot(x2,y2, marker='o', color='red', markersize=4)\n",
463 "plt.plot(x3,y3, marker='o', color='red', markersize=4)\n",
464 "plt.plot(x4,y4, marker='o', color='red', markersize=4)\n",
465 "\n",
466 "# find nearest L2 pixel\n",
467 "plt.plot(nearest[0],nearest[1], marker='o', color='purple', markersize=5,\n",
468 " label='Interpolated')\n",
469 "interp_val = data[nearest[1], nearest[0]]\n",
470 "plt.annotate(\"Interp. Value=\"+str(round(interp_val, 5)), xy=(x_l1, y_l1), xytext=(x_l1-1, y_l1-2), color='purple')\n",
471 "\n",
472 "\n",
473 "plt.title(\"L2 Nearest Neighbor\")\n",
474 "plt.legend()"
475 ]
476 },
477 {
478 "cell_type": "markdown",
479 "id": "a6c2987c-99ed-4313-85bc-47630f111775",
480 "metadata": {},
481 "source": [
482 "## Average Bracketing Points\n",
483 "\n",
484 "Given the 4 bracketing points, interpolate the mean value."
485 ]
486 },
487 {
488 "cell_type": "code",
489 "execution_count": null,
490 "id": "8d0124de-dd5a-4913-acb9-b2b81ba2bfcc",
491 "metadata": {},
492 "outputs": [],
493 "source": [
494 "plt.imshow(data)\n",
495 "\n",
496 "# Plot pixel translation with annotations\n",
497 "plt.plot(pixel_x,pixel_y, marker='o', color='blue', markersize=6,\n",
498 " label=f\"({round(pixel_x, 2)}, {round(pixel_y, 2)})\")\n",
499 "\n",
500 "# Zoom in, set limits from plotted pixel\n",
501 "offset = 5\n",
502 "plt.xlim(pixel_x-offset, pixel_x+offset) # Note different +/- bc origin is topright\n",
503 "plt.ylim(pixel_y+offset,pixel_y-offset)\n",
504 "\n",
505 "# Plot 4 points bracketing target pixel\n",
506 "x1, y1=int(np.floor(pixel_x)), int(np.floor(pixel_y))\n",
507 "x2, y2=int(np.floor(pixel_x)), int(np.ceil(pixel_y))\n",
508 "x3, y3=int(np.ceil(pixel_x)), int(np.floor(pixel_y))\n",
509 "x4, y4=int(np.ceil(pixel_x)), int(np.ceil(pixel_y))\n",
510 "\n",
511 "plt.plot(x1,y1, marker='o', color='red', markersize=4, label='Bracketing Points')\n",
512 "plt.plot(x2,y2, marker='o', color='red', markersize=4)\n",
513 "plt.plot(x3,y3, marker='o', color='red', markersize=4)\n",
514 "plt.plot(x4,y4, marker='o', color='red', markersize=4)\n",
515 "\n",
516 "interp_val = np.mean([data[y1, x1], data[y2, x2], data[y3, x3], data[y4, 4]])\n",
517 "plt.annotate(\"Interp. Value=\"+str(round(interp_val, 5)), xy=(pixel_x, pixel_y), xytext=(pixel_x-1, pixel_y-2), color='purple')\n",
518 "\n",
519 "\n",
520 "plt.title(\"Average 4 Nearest Neighbors\")\n",
521 "plt.legend()"
522 ]
523 },
524 {
525 "cell_type": "markdown",
526 "id": "80322dc5-efa2-4241-8f08-ffc58828f6f2",
527 "metadata": {},
528 "source": [
529 "## Scipy griddata interpolation package"
530 ]
531 },
532 {
533 "cell_type": "markdown",
534 "id": "76f5146e-b6ed-4b2b-8da6-927a747c0ed0",
535 "metadata": {},
536 "source": [
537 "https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html"
538 ]
539 },
540 {
541 "cell_type": "code",
542 "execution_count": null,
543 "id": "361bdffc-13e8-4366-8b4e-4d9292f135ae",
544 "metadata": {},
545 "outputs": [],
546 "source": [
547 "from scipy.interpolate import griddata"
548 ]
549 },
550 {
551 "cell_type": "markdown",
552 "id": "5e969ec7-dd62-4a70-bfaf-76c70e70a02d",
553 "metadata": {},
554 "source": [
555 "## WRF-Fire Method\n",
556 "\n",
557 "The WRF-Fire interpolation method combines a nearest neighbors approach and a linear interpolation. The steps are:\n",
558 "* Find nearest neighbor\n",
559 "* Construct 3x3 array centered at nearest neighbor\n",
560 "* Run linear interpolation (matlab's `scatteredInterpolant` method)\n",
561 "\n",
562 "The python equivalent of `scatteredInterpolant` is `scipy.interpolate.griddata` (according to StackOverflow and ChatGPT).\n",
563 "\n",
564 "WRF-Fire Sources:\n",
565 "\n",
566 "https://github.com/openwfm/wrf-fire-matlab/blob/master/vis/ts_at.m\n",
567 "\r\n",
568 "https://github.com/openwfm/wrf-fire-matlab/blob/master/vis/ts_at_test.m\r\n"
569 ]
570 },
571 {
572 "cell_type": "markdown",
573 "id": "6da0597f-b393-4d27-9e2a-e522bf8506f3",
574 "metadata": {},
575 "source": [
576 "## Validate Geolocation\n",
577 "\n",
578 "Goals:\n",
579 "\n",
580 "* create two fake geotiff\n",
581 "* create data array with known lon/lat values\n",
582 "* evaluate and you should get coord back exactly\n",
583 "\n",
584 "Steps here:\n",
585 "\n",
586 "* Convert pixel indices (e.g. (0,0), (0,1)...) to lat lon\n",
587 "* Save lats and lons in array of same dimension as raster band\n",
588 "* Save as geotiff\n",
589 "* Read that file back in to do steps above\n",
590 "\n",
591 "\n",
592 "Modified from Source: https://stackoverflow.com/questions/59052516/find-lat-long-coordinates-from-pixel-point-in-geotiff-using-python-and-gdaltly"
593 ]
594 },
595 {
596 "cell_type": "code",
597 "execution_count": null,
598 "id": "7d3560ad-6e14-4686-8ca3-cac6816b77e4",
599 "metadata": {},
600 "outputs": [],
601 "source": [
602 "# Generate Arrays of pixel indices\n",
603 "lons = np.zeros((1059, 1799))\n",
604 "lats = np.zeros((1059, 1799))"
605 ]
606 },
607 {
608 "cell_type": "code",
609 "execution_count": null,
610 "id": "b356ed24-b99c-4a3c-b783-9b708a34555c",
611 "metadata": {},
612 "outputs": [],
613 "source": [
614 "from osgeo import osr, ogr, gdal\n",
615 "\n",
616 "def pixel_to_world(geo_matrix, x, y):\n",
617 " # Given geotransform info of a geotiff file and an (x,y) pixel coord pair, return the coord pair that matches the geotiff in meters\n",
618 " # Inputs: \n",
619 " # geomatrix: output of ds.GetGeoTransform() for given geotiff file\n",
620 " # tuple of length 6 contains: \n",
621 " # A geotransform consists in a set of 6 coefficients\n",
622 " # GT(0) x-coordinate of the upper-left corner of the upper-left pixel.\n",
623 " # GT(1) w-e pixel resolution / pixel width.\n",
624 " # GT(2) row rotation (typically zero).\n",
625 " # GT(3) y-coordinate of the upper-left corner of the upper-left pixel.\n",
626 " # GT(4) column rotation (typically zero).\n",
627 " # GT(5) n-s pixel resolution / pixel height (negative value for a north-up image).\n",
628 " # x: pixel index x coord (1)\n",
629 " # y: pixel index y coord (1)\n",
630 " # Return: coordinates of same point as given x,y as offset from UL (m)\n",
631 " # Example: pixel_to_world(mat, 0, 0) returns UL x,y from geotiff\n",
632 " \n",
633 " ul_x = geo_matrix[0]\n",
634 " ul_y = geo_matrix[3]\n",
635 " x_dist = geo_matrix[1]\n",
636 " y_dist = geo_matrix[5]\n",
637 " _x = x * x_dist + ul_x\n",
638 " _y = y * y_dist + ul_y\n",
639 " return _x, _y\n",
640 "\n",
641 "\n",
642 "def build_transform_inverse(dataset, EPSG):\n",
643 " # Given gdal dataset and target EPSG, return transformation function that transforms meter coord pairs to pixel coord pairs \n",
644 " # Inputs:\n",
645 " # dataset: geotiff file\n",
646 " # EPSG: integer\n",
647 " source = osr.SpatialReference(wkt=dataset.GetProjection())\n",
648 " target = osr.SpatialReference()\n",
649 " target.ImportFromEPSG(EPSG)\n",
650 " return osr.CoordinateTransformation(source, target)\n",
651 "\n",
652 "def world_to_epsg(wx, wy, trans):\n",
653 " # Inputs:\n",
654 " # wx, wy: output of build_transform_inverse\n",
655 " # wx: x coordinate (m) related to geotiff reference point\n",
656 " # wy: y coordinate (m) related to geotiff reference point\n",
657 " # transform: function to transform to given epsg, function type is osgeo.osr.CoordinateTransformation\n",
658 " # Return: \n",
659 " # point from osgeo Geometry object\n",
660 " point = ogr.Geometry(ogr.wkbPoint)\n",
661 " point.AddPoint(wx, wy)\n",
662 " point.Transform(trans)\n",
663 " return point\n",
664 "\n",
665 "def find_spatial_coordinate_from_pixel(dataset, x, y, transform=None, epsg=4326):\n",
666 " # Given gdal dataset, target x y pixel pair, and EPSG, return the EPSG defined coordinate pair \n",
667 " # dataset: gdal dataset, from geotiff file\n",
668 " # x (int): pixel x index \n",
669 " # y (int): pixel y index \n",
670 " ## Upper left corner is often (0,0)\n",
671 " # transform: transform inverse. output of build_transform_inverse, default none and it calculates from epsg\n",
672 " # supply transform to save computational time\n",
673 " # epsg: default 4326 (WGS84)\n",
674 " # Return: coord pair in given epsg, eg lat/lon (floats)\n",
675 " if transform is None:\n",
676 " transform = build_transform_inverse(ds, epsg)\n",
677 " world_x, world_y = pixel_to_world(dataset.GetGeoTransform(), x, y)\n",
678 " point = world_to_epsg(world_x, world_y, transform)\n",
679 " return point.GetX(), point.GetY()\n",
680 "\n",
681 "ds = gdal.Open(tfile)"
682 ]
683 },
684 {
685 "cell_type": "markdown",
686 "id": "2f86b706-7eea-49d5-8e32-c5f332b7f439",
687 "metadata": {},
688 "source": [
689 "We want this process to match information from the gdalinfo in the geotiff file. The command line `gdalinfo hrrr.t00z.wrfprsf00.616.tif` returns info on the corner coordinates and center (not (0,0) for some reason). The following output should thus match the return of that command.\n",
690 "\n",
691 "`gdalinfo` returns coordinates in the \"degrees minute second\" format, so we need to convert decimal degrees to this format to make sure it matches. \n"
692 ]
693 },
694 {
695 "cell_type": "code",
696 "execution_count": null,
697 "id": "9a278659-2196-4055-8f0f-9f79196d8c26",
698 "metadata": {},
699 "outputs": [],
700 "source": [
701 "## Test world_to_epsg\n",
702 "## For 'center' in file (-520.143, -306.153), not zero for some reason\n",
703 "## Should return ( 97d30'21.52\"W, 38d29'50.09\"N)\n",
704 "trans = build_transform_inverse(ds, 4326)\n",
705 "pt = world_to_epsg(-520.143, -306.153, trans)\n",
706 "print(f\"({pt.GetX()},{pt.GetY()})\")"
707 ]
708 },
709 {
710 "cell_type": "code",
711 "execution_count": null,
712 "id": "58d1ff37-dbc2-45dd-9b09-c13b6e39a6c6",
713 "metadata": {},
714 "outputs": [],
715 "source": [
716 "## Test on a corner coords\n",
717 "print(find_spatial_coordinate_from_pixel(ds, 0, 0)) # upper left\n",
718 "print(find_spatial_coordinate_from_pixel(ds, 0, 1058)) # upper right\n",
719 "print(find_spatial_coordinate_from_pixel(ds, 1798, 0)) # lower left\n",
720 "print(find_spatial_coordinate_from_pixel(ds, 1798, 1058)) # lower right"
721 ]
722 },
723 {
724 "cell_type": "markdown",
725 "id": "80950022-acbb-4c79-a544-a5c805e28da4",
726 "metadata": {},
727 "source": [
728 "`gdalinfo` returns coordinates with the format \"degree minute second(direction)\", so we convert the values from `gdalinfo` to check...\n",
729 "\n",
730 "Source: stackexchange ___"
731 ]
732 },
733 {
734 "cell_type": "code",
735 "execution_count": null,
736 "id": "d0fc3849-1dac-4fc0-b2a6-9cf270932fff",
737 "metadata": {},
738 "outputs": [],
739 "source": [
740 "def decdeg2dms(dd, nsec=4):\n",
741 " # Given decimal degree, turn Degree Minute Second format\n",
742 " # dd: decimal lat/lon coord\n",
743 " # nsec: number of digits to round second\n",
744 " mult = -1 if dd < 0 else 1\n",
745 " mnt,sec = divmod(abs(dd)*3600, 60)\n",
746 " deg,mnt = divmod(mnt, 60)\n",
747 " deg=mult*deg\n",
748 " mnt=mult*mnt\n",
749 " sec=np.round(mult*sec, nsec)\n",
750 " \n",
751 " date_str = f\"{deg}d {mnt}'{sec}\"\n",
752 " return date_str"
753 ]
754 },
755 {
756 "cell_type": "markdown",
757 "id": "9050dd35-202d-44f8-ab6f-744910839515",
758 "metadata": {},
759 "source": [
760 "Check on coordinates listed in `gdalinfo`\n",
761 "\n",
762 "Using coords in meters, convert first with `world_to_epsg` then evaluate with `decdeg2dms`. Note: negative coordinates correspond to degrees West or degrees South from point of origin (in HRRR files the origin is in continental US)."
763 ]
764 },
765 {
766 "cell_type": "code",
767 "execution_count": null,
768 "id": "6ab633c2-bfa4-41b3-b7b0-1f8646261414",
769 "metadata": {},
770 "outputs": [],
771 "source": [
772 "# Print lon/lat coords\n",
773 "ul = world_to_epsg(-2699020.143, 1588193.847, trans) # Upper left\n",
774 "print(f\"({ul.GetX()},{ul.GetY()})\")\n",
775 "ll = world_to_epsg(-2699020.143,-1588806.153, trans) # Lower left\n",
776 "print(f\"({ll.GetX()},{ll.GetY()})\")\n",
777 "ur = world_to_epsg(2697979.857, 1588193.847, trans) # Upper right\n",
778 "print(f\"({ur.GetX()},{ur.GetY()})\")\n",
779 "lr = world_to_epsg(2697979.857,-1588806.153, trans) # Lower right\n",
780 "print(f\"({lr.GetX()},{lr.GetY()})\")\n",
781 "center = world_to_epsg(-520.143, -306.153, trans) # center\n",
782 "print(f\"({center.GetX()},{center.GetY()})\")"
783 ]
784 },
785 {
786 "cell_type": "code",
787 "execution_count": null,
788 "id": "ab76e16c-c5ad-420a-b43a-5c63d1b6a20e",
789 "metadata": {},
790 "outputs": [],
791 "source": [
792 "# Convert to dms\n",
793 "print(f\"({decdeg2dms(ul.GetY())}, {decdeg2dms(ul.GetX())})\")\n",
794 "print(f\"({decdeg2dms(ll.GetY())}, {decdeg2dms(ll.GetX())})\")\n",
795 "print(f\"({decdeg2dms(ur.GetY())}, {decdeg2dms(ur.GetX())})\")\n",
796 "print(f\"({decdeg2dms(lr.GetY())}, {decdeg2dms(lr.GetX())})\")\n",
797 "print(f\"({decdeg2dms(center.GetY())}, {decdeg2dms(center.GetX())})\")"
798 ]
799 },
800 {
801 "cell_type": "markdown",
802 "id": "8c769665-6014-4f7c-b83b-31fc67e9a658",
803 "metadata": {},
804 "source": [
805 "### Build dummy geotiff files\n",
806 "\n",
807 "Create 2 geotiff files with latitude/longitude values corresponding to grid pixel locations."
808 ]
809 },
810 {
811 "cell_type": "code",
812 "execution_count": null,
813 "id": "a9944588-50df-4f47-8af3-cdf4727cfd1a",
814 "metadata": {},
815 "outputs": [],
816 "source": [
817 "# dimensions of numpy ndarray stored as values of geotiff file\n",
818 "print(np.shape(data)) "
819 ]
820 },
821 {
822 "cell_type": "code",
823 "execution_count": null,
824 "id": "d084fb21-0de9-4b27-b7ce-aa9f70afba25",
825 "metadata": {},
826 "outputs": [],
827 "source": [
828 "# Initialize empty arrays\n",
829 "lons=np.zeros(np.shape(data))\n",
830 "lats=np.zeros(np.shape(data))"
831 ]
832 },
833 {
834 "cell_type": "code",
835 "execution_count": null,
836 "id": "d762a7d7-d322-4063-8704-2a5cd7ebb684",
837 "metadata": {},
838 "outputs": [],
839 "source": [
840 "# get transformation once and reuse\n",
841 "transform = build_transform_inverse(ds, EPSG=4326)\n",
842 "# Loop over indices and fill\n",
843 "for i in range(0, np.shape(lons)[0]): # iterate i over x coord (longitude)\n",
844 " for j in range(0, np.shape(lons)[1]): # iterate j over y coord (latitude)\n",
845 " coord = find_spatial_coordinate_from_pixel(ds, j, i, transform=transform) # note order flip is intentional\n",
846 " lats[i,j]=coord[0]\n",
847 " lons[i,j]=coord[1]"
848 ]
849 },
850 {
851 "cell_type": "markdown",
852 "id": "b4efef93-bee9-4c87-a353-f86eab767954",
853 "metadata": {},
854 "source": [
855 "Function to write geotiff file.\n",
856 "\n",
857 "source: https://here.isnew.info/how-to-save-a-numpy-array-as-a-geotiff-file-using-gdal.html"
858 ]
859 },
860 {
861 "cell_type": "code",
862 "execution_count": null,
863 "id": "927270f5-af27-4bbd-a9a9-8f4f73ddfc16",
864 "metadata": {},
865 "outputs": [],
866 "source": [
867 "def write_geotiff(filename, arr, in_ds): \n",
868 " # Given file name, data array, and reference gdal Dataset, write data array as geotiff file with geolocation from reference\n",
869 " # Inputs:\n",
870 " # filename: output file name, expecting .tif extension\n",
871 " # arr: numpy ndarray\n",
872 " # in_ds: gdal dataset with geolocation info\n",
873 "\n",
874 " print(\"Writing \"+filename)\n",
875 " \n",
876 " if arr.dtype == np.float32:\n",
877 " arr_type = gdal.GDT_Float32\n",
878 " else:\n",
879 " arr_type = gdal.GDT_Int32\n",
880 "\n",
881 " driver = gdal.GetDriverByName(\"GTiff\")\n",
882 " out_ds = driver.Create(filename, arr.shape[1], arr.shape[0], 1, arr_type)\n",
883 " out_ds.SetProjection(in_ds.GetProjection())\n",
884 " out_ds.SetGeoTransform(in_ds.GetGeoTransform())\n",
885 " band = out_ds.GetRasterBand(1)\n",
886 " band.WriteArray(arr)\n",
887 " band.FlushCache()\n",
888 " band.ComputeStatistics(False)"
889 ]
890 },
891 {
892 "cell_type": "code",
893 "execution_count": null,
894 "id": "c22e54fb-e65c-4cf7-87c8-0af326300943",
895 "metadata": {},
896 "outputs": [],
897 "source": [
898 "# Write output as geotiff file\n",
899 "# osp.join(os.getcwd(), \"lons.tif\")\n",
900 "# osp.join(os.getcwd(), \"lats.tif\")\n",
901 "write_geotiff(osp.join(os.getcwd(), \"lons.tif\"), lons, ds)\n",
902 "write_geotiff(osp.join(os.getcwd(), \"lats.tif\"), lats, ds)"
903 ]
904 },
905 {
906 "cell_type": "markdown",
907 "id": "6b2230db-0699-4152-aa53-fbc87a928139",
908 "metadata": {},
909 "source": [
910 "### Validate geolocation procedure\n",
911 "\n",
912 "Using geotiff files written above, use the procedure for plucking nearest points on the lat/lon pairs from RAWS stations. Given input lat/lon, expect the same values back."
913 ]
914 },
915 {
916 "cell_type": "code",
917 "execution_count": null,
918 "id": "a4a4e32a-9131-412d-8fde-3a251ba465dd",
919 "metadata": {},
920 "outputs": [],
921 "source": [
922 "ds_lon = gdal.Open(osp.join(os.getcwd(), \"lons.tif\"))\n",
923 "ds_lat = gdal.Open(osp.join(os.getcwd(), \"lats.tif\"))"
924 ]
925 },
926 {
927 "cell_type": "code",
928 "execution_count": null,
929 "id": "d2664341-88d4-4a26-a417-eeb3f9322306",
930 "metadata": {},
931 "outputs": [],
932 "source": [
933 "ds = gdal.Open(tfile)\n",
934 "gt = ds_lon.GetGeoTransform()\n",
935 "gp = ds_lon.GetProjection()\n",
936 "\n",
937 "# Check values the same for lat file\n",
938 "print('Geotrans matches: '+str(gt == ds_lon.GetGeoTransform()))\n",
939 "print('Proj matches: '+str(gp == ds_lat.GetProjection()))"
940 ]
941 },
942 {
943 "cell_type": "code",
944 "execution_count": null,
945 "id": "cb48ccd8-bb12-4e0d-81ea-48eeb29d7504",
946 "metadata": {},
947 "outputs": [],
948 "source": [
949 "# Get Projection info\n",
950 "point_srs = osr.SpatialReference()\n",
951 "point_srs.ImportFromEPSG(4326) # hardcode for lon/lat\n",
952 "# GDAL>=3: make sure it's x/y\n",
953 "# see https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn\n",
954 "point_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) \n",
955 "file_srs = osr.SpatialReference()\n",
956 "file_srs.ImportFromWkt(gp)\n",
957 "ct = osr.CoordinateTransformation(point_srs, file_srs)"
958 ]
959 },
960 {
961 "cell_type": "code",
962 "execution_count": null,
963 "id": "e65462e0-7a33-40c8-9eb7-4a3c793652aa",
964 "metadata": {},
965 "outputs": [],
966 "source": [
967 "lons = ds_lon.GetRasterBand(1).ReadAsArray()\n",
968 "lats = ds_lat.GetRasterBand(1).ReadAsArray()"
969 ]
970 },
971 {
972 "cell_type": "markdown",
973 "id": "8140fcb7-fc99-425d-84bf-74e0014731a6",
974 "metadata": {},
975 "source": [
976 "Loop through some RAWS sites. Print out known lat/lon coord and then the value from the geotiff files, expect the output to be the same."
977 ]
978 },
979 {
980 "cell_type": "code",
981 "execution_count": null,
982 "id": "8e093ec2-2db2-4d6c-bd1c-db906312c99a",
983 "metadata": {},
984 "outputs": [],
985 "source": [
986 "# Helper function to define nearest neighbor approach\n",
987 "def interp_l1(x, y):\n",
988 " x = round(x)\n",
989 " y = round(y)\n",
990 " return x, y\n",
991 "\n",
992 "# loop through first N station coords\n",
993 "for i in range(0, 10):\n",
994 " print(\"~\"*35)\n",
995 " point_x = sts['lon'][i] # lon\n",
996 " point_y = sts['lat'][i] # lat\n",
997 " print(f\"RAWS Station {sts['STID'][i]} lon/lat: ({point_x}, {point_y})\")\n",
998 "\n",
999 " mapx, mapy, z = ct.TransformPoint(np.float64(point_x), np.float64(point_y))\n",
1000 " gt_inv = gdal.InvGeoTransform(gt)\n",
1001 " pixel_x, pixel_y = gdal.ApplyGeoTransform(gt_inv, mapx, mapy)\n",
1002 " x, y = interp_l1(pixel_x, pixel_y)\n",
1003 " print(f\"Fitted: ({lons[y, x]}, {lats[y,x]})\")\n"
1004 ]
1005 },
1006 {
1007 "cell_type": "code",
1008 "execution_count": null,
1009 "id": "119b2302-1aef-4d6c-bb82-fb512a58b9a4",
1010 "metadata": {},
1011 "outputs": [],
1012 "source": []
1013 },
1014 {
1015 "cell_type": "code",
1016 "execution_count": null,
1017 "id": "6e40fd62-31e2-461f-b32a-bdc5a26a9e2b",
1018 "metadata": {},
1019 "outputs": [],
1020 "source": []
1021 },
1022 {
1023 "cell_type": "code",
1024 "execution_count": null,
1025 "id": "aa402e2c-318a-4041-a20b-61f3d941aa54",
1026 "metadata": {},
1027 "outputs": [],
1028 "source": []
1029 }
1030 ],
1031 "metadata": {
1032 "kernelspec": {
1033 "display_name": "Python 3 (ipykernel)",
1034 "language": "python",
1035 "name": "python3"
1036 },
1037 "language_info": {
1038 "codemirror_mode": {
1039 "name": "ipython",
1040 "version": 3
1041 },
1042 "file_extension": ".py",
1043 "mimetype": "text/x-python",
1044 "name": "python",
1045 "nbconvert_exporter": "python",
1046 "pygments_lexer": "ipython3",
1047 "version": "3.12.5"
1048 }
1049 },
1050 "nbformat": 4,
1051 "nbformat_minor": 5
1052 }