More Useful Libraries

[2320sharon]

Useful Libraries

GeeMaps for Data Extraction and Loading Tifs

Website
Youtube Channel

• As I have been digging into ipyleaflet for data extraction I have discovered this library called GeeMaps its built on top of ipyleaflet and has additional functionality like the ability to load in rasters (tif) and more fine control over drawing on the maps and much more.
• Here is a tutorial on how to load rasters into the map.
• His youtube channel is an absolute 🥇 gold mine
• Let me know what you think

[dbuscombe-usgs]

I agree. leafmap and geemap are maintained by the same person, and there is overlap in their functionality

I'm a little familiar with geemap and have explored the image downloading functionality ... that was how I downloaded sentinel-2 data for coast train, before I started using CoastSat. I regret not using CoastSat from the start, however, because it is better at keeping track of the satellite positional error, masks imagery for clouds, and it doesnt appear to have limits on the size of image you can download, unlike geemap

That said, I strongly agree geemap could be really useful and exactly what we need for webmapping, geojson ROI definition, etc

Also, I didn't exhaustively explore image download options using geemap, and things may have improved since the summer

Let's add geemap to the list of things to explore

[2320sharon]

After doing a little bit of research into the key features of both Geemap and Leafmap. Geemaps is optimized for interacting with EE while Leafmap is designed primarily around mapping, rasters, and adding more functionality to the ipyleaflet library.

I think I will try to make some functions work with Leafmap first, then if I need any special information from EE that CoastSat cannot provide I will dig deeper into Geemap's capabilities.

[dbuscombe-usgs]

Here is an example of how I was using geemap

import geemap,ee, reverse_geocoder
import os, json  #shutil
from ipyleaflet import GeoJSON
import numpy as np
from glob import glob
from area import area
from joblib import Parallel, delayed
from shapely.geometry import LineString, MultiPolygon, Polygon
from shapely.ops import split

assuming you have a geojson file containing multiple sites to download imagery, e.g. coastTrain_sites.geojson

with open('coastTrain_sites.geojson') as f:
    json_data = json.load(f)
features = json_data['features']

extract the geographic coordinates (wgs84, or longitudes and latitudes) from the data

L=[]; LL=[]
for feature,site in zip(features, np.arange(1,1+len(features))):
    coordinates = feature['geometry']['coordinates']
    lng, lat = np.mean(coordinates[0], axis=0)
    #print((lng, lat))
    LL.append(lng); L.append(lat)

coordinates = [tuple((a,b)) for a,b in zip(L,LL)]

define start and end date, and maximum cloud cover

start_date = '2010-01-01'
end_date = '2020-12-31'
MAXCLOUD = 5  #percentage maximum cloud cover tolerated

# define S2 collection of GEE
gee_collection = 'COPERNICUS/S2_SR'
OUT_RES_M = 10 #output pixel size in m (10m for sentinel visisble band)

## make a collection
collection = ee.ImageCollection(gee_collection)

### define polygons based on provided coordinates
polys = ee.Geometry.Polygon(coordinates)
centroid = polys.centroid()
lng, lat = centroid.getInfo()['coordinates']
# print("lng = {}, lat = {}".format(lng, lat))

## filter collection based on time and cloud cover
lng_lat = ee.Geometry.Point(lng, lat)
collection = collection.filterBounds(polys)
collection = collection.filterDate(start_date, end_date)
count = collection.size().getInfo()
print("Number of cloudy scenes: ", count)

collection = collection.filter('CLOUDY_PIXEL_PERCENTAGE < '+str(MAXCLOUD))

count = collection.size().getInfo()
print("Number of non-cloudy scenes: ", count)

next,make an image collection

img_lst = collection.toList(1000)
N = []
for i in range(0, count):
    image = ee.Image(img_lst.get(i))
    name = image.get('system:index').getInfo()
    N.append(name)

download data

for n in N:
    print(n)
    image = ee.Image(gee_collection+'/'+n).select(['B4', 'B3', 'B2'])
    try:
        os.remove(os.getcwd()+os.sep+n+'.tif')
    except:
        pass
    geemap.ee_export_image(image, os.getcwd()+os.sep+n+'.tif', scale=OUT_RES_M, region=polys) #, file_per_band=True)

    image = ee.Image(gee_collection+'/'+n).select(['B8'])
    try:
        os.remove(os.getcwd()+os.sep+n+'_ir.tif')
    except:
        pass
    geemap.ee_export_image(image, os.getcwd()+os.sep+n+'_ir.tif', scale=OUT_RES_M, region=polys) #, file_per_band=True)

finally, clean up images into folders

try:
    os.mkdir('sentinel2')
except:
    pass


for f in glob(os.getcwd()+os.sep+'*.tif'):
    print(f)
    shutil.move(f,os.getcwd()+os.sep+'sentinel2')

[2320sharon]

This code is great for getting a simplified understanding of how CoastSat works. Thank you for sharing 🏅

[dbuscombe-usgs]

Just for clarification, this is how you could use geemap to download imagery. It's an alternative to the simpler coastsat way

[dbuscombe-usgs]

I also have another version of this workflow that tries to download chunks of images if the given ROI polygon is too large, i.e. if the image area is larger than the GEE maximum download limit

(I dont yet know how CoastSat circumvents this limitation - I imagine it uses a similar approach to below)

