Accessing via Python¶

# 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')

# 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)

# 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('')

#  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 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 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 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 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 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('')

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

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)

# 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('')

# 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)

# 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('')

# 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('')

# 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('')

# Make GDAL input XML file like globe.xml below.
xml = xmltree.parse("python-examples/globe.xml")
pretty = xmltree.tostring(xml, encoding="unicode", pretty_print=True)
print(pretty)

# Run GDAL command.
cmd  = 'gdal_translate -of JPEG -outsize 1200 600 -projwin -180 90 180 -90 python-examples/globe.xml python-examples/globe.jpg'
os.system(cmd)

# Output image is globe.jpg.
img = plimg.open('python-examples/globe.jpg')

# Setup map size and projection.
fig = plt.figure(figsize = (5, 10), dpi = 100)
ax = plt.axes([1, 1, 1, 1], projection = ccrs.PlateCarree(central_longitude = 0))
ax.set_xlim([-180, 180])
ax.set_ylim([-90, 90])

# Display image on map.
imgExtent = (-180,180,-90,90)
cmp = plt.imshow(img, extent = imgExtent)

# Draw coastline.
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])
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'}

plt.title('WMTS Terra True Color Image By GDAL',\
          fontname = "Times New Roman", fontsize = 20, color = 'green')

plt.show()

print('')

layerName = "Agricultural_Lands_Croplands_2000"
legendURLHorizontal = None
legendURLVertical = None

# 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 == layerName:                    
                    # Style tag
                    e = layer.find("{http://www.opengis.net/wmts/1.0}Style")
                    if e is not None:
                        for legendTag in e.findall("{http://www.opengis.net/wmts/1.0}LegendURL"):
                            # Horizontal legend
                            if legendTag.attrib["{http://www.w3.org/1999/xlink}role"] == "http://earthdata.nasa.gov/gibs/legend-type/horizontal":
                                legendURLHorizontal = legendTag.attrib["{http://www.w3.org/1999/xlink}href"]
                            # Vertical legend
                            else:
                                legendURLVertical = legendTag.attrib["{http://www.w3.org/1999/xlink}href"]

print("Horizontal Legend:\n")
display(Image(url=legendURLHorizontal))
print("Vertical Legend:\n")
display(Image(url=legendURLVertical))

# Make GDAL input XML file like globe_cropland.xml below.
xml = xmltree.parse("python-examples/globe_cropland.xml")
pretty = xmltree.tostring(xml, encoding="unicode", pretty_print=True)
print(pretty)

# cairosvg will be needed to convert the legend from SVG to PNG
import cairosvg

# Run GDAL command.
cmd  = 'gdal_translate -of PNG -outsize 1200 600 -projwin -180 90 180 -90 python-examples/globe_cropland.xml python-examples/globe_cropland.png'
os.system(cmd)

# Output image is globe_cropland.png.
img = plimg.open('python-examples/globe_cropland.png')

# Setup map size and projection.
fig = plt.figure(figsize = (5, 10), dpi = 100)
ax = plt.axes([1, 1, 1, 1], projection = ccrs.PlateCarree(central_longitude = 0))
ax.set_xlim([-180, 180])
ax.set_ylim([-90, 90])

# Display image on map.
imgExtent = (-180,180,-90,90)
cmp = plt.imshow(img, extent = imgExtent)

# 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.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'}

plt.title('WMTS Croplands (Global Agricultural Lands, 2000) By GDAL',\
        fontname = "Times New Roman", fontsize = 20, color = 'green')

# Horizontal legend:
# Get the legend image from the URL and convert it to PNG
legend_png_h = cairosvg.svg2png(urllib.request.urlopen(legendURLHorizontal).read())
# Convert to numpy array for matplotlib
legendImgArr_h = np.array(plimg.open(BytesIO(legend_png_h)))

# use data coordinates to specify the position and dimensions of new inset axes
axin_h = ax.inset_axes([-125,-260,250,250],transform=ax.transData)
axin_h.imshow(legendImgArr_h)
axin_h.axis('off')

# Vertical legend:
# Get the legend image from the URL and convert it to PNG
legend_png_v = cairosvg.svg2png(urllib.request.urlopen(legendURLVertical).read())
# Convert to numpy array for matplotlib
legendImgArr_v = np.array(plimg.open(BytesIO(legend_png_v)))

