Accessing via Python¶
The following is a collection of Python examples demonstrating how to connect to GIBS access points and exercise various capabilities. Included are examples of how to visualize raster and vector-based data from GIBS, plot imagery on maps, list GIBS capabilities, access GIBS metadata, basic image analysis and more. Please scroll down or use the navigation bar to browse through the examples.
These examples are also downloadable as a Jupyter Notebook.Import Python Packages And Modules
Major packages are requests, xml, json, skiimage, matplotlib, cartopy and pillow image.
In [1]:
Copied!
# install necessary packages for imports
%pip install scikit-image
%pip install scikit-learn
%pip install matplotlib
%pip install cartopy
%pip install folium
%pip install mapbox_vector_tile
%pip install lxml
%pip install pandas
%pip install owslib
%pip install geopandas
%pip install rasterio
%pip install fiona
%pip install ipyleaflet
%pip install cairosvg # If needed, more specific install instructions for cairosvg: https://cairosvg.org/documentation/
# install necessary packages for imports
%pip install scikit-image
%pip install scikit-learn
%pip install matplotlib
%pip install cartopy
%pip install folium
%pip install mapbox_vector_tile
%pip install lxml
%pip install pandas
%pip install owslib
%pip install geopandas
%pip install rasterio
%pip install fiona
%pip install ipyleaflet
%pip install cairosvg # If needed, more specific install instructions for cairosvg: https://cairosvg.org/documentation/
In [2]:
Copied!
import os
from io import BytesIO
from skimage import io
import requests
import json
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import folium
import urllib.request
import urllib.parse
import mapbox_vector_tile
import xml.etree.ElementTree as xmlet
import lxml.etree as xmltree
from PIL import Image as plimg
from PIL import ImageDraw
import numpy as np
import pandas as pd
from owslib.wms import WebMapService
from IPython.display import Image, display
import geopandas as gpd
from shapely.geometry import box
import urllib.request
import rasterio
from rasterio.mask import mask
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.plot import show
import fiona
from datetime import datetime, timedelta
%matplotlib inline
import os
from io import BytesIO
from skimage import io
import requests
import json
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import folium
import urllib.request
import urllib.parse
import mapbox_vector_tile
import xml.etree.ElementTree as xmlet
import lxml.etree as xmltree
from PIL import Image as plimg
from PIL import ImageDraw
import numpy as np
import pandas as pd
from owslib.wms import WebMapService
from IPython.display import Image, display
import geopandas as gpd
from shapely.geometry import box
import urllib.request
import rasterio
from rasterio.mask import mask
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.plot import show
import fiona
from datetime import datetime, timedelta
%matplotlib inline
OGC Web Map Service (WMS)¶
Web Map Service (WMS) is the preferred method for accessing static imagery (whereas Web Map Tile Service WMTS is preferred for interactive web maps). For smaller-scale, single image requests, WMS is usually easier to configure than WMTS and can also perform server-side compositing of multiple layers (both vector and raster).
Basic WMS Connection¶
First we will connect to the GIBS WMS Service and visualize the MODIS_Terra_CorrectedReflectance_TrueColor layer.
In [2]:
Copied!
# Connect to GIBS WMS Service
wms = WebMapService('https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?', version='1.1.1')
# Configure request for MODIS_Terra_CorrectedReflectance_TrueColor
img = wms.getmap(layers=['MODIS_Terra_CorrectedReflectance_TrueColor'], # Layers
srs='epsg:4326', # Map projection
bbox=(-180,-90,180,90), # Bounds
size=(1200, 600), # Image size
time='2021-09-21', # Time of data
format='image/png', # Image format
transparent=True) # Nodata transparency
# Save output PNG to a file
out = open('python-examples/MODIS_Terra_CorrectedReflectance_TrueColor.png', 'wb')
out.write(img.read())
out.close()
# View image
Image('python-examples/MODIS_Terra_CorrectedReflectance_TrueColor.png')
# Connect to GIBS WMS Service
wms = WebMapService('https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?', version='1.1.1')
# Configure request for MODIS_Terra_CorrectedReflectance_TrueColor
img = wms.getmap(layers=['MODIS_Terra_CorrectedReflectance_TrueColor'], # Layers
srs='epsg:4326', # Map projection
bbox=(-180,-90,180,90), # Bounds
size=(1200, 600), # Image size
time='2021-09-21', # Time of data
format='image/png', # Image format
transparent=True) # Nodata transparency
# Save output PNG to a file
out = open('python-examples/MODIS_Terra_CorrectedReflectance_TrueColor.png', 'wb')
out.write(img.read())
out.close()
# View image
Image('python-examples/MODIS_Terra_CorrectedReflectance_TrueColor.png')
Out[2]:
Get WMS Capabilities¶
For WMS, first we want to access the "GetCapabilities" document . GIBS provides four map projections, so there are four WMS endpoints GetCapabilities:
Geographic - EPSG:4326: https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi
Web Mercator - EPSG:3857: https://gibs.earthdata.nasa.gov/wms/epsg3857/best/wms.cgi
Arctic polar stereographic - EPSG:3413: https://gibs.earthdata.nasa.gov/wms/epsg3413/best/wms.cgi
Antarctic polar stereographic - EPSG:3031: https://gibs.earthdata.nasa.gov/wms/epsg3031/best/wms.cgi The code below will show how to get capabilities.
In [3]:
Copied!
# Construct capability URL.
wmsUrl = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
SERVICE=WMS&REQUEST=GetCapabilities'
# Request WMS capabilities.
response = requests.get(wmsUrl)
# Display capabilities XML in original format. Tag and content in one line.
WmsXml = xmltree.fromstring(response.content)
# print(xmltree.tostring(WmsXml, pretty_print = True, encoding = str))
# Construct capability URL.
wmsUrl = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
SERVICE=WMS&REQUEST=GetCapabilities'
# Request WMS capabilities.
response = requests.get(wmsUrl)
# Display capabilities XML in original format. Tag and content in one line.
WmsXml = xmltree.fromstring(response.content)
# print(xmltree.tostring(WmsXml, pretty_print = True, encoding = str))
Display WMS All Layers¶
Parse WMS capabilities XML to get total number of layers and display all layer names.
In [4]:
Copied!
# Currently total layers are 1081.
# Coverts response to XML tree.
WmsTree = xmlet.fromstring(response.content)
alllayer = []
layerNumber = 0
# Parse XML.
for child in WmsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wms}Capability/{http://www.opengis.net/wms}Layer//*/"):
if layer.tag == '{http://www.opengis.net/wms}Layer':
f = layer.find("{http://www.opengis.net/wms}Name")
if f is not None:
alllayer.append(f.text)
layerNumber += 1
print('There are layers: ' + str(layerNumber))
for one in sorted(alllayer)[:5]:
print(one)
print('...')
for one in sorted(alllayer)[-5:]:
print(one)
# Currently total layers are 1081.
# Coverts response to XML tree.
WmsTree = xmlet.fromstring(response.content)
alllayer = []
layerNumber = 0
# Parse XML.
for child in WmsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wms}Capability/{http://www.opengis.net/wms}Layer//*/"):
if layer.tag == '{http://www.opengis.net/wms}Layer':
f = layer.find("{http://www.opengis.net/wms}Name")
if f is not None:
alllayer.append(f.text)
layerNumber += 1
print('There are layers: ' + str(layerNumber))
for one in sorted(alllayer)[:5]:
print(one)
print('...')
for one in sorted(alllayer)[-5:]:
print(one)
There are layers: 1175 AIRS_L2_Carbon_Monoxide_500hPa_Volume_Mixing_Ratio_Day AIRS_L2_Carbon_Monoxide_500hPa_Volume_Mixing_Ratio_Night AIRS_L2_Cloud_Top_Height_Day AIRS_L2_Cloud_Top_Height_Night AIRS_L2_Dust_Score_Day ... VIIRS_SNPP_L2_Sea_Surface_Temp_Day VIIRS_SNPP_L2_Sea_Surface_Temp_Night VIIRS_SNPP_Thermal_Anomalies_375m_All VIIRS_SNPP_Thermal_Anomalies_375m_Day VIIRS_SNPP_Thermal_Anomalies_375m_Night
Search WMS Layer And Its Attributes¶
Requesting WMS data needs layer name, bounding box, time, projection, data format and so on. Enter a layer name to search its attributes.
In [5]:
Copied!
# Define layername to use.
layerName = 'Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual'
# Get general information of WMS.
for child in WmsTree.iter():
if child.tag == '{http://www.opengis.net/wms}WMS_Capabilities':
print('Version: ' +child.get('version'))
if child.tag == '{http://www.opengis.net/wms}Service':
print('Service: ' +child.find("{http://www.opengis.net/wms}Name").text)
if child.tag == '{http://www.opengis.net/wms}Request':
print('Request: ')
for e in child:
print('\t ' + e.tag.partition('}')[2])
all = child.findall(".//{http://www.opengis.net/wms}Format")
if all is not None:
print("Format: ")
for g in all:
print("\t " + g.text)
for e in child.iter():
if e.tag == "{http://www.opengis.net/wms}OnlineResource":
print('URL: ' + e.get('{http://www.w3.org/1999/xlink}href'))
break
# Get layer attributes.
for child in WmsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wms}Capability/{http://www.opengis.net/wms}Layer//*/"):
if layer.tag == '{http://www.opengis.net/wms}Layer':
f = layer.find("{http://www.opengis.net/wms}Name")
if f is not None:
if f.text == layerName:
# Layer name.
print('Layer: ' + f.text)
# All elements and attributes:
# CRS
e = layer.find("{http://www.opengis.net/wms}CRS")
if e is not None:
print('\t CRS: ' + e.text)
# BoundingBox.
e = layer.find("{http://www.opengis.net/wms}EX_GeographicBoundingBox")
if e is not None:
print('\t LonMin: ' + e.find("{http://www.opengis.net/wms}westBoundLongitude").text)
print('\t LonMax: ' + e.find("{http://www.opengis.net/wms}eastBoundLongitude").text)
print('\t LatMin: ' + e.find("{http://www.opengis.net/wms}southBoundLatitude").text)
print('\t LatMax: ' + e.find("{http://www.opengis.net/wms}northBoundLatitude").text)
# Time extent.
e = layer.find("{http://www.opengis.net/wms}Dimension")
if e is not None:
print('\t TimeExtent: ' + e.text)
# Style.
e = layer.find("{http://www.opengis.net/wms}Style")
if e is not None:
f = e.find("{http://www.opengis.net/wms}Name")
if f is not None:
print('\t Style: ' + f.text)
print('')
# Define layername to use.
layerName = 'Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual'
# Get general information of WMS.
for child in WmsTree.iter():
if child.tag == '{http://www.opengis.net/wms}WMS_Capabilities':
print('Version: ' +child.get('version'))
if child.tag == '{http://www.opengis.net/wms}Service':
print('Service: ' +child.find("{http://www.opengis.net/wms}Name").text)
if child.tag == '{http://www.opengis.net/wms}Request':
print('Request: ')
for e in child:
print('\t ' + e.tag.partition('}')[2])
all = child.findall(".//{http://www.opengis.net/wms}Format")
if all is not None:
print("Format: ")
for g in all:
print("\t " + g.text)
for e in child.iter():
if e.tag == "{http://www.opengis.net/wms}OnlineResource":
print('URL: ' + e.get('{http://www.w3.org/1999/xlink}href'))
break
# Get layer attributes.
for child in WmsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wms}Capability/{http://www.opengis.net/wms}Layer//*/"):
if layer.tag == '{http://www.opengis.net/wms}Layer':
f = layer.find("{http://www.opengis.net/wms}Name")
if f is not None:
if f.text == layerName:
# Layer name.
print('Layer: ' + f.text)
# All elements and attributes:
# CRS
e = layer.find("{http://www.opengis.net/wms}CRS")
if e is not None:
print('\t CRS: ' + e.text)
# BoundingBox.
e = layer.find("{http://www.opengis.net/wms}EX_GeographicBoundingBox")
if e is not None:
print('\t LonMin: ' + e.find("{http://www.opengis.net/wms}westBoundLongitude").text)
print('\t LonMax: ' + e.find("{http://www.opengis.net/wms}eastBoundLongitude").text)
print('\t LatMin: ' + e.find("{http://www.opengis.net/wms}southBoundLatitude").text)
print('\t LatMax: ' + e.find("{http://www.opengis.net/wms}northBoundLatitude").text)
# Time extent.
e = layer.find("{http://www.opengis.net/wms}Dimension")
if e is not None:
print('\t TimeExtent: ' + e.text)
# Style.
e = layer.find("{http://www.opengis.net/wms}Style")
if e is not None:
f = e.find("{http://www.opengis.net/wms}Name")
if f is not None:
print('\t Style: ' + f.text)
print('')
Version: 1.3.0 Service: WMS Request: GetCapabilities GetMap Format: text/xml image/png application/vnd.google-earth.kml.xml application/vnd.google-earth.kmz image/jpeg image/png; mode=8bit image/vnd.jpeg-png image/vnd.jpeg-png8 image/tiff application/json URL: https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi? Layer: Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual CRS: EPSG:4326 LonMin: -180 LonMax: 180 LatMin: -90 LatMax: 90 TimeExtent: 1983-12-01/1985-12-01/P1Y,1988-12-01/1990-12-01/P1Y,1998-12-01/2000-12-01/P1Y
You can also use rasterio to access the properties of a geospatial raster file
In [2]:
Copied!
# Save a global extents tiff file
# Connect to GIBS WMS Service
wms = WebMapService('https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?', version='1.3.0')
# Configure request for MODIS_Terra_SurfaceReflectance_Bands143
img = wms.getmap(layers=['MODIS_Terra_SurfaceReflectance_Bands143'], # Layers
srs='epsg:4326', # Map projection
bbox=(-180,-90,180,90), # Bounds
size=(1024, 512), # Image size
time='2021-11-25', # Time of data
format='image/tiff', # Image format
transparent=True) # Nodata transparency
# Save output TIFF to a file
out = open('global.tiff', 'wb')
out.write(img.read())
out.close()
# Access properties of geospatial raster file
with rasterio.open("global.tiff") as src:
print(src.width, src.height)
print(src.crs)
print(src.transform)
print(src.count)
print(src.indexes)
# Save a global extents tiff file
# Connect to GIBS WMS Service
wms = WebMapService('https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?', version='1.3.0')
# Configure request for MODIS_Terra_SurfaceReflectance_Bands143
img = wms.getmap(layers=['MODIS_Terra_SurfaceReflectance_Bands143'], # Layers
srs='epsg:4326', # Map projection
bbox=(-180,-90,180,90), # Bounds
size=(1024, 512), # Image size
time='2021-11-25', # Time of data
format='image/tiff', # Image format
transparent=True) # Nodata transparency
# Save output TIFF to a file
out = open('global.tiff', 'wb')
out.write(img.read())
out.close()
# Access properties of geospatial raster file
with rasterio.open("global.tiff") as src:
print(src.width, src.height)
print(src.crs)
print(src.transform)
print(src.count)
print(src.indexes)
1024 512 EPSG:4326 | 0.35, 0.00,-180.00| | 0.00,-0.35, 90.00| | 0.00, 0.00, 1.00| 4 (1, 2, 3, 4)
Visualize WMS Raster Data In Geographic Projection¶
This example shows how to get a geographic projection (EPSG:4326) image. Use a layer name and its attributes to form a URL for requesting a WMS image. After an image is returned, display it on a map using matplotlib and cartopy.
In [6]:
Copied!
# Construct Geographic projection URL.
proj4326 = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
version=1.3.0&service=WMS&\
request=GetMap&format=image/png&STYLE=default&bbox=-40,-40,40,40&CRS=EPSG:4326&\
HEIGHT=600&WIDTH=600&TIME=2000-12-01&layers=Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual'
# Request image.
img = io.imread(proj4326)
# Display image on map.
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
extent = (-40, 40, -40, 40)
plt.imshow(img, transform = ccrs.PlateCarree(), extent = extent, origin = 'upper')
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 30, -30, 0])
gl.ylocator = mticker.FixedLocator([-30, 0, 30])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
plt.title('WMS Landsat Reflectance In Geographic Projection',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
# Construct Geographic projection URL.
proj4326 = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
version=1.3.0&service=WMS&\
request=GetMap&format=image/png&STYLE=default&bbox=-40,-40,40,40&CRS=EPSG:4326&\
HEIGHT=600&WIDTH=600&TIME=2000-12-01&layers=Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual'
# Request image.
img = io.imread(proj4326)
# Display image on map.
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
extent = (-40, 40, -40, 40)
plt.imshow(img, transform = ccrs.PlateCarree(), extent = extent, origin = 'upper')
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 30, -30, 0])
gl.ylocator = mticker.FixedLocator([-30, 0, 30])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
plt.title('WMS Landsat Reflectance In Geographic Projection',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
Visualize WMS Raster Data In Web Mercator Projection¶
This example shows how to get an image in WMS Web Mercator projection (EPSG:3857) and display it on map.
In [7]:
Copied!
# Construct Web Mercator projection URL.
proj3857 = 'https://gibs.earthdata.nasa.gov/wms/epsg3857/best/wms.cgi?\
version=1.3.0&service=WMS&\
request=GetMap&format=image/png&STYLE=default&bbox=-8000000,-8000000,8000000,8000000&\
CRS=EPSG:3857&HEIGHT=600&WIDTH=600&TIME=2000-12-01&layers=Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual'
# Request image.
img=io.imread(proj3857)
# Display image on map.
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection = ccrs.Mercator.GOOGLE)
extent = (-8000000, 8000000, -8000000, 8000000)
plt.imshow(img, transform = ccrs.Mercator.GOOGLE, extent = extent, origin = 'upper')
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 30, -30, 0])
gl.ylocator = mticker.FixedLocator([-30, 0, 30])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
plt.title('WMS Landsat Reflectance In Web Mercator Projection',
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
# Construct Web Mercator projection URL.
proj3857 = 'https://gibs.earthdata.nasa.gov/wms/epsg3857/best/wms.cgi?\
version=1.3.0&service=WMS&\
request=GetMap&format=image/png&STYLE=default&bbox=-8000000,-8000000,8000000,8000000&\
CRS=EPSG:3857&HEIGHT=600&WIDTH=600&TIME=2000-12-01&layers=Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual'
# Request image.
img=io.imread(proj3857)
# Display image on map.
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection = ccrs.Mercator.GOOGLE)
extent = (-8000000, 8000000, -8000000, 8000000)
plt.imshow(img, transform = ccrs.Mercator.GOOGLE, extent = extent, origin = 'upper')
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 30, -30, 0])
gl.ylocator = mticker.FixedLocator([-30, 0, 30])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
plt.title('WMS Landsat Reflectance In Web Mercator Projection',
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
Visualize WMS Raster Data In Arctic Polar Stereographic Projection¶
This example shows how to get WMS Arctic Polar Stereographic projection (EPSG:3413) image and display it on map.
In [8]:
Copied!
# Construct Arctic Polar Stereographic projection URL.
proj3413 = 'https://gibs.earthdata.nasa.gov/wms/epsg3413/best/wms.cgi?\
version=1.3.0&service=WMS&request=GetMap&\
format=image/png&STYLE=default&bbox=-4194300,-4194300,4194300,4194300&CRS=EPSG:3413&\
HEIGHT=512&WIDTH=512&TIME=2021-08-01&layers=MODIS_Terra_CorrectedReflectance_TrueColor'
# Request image.
img = io.imread(proj3413)
# Display image on map.
plt.figure(figsize=(5, 5))
ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=-45))
plt.imshow(img, extent = (-4194300,4194300,-4194300,4194300), origin = 'upper')
# Draw coastline and grid.
ax.coastlines(color='blue', linewidth=1)
ax.gridlines()
plt.title('WMS Terra True Color Image In Arctic Polar Stereographic',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
# Construct Arctic Polar Stereographic projection URL.
proj3413 = 'https://gibs.earthdata.nasa.gov/wms/epsg3413/best/wms.cgi?\
version=1.3.0&service=WMS&request=GetMap&\
format=image/png&STYLE=default&bbox=-4194300,-4194300,4194300,4194300&CRS=EPSG:3413&\
HEIGHT=512&WIDTH=512&TIME=2021-08-01&layers=MODIS_Terra_CorrectedReflectance_TrueColor'
# Request image.
img = io.imread(proj3413)
# Display image on map.
plt.figure(figsize=(5, 5))
ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=-45))
plt.imshow(img, extent = (-4194300,4194300,-4194300,4194300), origin = 'upper')
# Draw coastline and grid.
ax.coastlines(color='blue', linewidth=1)
ax.gridlines()
plt.title('WMS Terra True Color Image In Arctic Polar Stereographic',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
Visualize WMS Raster Data In Antarctic Polar Stereographic Projection¶
This example shows how to get WMS Antarctic Polar Stereographic projection (EPSG:3031) image and display it on map.
In [9]:
Copied!
# Construct Antarctic Polar Stereographic project
proj3031 = 'https://gibs.earthdata.nasa.gov/wms/epsg3031/best/wms.cgi?\
version=1.3.0&service=WMS&request=GetMap&\
format=image/png&STYLE=default&bbox=-4194300,-4194300,4194300,4194300&CRS=EPSG:3031&\
HEIGHT=512&WIDTH=512&TIME=2021-03-01&layers=MODIS_Terra_CorrectedReflectance_TrueColor'
# Request image.
img = io.imread(proj3031)
# Display image on map.
plt.figure(figsize=(5, 5))
ax = plt.axes(projection=ccrs.SouthPolarStereo())
plt.imshow(img, extent = (-4194300,4194300,-4194300,4194300), origin = 'upper')
# Draw coastline and grid.
ax.coastlines(color='blue', linewidth=1)
ax.gridlines()
plt.title('WMS Terra True Color Image In Antarctic Polar Stereographic',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
# Construct Antarctic Polar Stereographic project
proj3031 = 'https://gibs.earthdata.nasa.gov/wms/epsg3031/best/wms.cgi?\
version=1.3.0&service=WMS&request=GetMap&\
format=image/png&STYLE=default&bbox=-4194300,-4194300,4194300,4194300&CRS=EPSG:3031&\
HEIGHT=512&WIDTH=512&TIME=2021-03-01&layers=MODIS_Terra_CorrectedReflectance_TrueColor'
# Request image.
img = io.imread(proj3031)
# Display image on map.
plt.figure(figsize=(5, 5))
ax = plt.axes(projection=ccrs.SouthPolarStereo())
plt.imshow(img, extent = (-4194300,4194300,-4194300,4194300), origin = 'upper')
# Draw coastline and grid.
ax.coastlines(color='blue', linewidth=1)
ax.gridlines()
plt.title('WMS Terra True Color Image In Antarctic Polar Stereographic',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
Visualize WMS Raster Data In Any Projection¶
This example shows how to reproject WMS Raster Data to the projection of your choice using rasterio
In [3]:
Copied!
# Here you can set the new projection to whatever you like
dst_crs = 'EPSG:6933'
# Save a global extents tiff file
# Connect to GIBS WMS Service
wms = WebMapService('https://gibs.earthdata.nasa.gov/wms/epsg3857/best/wms.cgi?', version='1.3.0')
# Configure request for Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual
img = wms.getmap(layers=['Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual'], # Layers
srs='epsg:3857', # Map projection
bbox=(-8000000,-8000000,8000000,8000000), # Bounds
size=(600, 600), # Image size
time='2000-12-01', # Time of data
format='image/tiff', # Image format
transparent=True) # Nodata transparency
# Save output TIFF to a file
out = open('global_extents_3857.tiff', 'wb')
out.write(img.read())
out.close()
with rasterio.open('global_extents_3857.tiff') as src:
transform, width, height = calculate_default_transform(
src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
with rasterio.open(os.getcwd()+'/reprojectedImage.byte.tif', 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest)
reprojected_image = plimg.open('reprojectedImage.byte.tif', 'r')
plt.imshow(reprojected_image)
# Here you can set the new projection to whatever you like
dst_crs = 'EPSG:6933'
# Save a global extents tiff file
# Connect to GIBS WMS Service
wms = WebMapService('https://gibs.earthdata.nasa.gov/wms/epsg3857/best/wms.cgi?', version='1.3.0')
# Configure request for Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual
img = wms.getmap(layers=['Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual'], # Layers
srs='epsg:3857', # Map projection
bbox=(-8000000,-8000000,8000000,8000000), # Bounds
size=(600, 600), # Image size
time='2000-12-01', # Time of data
format='image/tiff', # Image format
transparent=True) # Nodata transparency
# Save output TIFF to a file
out = open('global_extents_3857.tiff', 'wb')
out.write(img.read())
out.close()
with rasterio.open('global_extents_3857.tiff') as src:
transform, width, height = calculate_default_transform(
src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
with rasterio.open(os.getcwd()+'/reprojectedImage.byte.tif', 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest)
reprojected_image = plimg.open('reprojectedImage.byte.tif', 'r')
plt.imshow(reprojected_image)
Out[3]:
<matplotlib.image.AxesImage at 0x1297e9940>
Visualize WMS Global Raster Data¶
This example shows how to get WMS global image and to display it on map.
In [10]:
Copied!
# Construct global image URL.
proj4326 = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
version=1.3.0&service=WMS&request=GetMap&\
format=image/jpeg&STYLE=default&bbox=-90,-180,90,180&CRS=EPSG:4326&\
HEIGHT=512&WIDTH=512&TIME=2021-11-25&layers=MODIS_Terra_SurfaceReflectance_Bands143'
# Request image.
img = io.imread(proj4326)
# Display image on map.
plt.figure(figsize = (9, 6))
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 0))
cmp = plt.imshow(img, transform = ccrs.PlateCarree(), extent = (-180,180,-90,90), origin = 'upper')
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 60, 120, 180, -120, -60, 0])
gl.ylocator = mticker.FixedLocator([-90, -60, -30, 0, 30, 60, 90])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
# Draw coastline.
ax.coastlines()
plt.title('WMS Terra MODIS Surface Reflectance',fontname="Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
# Construct global image URL.
proj4326 = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
version=1.3.0&service=WMS&request=GetMap&\
format=image/jpeg&STYLE=default&bbox=-90,-180,90,180&CRS=EPSG:4326&\
HEIGHT=512&WIDTH=512&TIME=2021-11-25&layers=MODIS_Terra_SurfaceReflectance_Bands143'
# Request image.
img = io.imread(proj4326)
# Display image on map.
plt.figure(figsize = (9, 6))
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 0))
cmp = plt.imshow(img, transform = ccrs.PlateCarree(), extent = (-180,180,-90,90), origin = 'upper')
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 60, 120, 180, -120, -60, 0])
gl.ylocator = mticker.FixedLocator([-90, -60, -30, 0, 30, 60, 90])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
# Draw coastline.
ax.coastlines()
plt.title('WMS Terra MODIS Surface Reflectance',fontname="Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
Visualize WMS Global Raster Data with Folium¶
This example shows how to display a WMS global image on an OpenStreetMap using Folium.
In [4]:
Copied!
# Other tile options are "CartoDB Positron" and "CartoDB Voyager"
m = folium.Map(location=[41, -70], zoom_start=4, tiles="OpenStreetMap")
folium.WmsTileLayer(
url="https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi",
name="WMS Landsat Reflectance",
fmt="image/png",
layers="Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual",
transparent=True,
overlay=True,
control=True,
).add_to(m)
folium.LayerControl().add_to(m)
m
# Other tile options are "CartoDB Positron" and "CartoDB Voyager"
m = folium.Map(location=[41, -70], zoom_start=4, tiles="OpenStreetMap")
folium.WmsTileLayer(
url="https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi",
name="WMS Landsat Reflectance",
fmt="image/png",
layers="Landsat_WELD_CorrectedReflectance_Bands157_Global_Annual",
transparent=True,
overlay=True,
control=True,
).add_to(m)
folium.LayerControl().add_to(m)
m
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Visualize WMS Global Vector Data¶
This example shows how to get WMS global vector data and to display it on map.
In [11]:
Copied!
# Construct WMS global vector URL.
wmsVector = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
TIME=2020-10-01T00:00:00Z&\
LAYERS=VIIRS_NOAA20_Thermal_Anomalies_375m_All&REQUEST=GetMap&SERVICE=WMS&\
FORMAT=image/png&WIDTH=480&HEIGHT=240&VERSION=1.1.1&SRS=epsg:4326&BBOX=-180,-90,180,90&TRANSPARENT=TRUE'
# Request image.
img = io.imread(wmsVector)
# Setup map size, projection and background.
fig = plt.figure(figsize = (10, 6))
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 0))
ax.set_facecolor("white")
ax.stock_img()
ax.coastlines()
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 60, 120, -120, -60, 0])
gl.ylocator = mticker.FixedLocator([-90, -60, -30, 0, 30, 60, 90])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
# Display image on map.
extent = (-180, 180, -90, 90)
plt.imshow(img, extent = extent)
plt.title('WMS Vector Data VIIRS Thermal Anomalies',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
print('')
# Construct WMS global vector URL.
wmsVector = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
TIME=2020-10-01T00:00:00Z&\
LAYERS=VIIRS_NOAA20_Thermal_Anomalies_375m_All&REQUEST=GetMap&SERVICE=WMS&\
FORMAT=image/png&WIDTH=480&HEIGHT=240&VERSION=1.1.1&SRS=epsg:4326&BBOX=-180,-90,180,90&TRANSPARENT=TRUE'
# Request image.
img = io.imread(wmsVector)
# Setup map size, projection and background.
fig = plt.figure(figsize = (10, 6))
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 0))
ax.set_facecolor("white")
ax.stock_img()
ax.coastlines()
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 60, 120, -120, -60, 0])
gl.ylocator = mticker.FixedLocator([-90, -60, -30, 0, 30, 60, 90])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
# Display image on map.
extent = (-180, 180, -90, 90)
plt.imshow(img, extent = extent)
plt.title('WMS Vector Data VIIRS Thermal Anomalies',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
print('')
Interactive Web Map with WMS¶
The next example shows how to display VIIRS_NOAA20_Thermal_Anomalies_375m_All layer in an interactive web map (may require additional Python libraries).
In [12]:
Copied!
from ipyleaflet import Map, WMSLayer, basemaps
# Make a WMS connection to a map layer
wms_layer = WMSLayer(url='https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?',
layers='VIIRS_NOAA20_Thermal_Anomalies_375m_All',
format='image/png',
transparent=True)
# Define map properties and add the WMS layer from above on top of basemap
m = Map(basemap=basemaps.NASAGIBS.BlueMarble, center=(0, -0), zoom=3)
m.add_layer(wms_layer)
# Display interactive web map
m
from ipyleaflet import Map, WMSLayer, basemaps
# Make a WMS connection to a map layer
wms_layer = WMSLayer(url='https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?',
layers='VIIRS_NOAA20_Thermal_Anomalies_375m_All',
format='image/png',
transparent=True)
# Define map properties and add the WMS layer from above on top of basemap
m = Map(basemap=basemaps.NASAGIBS.BlueMarble, center=(0, -0), zoom=3)
m.add_layer(wms_layer)
# Display interactive web map
m
Map(center=[0, 0], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title', 'zoom_out_text'…
Note that this will not render in our documents page, but will work if you try in your own notebook
Animated Web Map with WMS¶
The next example shows how to display IMERG_Precipitation_Rate layer in an animated web map (may require additional Python libraries).
In [5]:
Copied!
wms = WebMapService('https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?', version='1.3.0')
layers = ['MODIS_Aqua_CorrectedReflectance_TrueColor',
'IMERG_Precipitation_Rate',
'Reference_Features',
'Reference_Labels']
color = 'rgb(255,255,255)'
frames = []
start_date = datetime(2022, 9, 25)
end_date = datetime(2022, 10, 1)
dates = pd.date_range(start_date,end_date-timedelta(days=1),freq='d')
for day in dates:
datatime = day.strftime("%Y-%m-%d")
img = wms.getmap(layers=layers, # Layers
srs='epsg:4326', # Map projection
bbox=(-87, 18, -72, 35), # Bounds
size=(600,600), # Image size
time=datatime, # Time of data
format='image/png', # Image format
transparent=True) # Nodata transparency
image = plimg.open(img)
draw = ImageDraw.Draw(image)
(x, y) = (50, 20)
draw.text((x, y), f'IMERG Precipitation Rate - {datatime}', fill=color)
frames.append(image)
frames[0].save('IMERG_Precipitation_Rate_Ian.gif',
format='GIF',
append_images=frames,
save_all=True,
duration=1000,
loop=0)
Image('IMERG_Precipitation_Rate_Ian.gif')
wms = WebMapService('https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?', version='1.3.0')
layers = ['MODIS_Aqua_CorrectedReflectance_TrueColor',
'IMERG_Precipitation_Rate',
'Reference_Features',
'Reference_Labels']
color = 'rgb(255,255,255)'
frames = []
start_date = datetime(2022, 9, 25)
end_date = datetime(2022, 10, 1)
dates = pd.date_range(start_date,end_date-timedelta(days=1),freq='d')
for day in dates:
datatime = day.strftime("%Y-%m-%d")
img = wms.getmap(layers=layers, # Layers
srs='epsg:4326', # Map projection
bbox=(-87, 18, -72, 35), # Bounds
size=(600,600), # Image size
time=datatime, # Time of data
format='image/png', # Image format
transparent=True) # Nodata transparency
image = plimg.open(img)
draw = ImageDraw.Draw(image)
(x, y) = (50, 20)
draw.text((x, y), f'IMERG Precipitation Rate - {datatime}', fill=color)
frames.append(image)
frames[0].save('IMERG_Precipitation_Rate_Ian.gif',
format='GIF',
append_images=frames,
save_all=True,
duration=1000,
loop=0)
Image('IMERG_Precipitation_Rate_Ian.gif')
Out[5]:
<IPython.core.display.Image object>
Note that this will not render in our documents page, but will work if you try in your own notebook
Display a Legend for a WMS Layer¶
This example shows how the WMS GetCapabilities XML tree can be used to find and display a legend associated with a particular layer. For this example, we will use the "Croplands (Global Agricultural Lands, 2000)" layer.
We will use WmsTree, which was previously created in the Display WMS All Layers example from the XML tree returned by the WMS GetCapabilities request.
In [13]:
Copied!
layerName = "Agricultural_Lands_Croplands_2000"
legendImg = None
for child in WmsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wms}Capability/{http://www.opengis.net/wms}Layer//*/"):
if layer.tag == '{http://www.opengis.net/wms}Layer':
f = layer.find("{http://www.opengis.net/wms}Name")
if f is not None:
if f.text == layerName:
# Style.
e = layer.find(("{http://www.opengis.net/wms}Style/" +
"{http://www.opengis.net/wms}LegendURL/" +
"{http://www.opengis.net/wms}OnlineResource"))
if e is not None:
legendURL = e.attrib["{http://www.w3.org/1999/xlink}href"]
legendImg = Image(url=legendURL)
print("Legend URL:", legendURL)
display(legendImg)
layerName = "Agricultural_Lands_Croplands_2000"
legendImg = None
for child in WmsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wms}Capability/{http://www.opengis.net/wms}Layer//*/"):
if layer.tag == '{http://www.opengis.net/wms}Layer':
f = layer.find("{http://www.opengis.net/wms}Name")
if f is not None:
if f.text == layerName:
# Style.
e = layer.find(("{http://www.opengis.net/wms}Style/" +
"{http://www.opengis.net/wms}LegendURL/" +
"{http://www.opengis.net/wms}OnlineResource"))
if e is not None:
legendURL = e.attrib["{http://www.w3.org/1999/xlink}href"]
legendImg = Image(url=legendURL)
print("Legend URL:", legendURL)
display(legendImg)
Legend URL: https://gibs.earthdata.nasa.gov/legends/Agricultural_Lands_Croplands_2000_H.png
Visualize WMS Raster Data with a Legend¶
This example visualizes a WMS raster layer with its associated legend. It follows the procedure for visualizing the WMS layer established in the Visualize WMS Global Raster Data example and makes use of the legend URL obtained in the previous example (legendURL).
In [14]:
Copied!
# Construct global image URL.
proj4326 = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
version=1.3.0&service=WMS&request=GetMap&\
format=image/jpeg&STYLE=default&bbox=-90,-180,90,180&CRS=EPSG:4326&\
HEIGHT=512&WIDTH=512&TIME=2021-11-25&layers=Agricultural_Lands_Croplands_2000'
# Request image.
img = io.imread(proj4326)
# Display image on map.
plt.figure(figsize = (9, 6))
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 0))
cmp = plt.imshow(img, transform = ccrs.PlateCarree(), extent = (-180,180,-90,90), origin = 'upper')
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 60, 120, 180, -120, -60, 0])
gl.ylocator = mticker.FixedLocator([-90, -60, -30, 0, 30, 60, 90])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
# Draw coastline.
ax.coastlines()
plt.title('Croplands (Global Agricultural Lands, 2000)',fontname="Times New Roman", fontsize = 20, color = 'green')
# Get the legend image from the URL as a numpy array
legendImgArr = np.array(plimg.open(urllib.request.urlopen(legendURL)))
# use data coordinates to specify the position and dimensions of new inset axes
axin = ax.inset_axes([-125,-260,250,250],transform=ax.transData)
axin.imshow(legendImgArr)
axin.axis('off')
plt.show()
print('')
# Construct global image URL.
proj4326 = 'https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi?\
version=1.3.0&service=WMS&request=GetMap&\
format=image/jpeg&STYLE=default&bbox=-90,-180,90,180&CRS=EPSG:4326&\
HEIGHT=512&WIDTH=512&TIME=2021-11-25&layers=Agricultural_Lands_Croplands_2000'
# Request image.
img = io.imread(proj4326)
# Display image on map.
plt.figure(figsize = (9, 6))
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 0))
cmp = plt.imshow(img, transform = ccrs.PlateCarree(), extent = (-180,180,-90,90), origin = 'upper')
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth = 1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 60, 120, 180, -120, -60, 0])
gl.ylocator = mticker.FixedLocator([-90, -60, -30, 0, 30, 60, 90])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
# Draw coastline.
ax.coastlines()
plt.title('Croplands (Global Agricultural Lands, 2000)',fontname="Times New Roman", fontsize = 20, color = 'green')
# Get the legend image from the URL as a numpy array
legendImgArr = np.array(plimg.open(urllib.request.urlopen(legendURL)))
# use data coordinates to specify the position and dimensions of new inset axes
axin = ax.inset_axes([-125,-260,250,250],transform=ax.transData)
axin.imshow(legendImgArr)
axin.axis('off')
plt.show()
print('')
OGC Web Map Tile Service (WMTS)¶
Web Map Tile Service (WMTS) is normally used for interactive web mapping, but may be used for general visualizations and data analysis. WMTS is much more responsive for interactive maps and very scalable for generating large images or bulk downloads, but compared to WMS, it is more challenging to configure if you just need a single, reasonably-sized image.
Get WMTS Capabilities¶
This example shows how to get WMTS capabilities and display the GetCapabilities XML content.
In [15]:
Copied!
# Construct WMTS capability URL.
wmtsUrl = 'http://gibs.earthdata.nasa.gov/wmts/epsg4326/best/wmts.cgi?SERVICE=WMTS&REQUEST=GetCapabilities'
# Request capabilities.
response = requests.get(wmtsUrl)
# Display capability XML.
WmtsXml = xmltree.fromstring(response.content)
# Uncomment the following to display the large file:
# print(xmltree.tostring(WmtsXml, pretty_print = True, encoding = str))
# Construct WMTS capability URL.
wmtsUrl = 'http://gibs.earthdata.nasa.gov/wmts/epsg4326/best/wmts.cgi?SERVICE=WMTS&REQUEST=GetCapabilities'
# Request capabilities.
response = requests.get(wmtsUrl)
# Display capability XML.
WmtsXml = xmltree.fromstring(response.content)
# Uncomment the following to display the large file:
# print(xmltree.tostring(WmtsXml, pretty_print = True, encoding = str))
Display All Layers of WMTS Capabilities.¶
This example shows how to parse the WMTS GetCapabilities document and print the names of all of its layers.
In [16]:
Copied!
# Convert capability response to XML tree.
WmtsTree = xmlet.fromstring(response.content)
alllayer = []
layerNumber = 0
# Parse capability XML tree.
for child in WmtsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wmts/1.0}Layer"):
if '{http://www.opengis.net/wmts/1.0}Layer' == layer.tag:
f=layer.find("{http://www.opengis.net/ows/1.1}Identifier")
if f is not None:
alllayer.append(f.text)
layerNumber += 1
# Print the first five and last five layers.
print('Number of layers: ', layerNumber)
for one in sorted(alllayer)[:5]:
print(one)
print('...')
for one in sorted(alllayer)[-5:]:
print(one)
# Convert capability response to XML tree.
WmtsTree = xmlet.fromstring(response.content)
alllayer = []
layerNumber = 0
# Parse capability XML tree.
for child in WmtsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wmts/1.0}Layer"):
if '{http://www.opengis.net/wmts/1.0}Layer' == layer.tag:
f=layer.find("{http://www.opengis.net/ows/1.1}Identifier")
if f is not None:
alllayer.append(f.text)
layerNumber += 1
# Print the first five and last five layers.
print('Number of layers: ', layerNumber)
for one in sorted(alllayer)[:5]:
print(one)
print('...')
for one in sorted(alllayer)[-5:]:
print(one)
Number of layers: 1056 AIRS_L2_Carbon_Monoxide_500hPa_Volume_Mixing_Ratio_Day AIRS_L2_Carbon_Monoxide_500hPa_Volume_Mixing_Ratio_Night AIRS_L2_Cloud_Top_Height_Day AIRS_L2_Cloud_Top_Height_Night AIRS_L2_Dust_Score_Day ... VIIRS_SNPP_L2_Sea_Surface_Temp_Day VIIRS_SNPP_L2_Sea_Surface_Temp_Night VIIRS_SNPP_Thermal_Anomalies_375m_All VIIRS_SNPP_Thermal_Anomalies_375m_Day VIIRS_SNPP_Thermal_Anomalies_375m_Night
Search WMTS Vector Layer, Attributes And Vector Information¶
This example shows how to search a WMTS layer and to parse its attributes and vector information.
In [17]:
Copied!
# Get general information of WMTS from XML tree.
for child in WmtsTree.iter():
if child.tag == '{http://www.opengis.net/wmts/1.0}Capabilities':
print('Version: ' + child.get('version'))
if child.tag == '{http://www.opengis.net/ows/1.1}ServiceType':
print('Service: ' + child.text)
if child.tag == '{http://www.opengis.net/ows/1.1}OperationsMetadata':
print('Request: ')
for e in child:
print('\t ' + e.get('name'))
# Parse layer attributes and vector information.
for child in WmtsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wmts/1.0}Layer"):
if '{http://www.opengis.net/wmts/1.0}Layer' == layer.tag:
f = layer.find("{http://www.opengis.net/ows/1.1}Identifier")
if f is not None:
if f.text == 'MODIS_Aqua_Thermal_Anomalies_All':
# Layer name.
print('Layer: ' + f.text)
# All elements and attributes:
# BoundingBox.
e = layer.find("{http://www.opengis.net/ows/1.1}WGS84BoundingBox")
if e is not None:
print("\t crs: " + e.get('crs'))
print("\t UpperCorner: " + e.find("{http://www.opengis.net/ows/1.1}UpperCorner").text)
print("\t LowerCorner: " + e.find("{http://www.opengis.net/ows/1.1}LowerCorner").text)
# TileMatrixSet.
e = layer.find("{http://www.opengis.net/wmts/1.0}TileMatrixSetLink")
if e is not None:
print("\t TileMatrixSet: " + e.find("{http://www.opengis.net/wmts/1.0}TileMatrixSet").text)
# Time extent.
e = layer.find("{http://www.opengis.net/wmts/1.0}Dimension")
if e is not None:
all = e.findall("{http://www.opengis.net/wmts/1.0}Value")
if all is not None:
print("\t TimeExtent: ")
for g in all:
print("\t\t " + g.text)
# Format.
e = layer.find("{http://www.opengis.net/wmts/1.0}Format")
if e is not None:
print("\t Format: " + e.text)
# Style.
e = layer.find("{http://www.opengis.net/wmts/1.0}Style")
if e is not None:
g=e.find("{http://www.opengis.net/ows/1.1}Identifier")
if g is not None:
print("\t Style: " + g.text)
# Template.
e = layer.find("{http://www.opengis.net/wmts/1.0}ResourceURL")
if e is not None:
print("\t Template: " + e.get('template'))
# Vector metadata.
for e in layer.findall("{http://www.opengis.net/ows/1.1}Metadata"):
if "vector-metadata" in e.get("{http://www.w3.org/1999/xlink}href"):
vectorMetadata=e.get("{http://www.w3.org/1999/xlink}href")
print('\t Vector metadata: ' + vectorMetadata)
response = urllib.request.urlopen(vectorMetadata)
# Load to json.
data = json.loads(response.read())
# Parse json.
for p in data['mvt_properties']:
keys = list(p.keys())
if 'Identifier' in keys:
print('\t\t Identifier: ' + p['Identifier'])
if 'Title' in keys:
print('\t\t Title: ' + p['Title'])
if 'Description' in keys:
print('\t\t Description: ' + p['Description'])
if 'Units' in keys:
print('\t\t Units: ' + p['Units'])
if 'DataType' in keys:
print('\t\t DataType: ' + p['DataType'])
if 'ValueRanges' in keys:
print('\t\t ValueRanges: ' + str(p['ValueRanges']))
if 'ValueMap' in keys:
print('\t\t ValueMap: ' + str(p['ValueMap']))
if 'Function' in keys:
print('\t\t Function: ' + p['Function'])
if 'IsOptional' in keys:
print('\t\t IsOptional: ' + str(p['IsOptional']))
if 'IsLabel' in keys:
print('\t\t IsLabel: ' + str(p['IsLabel']))
print('\n')
# There two vector metadata. Only need one, so break.
break
print('')
# Get general information of WMTS from XML tree.
for child in WmtsTree.iter():
if child.tag == '{http://www.opengis.net/wmts/1.0}Capabilities':
print('Version: ' + child.get('version'))
if child.tag == '{http://www.opengis.net/ows/1.1}ServiceType':
print('Service: ' + child.text)
if child.tag == '{http://www.opengis.net/ows/1.1}OperationsMetadata':
print('Request: ')
for e in child:
print('\t ' + e.get('name'))
# Parse layer attributes and vector information.
for child in WmtsTree.iter():
for layer in child.findall("./{http://www.opengis.net/wmts/1.0}Layer"):
if '{http://www.opengis.net/wmts/1.0}Layer' == layer.tag:
f = layer.find("{http://www.opengis.net/ows/1.1}Identifier")
if f is not None:
if f.text == 'MODIS_Aqua_Thermal_Anomalies_All':
# Layer name.
print('Layer: ' + f.text)
# All elements and attributes:
# BoundingBox.
e = layer.find("{http://www.opengis.net/ows/1.1}WGS84BoundingBox")
if e is not None:
print("\t crs: " + e.get('crs'))
print("\t UpperCorner: " + e.find("{http://www.opengis.net/ows/1.1}UpperCorner").text)
print("\t LowerCorner: " + e.find("{http://www.opengis.net/ows/1.1}LowerCorner").text)
# TileMatrixSet.
e = layer.find("{http://www.opengis.net/wmts/1.0}TileMatrixSetLink")
if e is not None:
print("\t TileMatrixSet: " + e.find("{http://www.opengis.net/wmts/1.0}TileMatrixSet").text)
# Time extent.
e = layer.find("{http://www.opengis.net/wmts/1.0}Dimension")
if e is not None:
all = e.findall("{http://www.opengis.net/wmts/1.0}Value")
if all is not None:
print("\t TimeExtent: ")
for g in all:
print("\t\t " + g.text)
# Format.
e = layer.find("{http://www.opengis.net/wmts/1.0}Format")
if e is not None:
print("\t Format: " + e.text)
# Style.
e = layer.find("{http://www.opengis.net/wmts/1.0}Style")
if e is not None:
g=e.find("{http://www.opengis.net/ows/1.1}Identifier")
if g is not None:
print("\t Style: " + g.text)
# Template.
e = layer.find("{http://www.opengis.net/wmts/1.0}ResourceURL")
if e is not None:
print("\t Template: " + e.get('template'))
# Vector metadata.
for e in layer.findall("{http://www.opengis.net/ows/1.1}Metadata"):
if "vector-metadata" in e.get("{http://www.w3.org/1999/xlink}href"):
vectorMetadata=e.get("{http://www.w3.org/1999/xlink}href")
print('\t Vector metadata: ' + vectorMetadata)
response = urllib.request.urlopen(vectorMetadata)
# Load to json.
data = json.loads(response.read())
# Parse json.
for p in data['mvt_properties']:
keys = list(p.keys())
if 'Identifier' in keys:
print('\t\t Identifier: ' + p['Identifier'])
if 'Title' in keys:
print('\t\t Title: ' + p['Title'])
if 'Description' in keys:
print('\t\t Description: ' + p['Description'])
if 'Units' in keys:
print('\t\t Units: ' + p['Units'])
if 'DataType' in keys:
print('\t\t DataType: ' + p['DataType'])
if 'ValueRanges' in keys:
print('\t\t ValueRanges: ' + str(p['ValueRanges']))
if 'ValueMap' in keys:
print('\t\t ValueMap: ' + str(p['ValueMap']))
if 'Function' in keys:
print('\t\t Function: ' + p['Function'])
if 'IsOptional' in keys:
print('\t\t IsOptional: ' + str(p['IsOptional']))
if 'IsLabel' in keys:
print('\t\t IsLabel: ' + str(p['IsLabel']))
print('\n')
# There two vector metadata. Only need one, so break.
break
print('')
Version: 1.0.0
Service: OGC WMTS
Request:
GetCapabilities
GetTile
Layer: MODIS_Aqua_Thermal_Anomalies_All
crs: urn:ogc:def:crs:OGC:2:84
UpperCorner: 180 90
LowerCorner: -180 -90
TileMatrixSet: 1km
TimeExtent:
2002-07-04/2002-07-29/P1D
2002-08-08/2002-09-12/P1D
2002-09-14/2020-08-16/P1D
2020-09-02/2023-01-31/P1D
Format: application/vnd.mapbox-vector-tile
Style: default
Template: https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/MODIS_Aqua_Thermal_Anomalies_All/default/{Time}/{TileMatrixSet}/{TileMatrix}/{TileRow}/{TileCol}.mvt
Vector metadata: https://gibs.earthdata.nasa.gov/vector-metadata/v1.0/FIRMS_MODIS_Thermal_Anomalies.json
Identifier: LATITUDE
Title: Latitude
Description: Latitude in Decimal Degrees
Units: °
DataType: float
Function: Describe
IsOptional: False
IsLabel: False
Identifier: LONGITUDE
Title: Longitude
Description: Longitude in Decimal Degrees
Units: °
DataType: float
Function: Describe
IsOptional: False
IsLabel: False
Identifier: BRIGHTNESS
Title: Brightness Temperature (Channel 21/22)
Description: Channel 21/22 brightness temperature of the fire pixel, measured in Kelvin.
Units: Kelvin
DataType: float
ValueRanges: [{'Min': 0, 'Max': 500}]
Function: Style
IsOptional: False
IsLabel: False
Identifier: BRIGHT_T31
Title: Brightness Temperature (Channel 31)
Description: Channel 31 brightness temperature of the fire pixel, measured in Kelvin.
Units: Kelvin
DataType: float
ValueRanges: [{'Min': 0, 'Max': 500}]
Function: Style
IsOptional: False
IsLabel: False
Identifier: FRP
Title: Fire Radiative Power
Description: The Fire Radiative Power (FRP) is a measure of the rate of radiant heat output from a fire. It has been demonstrated in small-scale experimental fires that the FRP of a fire is related to the rate at which fuel is being consumed (Wooster et al., 2005) and smoke emissions released (Freeborn et al., 2008).
Units: MW
DataType: float
ValueRanges: [{'Min': 0, 'Max': 20000}]
Function: Style
IsOptional: False
IsLabel: False
Identifier: CONFIDENCE
Title: Detection Confidence
Description: This value is based on a collection of intermediate algorithm quantities used in the detection process. It is intended to help users gauge the quality of individual hotspot/fire pixels. Confidence estimates range between 0 and 100%.
Units: %
DataType: int
ValueRanges: [{'Min': 0, 'Max': 100}]
Function: Style
IsOptional: False
IsLabel: False
Identifier: DAYNIGHT
Title: Day/Night Flag
Description: Indicates whether the fire point was observed during the ‘day’ or ‘night’.
DataType: string
ValueMap: {'D': 'Daytime Fire', 'N': 'Nighttime Fire'}
Function: Describe
IsOptional: False
IsLabel: False
Identifier: SCAN
Title: Along-Scan Pixel Size
Description: The algorithm produces 1km pixels at nadir, but pixels get bigger toward the edge of scan. This value reflects the actual along-scan pixel size.
Units: km
DataType: float
ValueRanges: [{'Min': 1.0, 'Max': 5.0}]
Function: Style
IsOptional: False
IsLabel: False
Identifier: TRACK
Title: Along-Track Pixel Size
Description: The algorithm produces 1km pixels at nadir, but pixels get bigger toward the edge of scan. This value reflects the actual along-track pixel size.
Units: km
DataType: float
ValueRanges: [{'Min': 1.0, 'Max': 2.0}]
Function: Style
IsOptional: False
IsLabel: False
Identifier: ACQ_DATE
Title: Acquisition Date
Description: The date of acquisition for this fire point. (YYYY-MM-DD)
DataType: string
Function: Describe
IsOptional: False
IsLabel: False
Identifier: ACQ_TIME
Title: Acquisition Time
Description: The time of acquisition for this fire point, in UTC. (HHMM)
DataType: string
Function: Describe
IsOptional: False
IsLabel: True
Identifier: SATELLITE
Title: Satellite
Description: Satellite from which the fire is observed.
DataType: string
ValueMap: {'A': 'Aqua', 'T': 'Terra', 'Aqua': 'Aqua', 'Terra': 'Terra'}
Function: Describe
IsOptional: False
IsLabel: False
Identifier: VERSION
Title: Collection and Source
Description: The collection (e.g. MODIS Collection 6) and source of data processing: Near Real-Time (NRT suffix added to collection) or Standard Processing (collection only).
DataType: string
ValueMap: {'6.1NRT': 'Collection 6.1 Near Real-Time processing', '6.1': 'Collection 6.1 Standard processing', '6.0NRT': 'Collection 6 Near Real-Time processing', '6.0': 'Collection 6 Standard processing', '6.02': 'Collection 6 Standard processing', '6.03': 'Collection 6 Standard processing'}
Function: Describe
IsOptional: False
IsLabel: False
Identifier: UID
Title: Unique Identifier
Description: Unique identifier of the data point.
DataType: int
Function: Identify
IsOptional: False
IsLabel: False
Read WMTS Vector Data¶
This example shows how to get WMTS vector data from a Mapbox Vector Tile (MVT). Also shows how to parse vector data values.
In [18]:
Copied!
# Vector data format.
'''
{
'MODIS_Aqua_Thermal_Anomalies_All':
{
'extent': 4096,
'version': 1,
'features':
[
{'geometry':
{'type': 'Point',
'coordinates': [4028, 3959]},
'properties': {'LATITUDE': 35.397,
'LONGITUDE': -90.3,
'BRIGHTNESS': 307.3,
'SCAN': 3.2,
'TRACK': 1.7,
'ACQ_DATE': '2020-10-01',
'ACQ_TIME': '18:30',
'SATELLITE': 'A',
'CONFIDENCE': 48,
'VERSION': '6.0NRT',
'BRIGHT_T31': 296.0,
'FRP': 21.4,
'DAYNIGHT': 'D',
'UID': 13159},
'id': 0,
'type': 1
}
}
}
,,,
]
}
'''
# Below both kvp and restful methods work.
'''
kvp = 'https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/wmts.cgi?\
TIME=2020-10-01T00:00:00Z&FORMAT=application/vnd.mapbox-vector-tile&\
layer=MODIS_Aqua_Thermal_Anomalies_All&tilematrixset=1km&\
Service=WMTS&Request=GetTile&Version=1.0.0&TileMatrix=4&TileCol=3&TileRow=3'
response = requests.get(kvp)
'''
restful = 'https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/MODIS_Aqua_Thermal_Anomalies_All\
/default/2020-10-01T00:00:00Z/1km/4/3/4.mvt'
# Request data.
response = requests.get(restful)
# Parse vector values.
data = response.content
dataDictionary = mapbox_vector_tile.decode(data)
for key in dataDictionary.keys():
parameterDictionary = dataDictionary[key]
features = parameterDictionary['features']
# Print vector data format.
#print(features)
lat = []
lon = []
brightness = []
for f in features:
p = f['properties']
lat.append(p['LATITUDE'])
lon.append(p['LONGITUDE'])
brightness.append(p['BRIGHTNESS'])
print('lat number: ' + str(len(lat)))
print(str(lat))
print('lon number: ' + str(len(lon)))
print(str(lon))
print('brightness number: ' + str(len(brightness)))
print('brightness min: ' + str(min(brightness)))
print('brightness min: ' + str(max(brightness)))
print(str(brightness))
print('')
# Vector data format.
'''
{
'MODIS_Aqua_Thermal_Anomalies_All':
{
'extent': 4096,
'version': 1,
'features':
[
{'geometry':
{'type': 'Point',
'coordinates': [4028, 3959]},
'properties': {'LATITUDE': 35.397,
'LONGITUDE': -90.3,
'BRIGHTNESS': 307.3,
'SCAN': 3.2,
'TRACK': 1.7,
'ACQ_DATE': '2020-10-01',
'ACQ_TIME': '18:30',
'SATELLITE': 'A',
'CONFIDENCE': 48,
'VERSION': '6.0NRT',
'BRIGHT_T31': 296.0,
'FRP': 21.4,
'DAYNIGHT': 'D',
'UID': 13159},
'id': 0,
'type': 1
}
}
}
,,,
]
}
'''
# Below both kvp and restful methods work.
'''
kvp = 'https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/wmts.cgi?\
TIME=2020-10-01T00:00:00Z&FORMAT=application/vnd.mapbox-vector-tile&\
layer=MODIS_Aqua_Thermal_Anomalies_All&tilematrixset=1km&\
Service=WMTS&Request=GetTile&Version=1.0.0&TileMatrix=4&TileCol=3&TileRow=3'
response = requests.get(kvp)
'''
restful = 'https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/MODIS_Aqua_Thermal_Anomalies_All\
/default/2020-10-01T00:00:00Z/1km/4/3/4.mvt'
# Request data.
response = requests.get(restful)
# Parse vector values.
data = response.content
dataDictionary = mapbox_vector_tile.decode(data)
for key in dataDictionary.keys():
parameterDictionary = dataDictionary[key]
features = parameterDictionary['features']
# Print vector data format.
#print(features)
lat = []
lon = []
brightness = []
for f in features:
p = f['properties']
lat.append(p['LATITUDE'])
lon.append(p['LONGITUDE'])
brightness.append(p['BRIGHTNESS'])
print('lat number: ' + str(len(lat)))
print(str(lat))
print('lon number: ' + str(len(lon)))
print(str(lon))
print('brightness number: ' + str(len(brightness)))
print('brightness min: ' + str(min(brightness)))
print('brightness min: ' + str(max(brightness)))
print(str(brightness))
print('')
lat number: 81 [35.397, 35.403, 35.405, 35.446, 35.45, 35.454, 35.467, 35.47, 35.484, 36.053, 33.225, 33.449, 33.45, 33.451, 33.451, 33.625, 33.627, 33.755, 33.756, 33.77, 33.8, 33.803, 33.853, 33.853, 33.855, 33.867, 34.076, 34.346, 34.356, 34.457, 31.032, 31.034, 31.046, 31.048, 31.365, 31.608, 31.89, 31.892, 31.899, 18.316, 19.026, 19.091, 19.094, 19.232, 19.592, 19.594, 19.639, 19.653, 19.8, 19.808, 19.961, 20.518, 20.601, 20.611, 20.723, 20.725, 21.151, 21.726, 21.728, 21.911, 22.04, 22.917, 23.775, 25.189, 25.191, 25.228, 25.723, 25.843, 27.436, 27.447, 28.799, 28.801, 28.914, 29.095, 29.096, 29.527, 29.529, 29.529, 29.531, 29.561, 29.805] lon number: 81 [-90.3, -90.272, -90.266, -90.676, -90.67, -90.641, -92.205, -92.215, -92.209, -89.904, -91.815, -94.62, -94.167, -94.593, -94.589, -93.99, -93.982, -94.507, -94.516, -94.509, -93.806, -93.798, -94.621, -94.626, -94.594, -94.627, -96.957, -91.197, -91.157, -91.021, -95.208, -95.182, -95.209, -95.183, -98.348, -95.122, -90.842, -90.851, -93.643, -100.281, -102.094, -102.919, -101.519, -101.459, -102.482, -102.472, -102.49, -101.918, -102.7, -102.712, -101.019, -100.926, -101.311, -101.245, -103.487, -103.493, -102.555, -102.297, -102.287, -104.148, -99.446, -99.056, -105.361, -99.646, -99.635, -98.347, -103.445, -98.224, -97.681, -97.682, -100.546, -100.535, -98.055, -98.17, -98.155, -97.284, -97.277, -97.267, -97.26, -98.374, -104.591] brightness number: 81 brightness min: 304.2 brightness min: 346.3 [307.3, 312.5, 309.7, 325.2, 335.2, 306.3, 331.4, 325.5, 305.6, 307.8, 304.9, 310.6, 305.1, 317.2, 312.8, 307.3, 310.3, 304.5, 309.3, 311.7, 304.2, 304.8, 343.0, 339.2, 307.3, 325.3, 308.9, 320.4, 341.2, 306.5, 319.4, 332.3, 328.0, 346.3, 317.0, 320.6, 305.4, 305.1, 307.0, 315.7, 319.5, 315.8, 311.3, 313.7, 315.6, 316.2, 312.0, 308.6, 313.6, 313.6, 326.1, 323.6, 318.2, 322.2, 320.2, 318.6, 316.3, 320.1, 318.7, 313.2, 316.2, 318.1, 311.3, 324.7, 326.4, 327.3, 325.8, 327.6, 327.4, 330.9, 326.0, 335.6, 317.2, 319.0, 324.2, 323.9, 333.0, 327.9, 316.7, 322.5, 326.2]
Display WMTS Vector Data¶
This example shows how to overlay WMTS vector data values from last cell on a map with a legend.
In [19]:
Copied!
# Setup map size and projection.
fig = plt.figure(figsize = (8, 5))
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 0))
# x min, x max, y min, y max.
extent = (-130,-30,-5,40)
ax.set_extent(extent)
# Plot lat, lon and brightness.
cmp = ax.scatter(lon, lat, c = brightness, cmap = 'hot')
# Plot legend.
cb = plt.colorbar(cmp, orientation='vertical',
fraction = 0.1, pad = 0.05, shrink = 0.8, label = 'Brightness'
).outline.set_visible(True)
# Draw background.
ax.stock_img()
ax.coastlines()
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth=1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 60, 120, -120, -60, 0])
gl.ylocator = mticker.FixedLocator([-60, -30, 0, 30, 60])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
plt.title('WMTS Vector Data Brightness',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
# Setup map size and projection.
fig = plt.figure(figsize = (8, 5))
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 0))
# x min, x max, y min, y max.
extent = (-130,-30,-5,40)
ax.set_extent(extent)
# Plot lat, lon and brightness.
cmp = ax.scatter(lon, lat, c = brightness, cmap = 'hot')
# Plot legend.
cb = plt.colorbar(cmp, orientation='vertical',
fraction = 0.1, pad = 0.05, shrink = 0.8, label = 'Brightness'
).outline.set_visible(True)
# Draw background.
ax.stock_img()
ax.coastlines()
# Draw grid.
gl = ax.gridlines(ccrs.PlateCarree(), linewidth=1, color = 'blue', alpha = 0.3, draw_labels = True)
gl.top_labels = False
gl.right_labels = False
gl.xlines = True
gl.ylines = True
gl.xlocator = mticker.FixedLocator([0, 60, 120, -120, -60, 0])
gl.ylocator = mticker.FixedLocator([-60, -30, 0, 30, 60])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'color': 'blue'}
gl.ylabel_style = {'color': 'blue'}
plt.title('WMTS Vector Data Brightness',\
fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.show()
print('')
Visualize WMTS Raster Data By Cartopy¶
This example shows how to display WMTS raster data on a map using Cartopy.
In [8]:
Copied!
# Define the WMTS URL
wmts_url = "https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/wmts.cgi"
# Create a map with PlateCarree projection
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([-180, 180, -90, 90]) # World extent
# Add WMTS layer with a specific date
layer = "VIIRS_SNPP_SurfaceReflectance_BandsM11-M7-M5"
time = "2023-01-01"
ax.add_wmts(wmts_url, layer_name=layer, wmts_kwargs={"time": time})
plt.title('Land Surface Reflectance (Bands M11-M7-M5, Best Available, VIIRS, Suomi NPP)')
plt.show()
# Define the WMTS URL
wmts_url = "https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/wmts.cgi"
# Create a map with PlateCarree projection
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([-180, 180, -90, 90]) # World extent
# Add WMTS layer with a specific date
layer = "VIIRS_SNPP_SurfaceReflectance_BandsM11-M7-M5"
time = "2023-01-01"
ax.add_wmts(wmts_url, layer_name=layer, wmts_kwargs={"time": time})
plt.title('Land Surface Reflectance (Bands M11-M7-M5, Best Available, VIIRS, Suomi NPP)')
plt.show()