in this version, the polygon is further subdivided into 4 x 4 = 16 equal pieces

from area import area

    ee.Initialize()

    coordinates = feature['geometry']['coordinates']
    lng, lat = np.mean(coordinates[0], axis=0)
    #print((lng, lat))
    area_sqkm = area(feature['geometry'])/1e6
    #print(area_sqkm)

    collection = ee.ImageCollection(gee_collection)

    polygon = Polygon([tuple(c) for c in coordinates[0]])
    minx, miny, maxx, maxy = polygon.bounds

    nx, ny = 4, 4  # number of columns and rows

    dx = (maxx - minx) / nx  # width of a small part
    dy = (maxy - miny) / ny  # height of a small part
    horizontal_splitters = [LineString([(minx, miny + i*dy), (maxx, miny + i*dy)]) for i in range(ny)]
    vertical_splitters = [LineString([(minx + i*dx, miny), (minx + i*dx, maxy)]) for i in range(nx)]
    splitters = horizontal_splitters + vertical_splitters

those pieces are then combined as a Multipolygon object

    result = polygon
    for splitter in splitters:
        result = MultiPolygon(split(result, splitter))
    parts = [list(part.exterior.coords) for part in result.geoms]

then work on each part in turn

    counter = 1
    for part in parts:

        #print(area({'type':'Polygon','coordinates': [[list(p) for p in part]]})/1e6)
        collection = ee.ImageCollection(gee_collection)

        polys = ee.Geometry.Polygon(part)

        centroid = polys.centroid()
        lng, lat = centroid.getInfo()['coordinates']
        #print("lng = {}, lat = {}".format(lng, lat))

        #lng_lat = ee.Geometry.Point(lng, lat)
        collection = collection.filterBounds(polys)
        collection = collection.filterDate(start_date, end_date).sort('system:time_start', True)
        count = collection.size().getInfo()
        #print("Number of cloudy scenes: ", count)
        img_lst = collection.toList(1000)

        N = []
        for i in range(0, count):
            image = ee.Image(img_lst.get(i))
            name = image.get('system:index').getInfo()
            N.append(name)

        for n in N:
            image = ee.Image('USDA/NAIP/DOQQ/'+n)
            #geemap.ee_export_image(image, os.getcwd()+os.sep+'site_'+str(site)+'_'+n+'_part_'+str(counter)+'_of_'+str(len(parts))+'.tif', scale=OUT_RES_M, region=polys, file_per_band=True)
            geemap.ee_export_image(image, os.getcwd()+os.sep+'naip'+os.sep+'site'+str(site)+os.sep+'chunk'+str(counter)+os.sep+'chunk'+str(counter)+'_'+n+'.tif', scale=OUT_RES_M, region=polys, file_per_band=True)
            geemap.ee_export_image(image, os.getcwd()+os.sep+'naip'+os.sep+'site'+str(site)+os.sep+'chunk'+str(counter)+os.sep+'chunk'+str(counter)+'_'+n+'_multiband.tif', scale=OUT_RES_M, region=polys, file_per_band=False)

        counter += 1

as you can see, this is a much more complicated workflow than CoastSat ...

[dbuscombe-usgs]

There are a few potential limitations with these workflows compared to CoastSat (need to be further explored):

1. a hard limit on the size of file you can download through geemap (I forget what it is, perhaps only a few MB)
2. geemap had no cloud masking when I tried it much earlier this year, perhaps it does now?
3. coastsat keeps track of EPSG and precision in the metadata
4. I like the way CoastSat pansharpens the imagery using the PCA method

However, that said, everything is currently on the table. We should explore geemap, wxee, and Coastsat as a mechanisms to download imagery - what are the relative merits?

Geemap could be EXACTLY what we need as a graphical frontend

[2320sharon]

I agree that CoastSat is better optimized for the task of data retrieval due to it having the cloud cover masking, tracking the metadata, amongst many other tricks. I was thinking of using Geemaps in combination with CoastSat as it could as a GUI (like you said) with some great drawing and customization options to make it even easier to download data with CoastSat.

TDLR Geemaps (or whatever our frontend is) could feed the geojson to CoastSat which would then do the data retrieval.

[dbuscombe-usgs]

I would agree with your assessment. Could you start exploring geemap as a mechanism for easy drawing of ROIs, and storing as geojson format?

[2320sharon]

I'll explore Geemaps and Leafmap as options for easy drawing of ROIs, and storing them as a geojson format. I'll make a comparison between them and see which one better fits the project's goal in another discussion.

[dbuscombe-usgs]

Adding this here. So many ways to query and download satellite data. This uses a STAC https://odc-stac.readthedocs.io/en/latest/notebooks/stac-load-e84-aws.html?s=09

[dbuscombe-usgs]

This combines webmapping with xarray. I really like the idea of storing our data cubes in xarray format... https://github.com/davidbrochart/xarray_leaflet

[2320sharon]

I'll dig into ipyleaflet and xarray integration.

[2320sharon]

xarray_leaflet is known to have issues with Windows machines based on the Github issues. I'll test it on my machine, but we may not be able to use it if we want our package to work on multiple kinds of machines.

[dbuscombe-usgs]

On xarray, this is a little old but mostly still valid https://rabernat.github.io/research_computing_2018/xarray.html

[dbuscombe-usgs]

See also https://github.com/jbusecke/xmovie