# use data coordinates to specify the position and dimensions of new inset axes
axin_v = ax.inset_axes([135,-115,200,200],transform=ax.transData)
axin_v.imshow(legendImgArr_v)
axin_v.axis('off')

plt.show()

print('')

from collections import Counter, OrderedDict
import re
import math

layerName = "Agricultural_Lands_Croplands_2000"
colormapURL = None

# Get the colormap URL by searching the GetCapabilities 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:
                if f.text == layerName:     
                    # Metadata tags
                    e = layer.findall("{http://www.opengis.net/ows/1.1}Metadata")
                    for metadataTag in e:
                        if metadataTag.attrib["{http://www.w3.org/1999/xlink}role"] == "http://earthdata.nasa.gov/gibs/metadata-type/colormap/1.3":
                            colormapURL = metadataTag.attrib["{http://www.w3.org/1999/xlink}href"]
print("ColorMap URL:", colormapURL)

# Obtain and read the colormap XML
response = requests.get(colormapURL)
colormapXML = xmltree.fromstring(response.content)

# Create a dictionary to be used to map RGB values to their bin names
rgb_bin_map = {}

colormaps = colormapXML.getroottree().findall("ColorMap")
for cmap in colormaps:
    if cmap.attrib["title"] == "Croplands Land Coverage":
        entries = cmap.find("Entries").findall("ColorMapEntry")
        for entry in entries:
            # Parse RGB values from string to 3-tuple of ints
            rgb_tuple = tuple(map(int, entry.attrib['rgb'].split(',')))
            # Parse the value's range string in the format [<lower_decimal>,<upper_decimal>)
            value_match = re.match(r'\[(\d+\.\d+),(\d+\.\d+)\)', entry.attrib['value'])
            if value_match:
                # Our bins will have ranges of 0.1 because our dataset goes from 0.0 to 1.0.
                # Construct the value's bin name:
                # (1) Floor the decimal component of the lower bound
                lowerbound_decimal = float(value_match.group(1).split('.')[1])
                lowerbound_floored = math.floor(lowerbound_decimal * 0.1)
                # (2) Create a string for the lower bound
                lowerbin_str = str(lowerbound_floored * 10)
                # (3) Obtain the upper bound as a string
                upperbin_str = str(int(lowerbin_str) + 9)

                # Assign the bin name to the RGB value in the dictionary
                bin_name = "{0}% - {1}%".format(lowerbin_str, upperbin_str)
                rgb_bin_map[rgb_tuple] = bin_name
                

# Create a list keeping track of each pixel's bin
percents = []

# Search the image for occurrences of each color
# First, open the image that GDAL created in the "Visualize WMTS Raster Data with a Legend" example
map_img = plimg.open("python-examples/globe_cropland.png")

# Iterate over the image and put pixels in bins
for x in range(map_img.width):
    for y in range(map_img.height):
        # Use rgb_bins to map the RGB value to its corresponding bin name
        try:
            bin_name = rgb_bin_map[map_img.getpixel((x,y))]
            percents.append(bin_name)
        # "No Data" values are ignored
        except KeyError:
            pass

# Sort the list of bin names to guarentee that the bar chart x-axis will be properly ordered
bin_names = list(rgb_bin_map.values())
bin_names.sort()

# Create an OrderedDict to maintain the sorted order with a default count of 0 for each bin.
# The default count of 0 guarentees each bin a spot on the chart even if it's empty.
rgbs_count = OrderedDict.fromkeys(bin_names, 0)
# Calculate the count of each bin's contents and merge into the OrderedDict.
rgbs_count.update(dict(Counter(percents)))

