Skip to content

Commit 4bac637

Browse files
weiji14weiji14
authored
Finishing touches (#10)
Add boto3 and scikit-image libraries for the extra Landsat 8 lessons and exploration time in Tutorial 2. Also added exploration time section to Tutorial 3 and a schedule timeline on what we'll be covering in the workshop!
1 parent d6f5b86 commit 4bac637

File tree

4 files changed

+375
-6
lines changed

4 files changed

+375
-6
lines changed

‎2_pydata_and_pygmt.ipynb‎

Lines changed: 247 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,11 @@
2424
"We'll start by importing the Python libraries we'll be using.\n",
2525
"Following common conventions, we'll abbreviate\n",
2626
"[pandas](https://pandas.pydata.org) as `pd` and\n",
27-
"[xarray](http://xarray.pydata.org) as `xr`."
27+
"[xarray](http://xarray.pydata.org) as `xr`.\n",
28+
"A little later into this tutorial,\n",
29+
"we'll also make use of\n",
30+
"[rioxarray](https://corteva.github.io/rioxarray) and\n",
31+
"[scikit-image](https://scikit-image.org)."
2832
]
2933
},
3034
{
@@ -35,6 +39,8 @@
3539
"source": [
3640
"import pandas as pd\n",
3741
"import pygmt\n",
42+
"import rioxarray\n",
43+
"import skimage.exposure\n",
3844
"import xarray as xr"
3945
]
4046
},
@@ -301,10 +307,8 @@
301307
"cell_type": "markdown",
302308
"metadata": {},
303309
"source": [
304-
"The next steps are just some additional preprocessing.\n",
305-
"We'll remove the NaN values so that they don't get plotted,\n",
306-
"and also sort our grid so that the x and y dimensions\n",
307-
"go in ascending order (i.e. West to East, South to North)."
310+
"The next step is just some additional preprocessing.\n",
311+
"We'll remove the NaN values so that they don't get plotted."
308312
]
309313
},
310314
{
@@ -314,7 +318,6 @@
314318
"outputs": [],
315319
"source": [
316320
"grid = grid.where(grid != grid.nodatavals, drop=True)\n",
317-
"grid = grid.sortby(variables=list(grid.dims))\n",
318321
"grid"
319322
]
320323
},
@@ -359,6 +362,244 @@
359362
" "
360363
]
361364
},
365+
{
366+
"cell_type": "markdown",
367+
"metadata": {},
368+
"source": [
369+
"## **2.3 Looking at true colour bands**\n",
370+
"\n",
371+
"What if your raster grid is not a DEM,\n",
372+
"but a multi-band satellite image instead?\n",
373+
"Well, grdimage also accepts RGB bands,\n",
374+
"but there are a couple of tricks we'll need to learn for things to work properly.\n",
375+
"\n",
376+
"Let's have a go at plotting a Landsat 8 Image of good ol' Wellington (and Cook Strait) this time.\n",
377+
"To give you a quick sense of the place from space, we'll do a quick preview of the JPG image."
378+
]
379+
},
380+
{
381+
"cell_type": "code",
382+
"execution_count": null,
383+
"metadata": {},
384+
"outputs": [],
385+
"source": [
386+
"fig = pygmt.Figure()\n",
387+
"fig.grdimage(\n",
388+
" grid=\"https://landsat-pds.s3.amazonaws.com/c1/L8/073/089/LC08_L1TP_073089_20191021_20191030_01_T1/LC08_L1TP_073089_20191021_20191030_01_T1_thumb_large.jpg\",\n",
389+
" frame=True,\n",
390+
")\n",
391+
"fig.show()"
392+
]
393+
},
394+
{
395+
"cell_type": "markdown",
396+
"metadata": {},
397+
"source": [
398+
"### **Getting some cloud optimized GeoTIFFs**\n",
399+
"\n",
400+
"Landsat 8 data is provided as [cloud optimized GeoTIFFs (COGs)](https://www.cogeo.org/)\n",
401+
"which makes it easy to pull down just what we need from the web.\n",
402+
"I recommend using [Remote Pixel's Satellite Search API](https://search.remotepixel.ca) to search for images.\n",
403+
"\n",
404+
"We'll switch from xarray to [rioxarray](https://corteva.github.io/rioxarray) this time,\n",
405+
"just because it's open_rasterio method works better with COGs.\n",
406+
"Instead of pulling down the full 30m resolution dataset,\n",
407+
"let's keep it easy on the servers and use overview_level 1."
408+
]
409+
},
410+
{
411+
"cell_type": "code",
412+
"execution_count": null,
413+
"metadata": {},
414+
"outputs": [],
415+
"source": [
416+
"band = {}\n",
417+
"for band_number in [4, 3, 2]: # Landsat 8's Red, Green and Blue bands\n",
418+
" band[band_number] = rioxarray.open_rasterio(\n",
419+
" f\"https://landsat-pds.s3.amazonaws.com/c1/L8/073/089/LC08_L1TP_073089_20191005_20191018_01_T1/LC08_L1TP_073089_20191005_20191018_01_T1_B{band_number}.TIF\",\n",
420+
" overview_level=1,\n",
421+
" masked=True,\n",
422+
" ).squeeze()\n",
423+
"band[3] # preview the green band"
424+
]
425+
},
426+
{
427+
"cell_type": "markdown",
428+
"metadata": {},
429+
"source": [
430+
"We'll concatenate all the bands into a grid to save us having to write for-loops later."
431+
]
432+
},
433+
{
434+
"cell_type": "code",
435+
"execution_count": null,
436+
"metadata": {},
437+
"outputs": [],
438+
"source": [
439+
"grid = xr.concat(\n",
440+
" objs=[band[4], band[3], band[2]],\n",
441+
" dim=pd.Index(data=[4, 3, 2], name=\"band\")\n",
442+
")"
443+
]
444+
},
445+
{
446+
"cell_type": "markdown",
447+
"metadata": {},
448+
"source": [
449+
"We can try to plot an RGB image by passing a list of the bands, but it won't turn out so good."
450+
]
451+
},
452+
{
453+
"cell_type": "code",
454+
"execution_count": null,
455+
"metadata": {},
456+
"outputs": [],
457+
"source": [
458+
"fig = pygmt.Figure()\n",
459+
"fig.grdimage(grid=[grid.sel(band=4), grid.sel(band=3), grid.sel(band=2)], frame=True)\n",
460+
"fig.show()"
461+
]
462+
},
463+
{
464+
"cell_type": "markdown",
465+
"metadata": {},
466+
"source": [
467+
"### **Scaling to 8 bits (0-255)**\n",
468+
"\n",
469+
"Each band in a Landsat 8 GeoTIFF file has a dynamic range of 16 bits (0-65535).\n",
470+
"However, in order to plot an RGB (red, green, blue) image in GMT,\n",
471+
"each individual band needs to have their values scaled to between 0-255 (i.e. 8 bits).\n",
472+
"\n",
473+
"GMT does have a [grdhisteq](https://docs.generic-mapping-tools.org/latest/grdhisteq.html) function\n",
474+
"for histogram equalization,\n",
475+
"but it is not available in PyGMT yet.\n",
476+
"We'll use scikit-image's\n",
477+
"[equalize_hist](https://scikit-image.org/docs/dev/api/skimage.exposure#skimage.exposure.equalize_hist)\n",
478+
"function instead for now."
479+
]
480+
},
481+
{
482+
"cell_type": "code",
483+
"execution_count": null,
484+
"metadata": {},
485+
"outputs": [],
486+
"source": [
487+
"grid.values = 256 * skimage.exposure.equalize_hist(\n",
488+
" image=grid.values,\n",
489+
" nbins=2**16, # 16 bit of Landsat 8\n",
490+
" mask=(grid != 0).values, # don't include zero values in histogram stretch\n",
491+
")"
492+
]
493+
},
494+
{
495+
"cell_type": "markdown",
496+
"metadata": {},
497+
"source": [
498+
"### **Reproject and Plot**"
499+
]
500+
},
501+
{
502+
"cell_type": "markdown",
503+
"metadata": {},
504+
"source": [
505+
"Let's reproject our grid so that we can more easily reference coordinates in longitude and latitude later."
506+
]
507+
},
508+
{
509+
"cell_type": "code",
510+
"execution_count": null,
511+
"metadata": {},
512+
"outputs": [],
513+
"source": [
514+
"grid = grid.rio.reproject(\n",
515+
" \"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\"\n",
516+
")"
517+
]
518+
},
519+
{
520+
"cell_type": "markdown",
521+
"metadata": {},
522+
"source": [
523+
"Now the image should look much nicer."
524+
]
525+
},
526+
{
527+
"cell_type": "code",
528+
"execution_count": null,
529+
"metadata": {},
530+
"outputs": [],
531+
"source": [
532+
"fig = pygmt.Figure()\n",
533+
"fig.grdimage(grid=[grid.sel(band=4), grid.sel(band=3), grid.sel(band=2)], frame=True)\n",
534+
"fig.show()"
535+
]
536+
},
537+
{
538+
"cell_type": "markdown",
539+
"metadata": {},
540+
"source": [
541+
"We can also subset our map to that of Wellington using the region argument as usual."
542+
]
543+
},
544+
{
545+
"cell_type": "code",
546+
"execution_count": null,
547+
"metadata": {},
548+
"outputs": [],
549+
"source": [
550+
"fig = pygmt.Figure()\n",
551+
"fig.grdimage(\n",
552+
" grid=[grid.sel(band=4), grid.sel(band=3), grid.sel(band=2)],\n",
553+
" region=[174.6, 174.9, -41.4, -41.1],\n",
554+
" frame=True\n",
555+
")\n",
556+
"fig.show()"
557+
]
558+
},
559+
{
560+
"cell_type": "markdown",
561+
"metadata": {},
562+
"source": [
563+
"## **2.4 Exploration time - Points or Grids, your pick**\n",
564+
"\n",
565+
"It's time to make a map on your own again!\n",
566+
"Depending on what you're interested in,\n",
567+
"use one of the Search APIs linked above\n",
568+
"(either the [GeoNet Quake Search](https://quakesearch.geonet.org.nz)\n",
569+
"or [Remote Pixel's Satellite Search API](https://search.remotepixel.ca))\n",
570+
"to pull in a good chunk of data for your city, town or region.\n",
571+
"\n",
572+
"If you're plotting earthquake points,\n",
573+
"try to play with different colormaps or change the size of the symbols.\n",
574+
"If you're plotting Landsat 8 grids,\n",
575+
"try to plot a false colour map instead or calculate\n",
576+
"[NDVI](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index) values\n",
577+
"and plot it with a colorbar.\n",
578+
"\n",
579+
"Go!"
580+
]
581+
},
582+
{
583+
"cell_type": "code",
584+
"execution_count": null,
585+
"metadata": {},
586+
"outputs": [],
587+
"source": []
588+
},
589+
{
590+
"cell_type": "code",
591+
"execution_count": null,
592+
"metadata": {},
593+
"outputs": [],
594+
"source": []
595+
},
596+
{
597+
"cell_type": "code",
598+
"execution_count": null,
599+
"metadata": {},
600+
"outputs": [],
601+
"source": []
602+
},
362603
{
363604
"cell_type": "markdown",
364605
"metadata": {},

0 commit comments

Comments
 (0)