# Create the plot
fig, axes = plt.subplots()
axes.bar(rgbs_count.keys(), rgbs_count.values())
plt.title('Frequencies of Percent Cropland',\
        fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.xlabel("Percentages")
plt.ylabel("Number of Occurences")
plt.setp(axes.get_xticklabels(), rotation=30, horizontalalignment='right');

# Make GDAL input XML file like globe_land_cover.xml below.
xml = xmltree.parse("python-examples/globe_land_cover.xml")
pretty = xmltree.tostring(xml, encoding="unicode", pretty_print=True)
print(pretty)

# Run GDAL command.
cmd  = 'gdal_translate -of PNG -outsize 1200 600 -projwin -180 90 180 -90 python-examples/globe_land_cover.xml python-examples/globe_land_cover.png'
os.system(cmd)

# Output image is globe_land_cover.png.
img = plimg.open('python-examples/globe_land_cover.png')

# Setup map size and projection.
fig = plt.figure(figsize = (5, 10), dpi = 100)
ax = plt.axes([1, 1, 1, 1], projection = ccrs.PlateCarree(central_longitude = 0))
ax.set_xlim([-180, 180])
ax.set_ylim([-90, 90])

# Display image on map.
imgExtent = (-180,180,-90,90)
cmp = plt.imshow(img, extent = imgExtent)

# 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.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'}

plt.title('WMTS Land Cover By GDAL',\
        fontname = "Times New Roman", fontsize = 20, color = 'green')

# Obtain and display Legend:

# First obtain the legend's URL
layerName = "MODIS_Combined_L3_IGBP_Land_Cover_Type_Annual"
legendURL = None

# 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 == layerName:                    
                    # Style tag
                    e = layer.find("{http://www.opengis.net/wmts/1.0}Style")
                    if e is not None:
                        for legendTag in e.findall("{http://www.opengis.net/wmts/1.0}LegendURL"):
                            legendURL = legendTag.attrib["{http://www.w3.org/1999/xlink}href"]

# Next, get the legend image from the URL and convert it to PNG
legend_png= cairosvg.svg2png(urllib.request.urlopen(legendURL).read())
# Convert to numpy array for matplotlib
legendImgArr = np.array(plimg.open(BytesIO(legend_png)))

# use data coordinates to specify the position and dimensions of new inset axes
axin = ax.inset_axes([200,-115,200,200],transform=ax.transData)
axin.imshow(legendImgArr)
axin.axis('off')

plt.show()

print('')

# First obtain the colormap's URL
layerName = "MODIS_Combined_L3_IGBP_Land_Cover_Type_Annual"
colormapURL = None

# Get the colormap URL by searching the GetCapabilities 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:
                if f.text == layerName:     
                    # Metadata tags
                    e = layer.findall("{http://www.opengis.net/ows/1.1}Metadata")
                    for metadataTag in e:
                        if metadataTag.attrib["{http://www.w3.org/1999/xlink}role"] == "http://earthdata.nasa.gov/gibs/metadata-type/colormap/1.3":
                            colormapURL = metadataTag.attrib["{http://www.w3.org/1999/xlink}href"]
print("ColorMap URL:", colormapURL)

# Obtain and read the colormap XML
response = requests.get(colormapURL)
colormapXML = xmltree.fromstring(response.content)

# Search the image for occurrences of each color
# First, open the image that GDAL created earlier
map_img = plimg.open("python-examples/globe_land_cover.png")

# Create a dictionary to keep track of the number of occurences of each pixel
landCoverRGBCounts = {}

# Iterate over the image and count the number of occurences of each pixel
for x in range(map_img.width):
    for y in range(map_img.height):
        landCoverRGBCounts[map_img.getpixel((x,y))] = landCoverRGBCounts.get(map_img.getpixel((x,y)), 0) + 1

# Map the RGB values to their corresponding land cover type
# Create a dictionary to be used to map RGB values and their counts to their land cover types
landCoverTypeCounts = {}

colormaps = colormapXML.getroottree().findall("ColorMap")
for cmap in colormaps:
    if cmap.attrib["title"] == "Classifications":
        entries = cmap.find("Legend").findall("LegendEntry")
        for entry in entries:
            # Parse RGB values from string to 3-tuple of ints
            rgb_tuple = tuple(map(int, entry.attrib['rgb'].split(',')))
            landCoverTypeCounts[entry.attrib['tooltip']] = landCoverRGBCounts[rgb_tuple]

# Sort the result using an ordered dictionary
landCoverSorted = OrderedDict(sorted(landCoverTypeCounts.items(), key=lambda x: x[1], reverse=True))

# Create the plot
fig, axes = plt.subplots()
axes.bar(landCoverSorted.keys(), landCoverSorted.values())
plt.title('Frequencies of Land Cover Types',\
        fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.xlabel("Land Cover Type")
plt.ylabel("Number of Pixel Occurences")
fig.set_figwidth(8)
plt.setp(axes.get_xticklabels(), rotation=30, horizontalalignment='right');

# Remove "Water Bodies"
landCoverWithoutWater = landCoverSorted.copy()
landCoverWithoutWater.pop("Water Bodies")

landCoverWithoutWaterPercentages = [x / sum(landCoverWithoutWater.values()) * 100 for x in landCoverWithoutWater.values()]

fig, axes = plt.subplots()
axes.bar(landCoverWithoutWater.keys(), landCoverWithoutWaterPercentages)
plt.title('Percentages of Land Cover Types',\
        fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.xlabel("Land Cover Type")
plt.ylabel("% Land Coverage")
fig.set_figwidth(8)
plt.setp(axes.get_xticklabels(), rotation=30, horizontalalignment='right');

# Make GDAL input XML file like globe_bands367.xml below.
xml = xmltree.parse("python-examples/globe_bands367.xml")
pretty = xmltree.tostring(xml, encoding="unicode", pretty_print=True)
print(pretty)

# Run GDAL command.
cmd  = 'gdal_translate -of PNG -outsize 1200 600 -projwin -180 90 180 -90 python-examples/globe_bands367.xml python-examples/globe_bands367.png'
os.system(cmd)

# Output image is globe_bands367.png.
img = plimg.open('python-examples/globe_bands367.png')

# Setup map size and projection.
fig = plt.figure(figsize = (5, 10), dpi = 100)
ax = plt.axes([1, 1, 1, 1], projection = ccrs.PlateCarree(central_longitude = 0))
ax.set_xlim([-180, 180])
ax.set_ylim([-90, 90])

# Display image on map.
imgExtent = (-180,180,-90,90)
cmp = plt.imshow(img, extent = imgExtent)

# 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.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'}

plt.title('WMTS Corrected Reflectance (Bands 7-2-1) By GDAL',\
        fontname = "Times New Roman", fontsize = 20, color = 'green')

plt.show();

# We'll be using KMeans clustering to perform the analysis
from sklearn.cluster import KMeans
from collections import Counter

# First, open the image that GDAL created in the example above as a numpy array
map_img = np.asarray(plimg.open("python-examples/globe_bands367.png"))
# map_img is currently a 3D array, with dimensions of row, column, and pixel values.
# Reshape it to form a 2D array of pixel values, with row and column positions disregarded.
map_pixels = map_img.reshape(-1, 3)

# Next, use KMeans clustering to find the most common colors in the image
clustering = KMeans(n_clusters = 4, random_state=1)
clustering.fit(map_pixels)

# Calculate the percentage of pixels that each cluster contains
total_pixels = len(clustering.labels_)
counter = Counter(clustering.labels_)
color_percents = []
for i in counter:
    # Create a list of (color, percentage) tuples
    color_percents.append((tuple(clustering.cluster_centers_[i]), np.round(counter[i]/total_pixels, 3) * 100))

# Sort the colors and the percents together based on the percents
color_percents_sorted = sorted(color_percents, key=lambda x: x[1], reverse=False)
# Obtain the sorted list of colors
colors = [x[0] for x in color_percents_sorted]
# Obtain the sorted list of percents
percents = [x[1] for x in color_percents_sorted]

# Obtain the RGB values on scale of [0,1] for Matplotlib to assign to each bar in the chart
colors_01 = [tuple(map(lambda y: y / 255, x)) for x in colors]

# Obtain properly formatted labels for each bar
labels = list(map(lambda x: "RGB=({:.1f}, {:.1f}, {:.1f})".format(x[0],x[1],x[2]), colors))

# Create the plot
fig, axes = plt.subplots()
bar = axes.barh(labels, percents, color = colors_01)
axes.bar_label(bar)
plt.title('Percentages of Colors for Corrected Reflectance (Bands 3-6-7)',\
        fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.ylabel("Color")
plt.xlabel("Percentage")
fig.set_figwidth(8)
plt.setp(axes.get_xticklabels(), rotation=30, horizontalalignment='right');

img = plimg.open("python-examples/globe_bands367.png")
outimg = img.copy()
# Reshape the list of labels from a 1-D array to match the shape of the image's dimensions
labels_matrix = clustering.labels_.reshape(np.asarray(img).shape[:2])
# Assign each pixel the color specified by its corresponding label
for y in range(outimg.height):
    for x in range(outimg.width):
        outimg.putpixel((x,y), tuple(clustering.cluster_centers_[labels_matrix[y, x]].astype(int)))
display(outimg)

# 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_Bands367 and MODIS_Terra_Data_No_Data
img = wms.getmap(layers=['MODIS_Terra_CorrectedReflectance_Bands367', 'MODIS_Terra_Data_No_Data'],  # Layers
                 srs='epsg:4326',  # Map projection
                 bbox=(-180,-90,180,90),  # Bounds
                 size=(1200, 600),  # Image size
                 time='2021-09-25',  # Time of data
                 format='image/png',  # Image format
                 transparent=True)  # Nodata transparency

Image(img.read())

# Make GDAL input XML files like modis_terra_nodata_mask.xml and osm_land_mask.xml below.
nodata_xml = xmltree.parse("python-examples/modis_terra_nodata_mask.xml")
nodata_pretty = xmltree.tostring(nodata_xml, encoding="unicode", pretty_print=True)
print("modis_terra_nodata_mask.xml:")
print(nodata_pretty)

# Run GDAL commands to download the images.
cmd  = 'gdal_translate -of PNG -outsize 1200 600 -projwin -180 90 180 -90 python-examples/globe_bands367.xml python-examples/globe_bands367.png'
os.system(cmd)

cmd  = 'gdal_translate -of PNG -outsize 1200 600 -projwin -180 90 180 -90 python-examples/modis_terra_nodata_mask.xml python-examples/modis_terra_nodata_mask.png'
os.system(cmd)

# First, open each image that GDAL created above as a numpy array
map_img = np.asarray(plimg.open("python-examples/globe_bands367.png"))
nodata_img = np.asarray(plimg.open("python-examples/modis_terra_nodata_mask.png"))
# map_img and nodata_img are currently 3D arrays, with dimensions of row, column, and pixel values.
# Reshape it to form a 2D array of pixel values, with row and column positions disregarded.
map_pixels = map_img.reshape(-1, 3)
nodata_pixels = nodata_img.reshape(-1, 3)

# Filter map_pixels such that any pixels that do not have values of 0,0,0 in landmask_img or nodata_img are eliminated.
# Create a filter array for each mask
make_filter = lambda p: np.array_equal(p,np.array([0,0,0]))
nodata_filter = np.apply_along_axis(make_filter, 1, nodata_pixels)
# filter the pixels
filtered_pixels = map_pixels[nodata_filter]

# Use KMeans clustering to find the most common colors in the image
clustering = KMeans(n_clusters = 4, random_state=1)
clustering.fit(filtered_pixels)

# Calculate the percentage of pixels that each cluster contains
total_pixels = len(clustering.labels_)
counter = Counter(clustering.labels_)
color_percents = []
for i in counter:
    # Create a list of (color, percentage) tuples
    color_percents.append((tuple(clustering.cluster_centers_[i]), np.round(counter[i]/total_pixels, 3) * 100))

# Sort the colors and the percents together based on the percents
color_percents_sorted = sorted(color_percents, key=lambda x: x[1], reverse=False)
# Obtain the sorted list of colors
colors = [x[0] for x in color_percents_sorted]
# Obtain the sorted list of percents
percents = [x[1] for x in color_percents_sorted]

# Obtain the RGB values on scale of [0,1] for Matplotlib to assign to each bar in the chart
colors_01 = [tuple(map(lambda y: y / 255, x)) for x in colors]

# Obtain properly formatted labels for each bar
labels = list(map(lambda x: "RGB=({:.1f}, {:.1f}, {:.1f})".format(x[0],x[1],x[2]), colors))

# Create the plot
fig, axes = plt.subplots()
bar = axes.barh(labels, percents, color = colors_01)
axes.bar_label(bar)
plt.title('Percentages of Colors in Corrected Reflectance (Bands 3-6-7)',\
        fontname = "Times New Roman", fontsize = 20, color = 'green')
plt.ylabel("Color")
plt.xlabel("Percentage")
fig.set_figwidth(8)
plt.setp(axes.get_xticklabels(), rotation=30, horizontalalignment='right');