Skip to content

Instantly share code, notes, and snippets.

@jduckles
Last active May 2, 2022
Embed
What would you like to do?

Select Worldview Snapshot

Put this together based on the following twitter prompt

Draw a box you want at World View Snapshosts I tend to set the spatial resolution to 250m, but that will limit you to smaller total areas.

Click Preview

In the preview click the checkbox "Always show today"

Place that link in download_modis.sh as SNAPSHOT_URL below and setup a cron job, noon local time or a bit after is usually a safe bet, set appropriate target paths to keep things organized if you like.

Once you have several images downloaded, say 1-week, 1-month, 1-year. You can try to do some cloud free compositing. What has worked best for me in NZ is to do a lower quantile value, so the lowest 10% across the pixels for a week/month. NZ is often covered in clouds, so this helps ensure that we don't pick a value from the series that has clouds in it. There are a few otehr methods that you can try coded in, but I really don't use them, they were just for experimental purposes.

Example of running the python script.

python3 stack_images.py nz_2021-09 quantile

drops an image called nz-2021-09_quantile_merged.jpg The nz_2021-09 bit does a pattern match on the files in the current directory, so make sure the script can find the files. See L#12 in the stack_images.py

#!/bin/sh
SNAPSHOT_URL="https://wvs.earthdata.nasa.gov/api/v1/snapshot?REQUEST=GetSnapshot&LAYERS=MODIS_Terra_CorrectedReflectance_TrueColor&CRS=EPSG:4326&TIME=P0D&WRAP=DAY&BBOX=-47.88515,166.18059,-35.578732,178.970474&FORMAT=image/jpeg&WIDTH=5821&HEIGHT=5601&AUTOSCALE=TRUE&ts=1600572288374"
wget "${SNAPSHOT_URL}" -O nz_$(date +%Y-%m-%d).jpg
ln -sf nz_$(date +%Y-%m-%d).jpg nz_latest.jpg # Optional but handy, this feeds a home-automation dashboard to show me a fresh pic everyday.
#!/usr/bin/python3
import os
import sys
from PIL import Image
import numpy as np
print("\nstart stacking images...")
prefix = sys.argv[1]
method = sys.argv[2]
print(prefix)
images = [ file for file in os.listdir() if file[0:len(prefix)] == prefix ]
red = []
green = []
blue = []
for image in images:
print('.')
with Image.open(image) as img:
if image == images[0]:
img_size = img.size
img_array = np.asarray(img)
print('Initial shape {}'.format(img_array.shape))
red.append(img_array[:,:,0])
green.append(img_array[:,:,1])
blue.append(img_array[:,:,2])
rgbArray = np.zeros(img_array.shape, 'uint8')
if method == 'quantile':
rgbArray[...,0] = np.quantile(np.stack(red, -1), q=[0.10], axis=2)
rgbArray[...,1] = np.quantile(np.stack(green, -1), q=[0.10], axis=2)
rgbArray[...,2] = np.quantile(np.stack(blue, -1), q=[0.10], axis=2)
elif method == 'median':
rgbArray[...,0] = np.median(np.stack(red, -1), axis=2)
rgbArray[...,1] = np.median(np.stack(green, -1), axis=2)
rgbArray[...,2] = np.median(np.stack(blue, -1), axis=2)
elif method == 'min':
rgbArray[...,0] = np.min(np.stack(red, -1), axis=2)
rgbArray[...,1] = np.min(np.stack(green, -1), axis=2)
rgbArray[...,2] = np.min(np.stack(blue, -1), axis=2)
elif method == 'max':
rgbArray[...,0] = np.max(np.stack(red, -1), axis=2)
rgbArray[...,1] = np.max(np.stack(green, -1), axis=2)
rgbArray[...,2] = np.max(np.stack(blue, -1), axis=2)
elif method == 'stdev':
std_red = np.std(np.stack(red, -1), axis=2)
std_green = np.std(np.stack(green, -1), axis=2)
std_blue = np.std(np.stack(blue, -1), axis=2)
mean_red = np.mean(np.stack(red, -1), axis=2)
mean_green = np.mean(np.stack(green, -1), axis=2)
mean_blue = np.mean(np.stack(blue, -1), axis=2)
rgbArray[...,0] = mean_red - std_red
rgbArray[...,1] = mean_green - std_green
rgbArray[...,2] = mean_blue - std_blue
pillImage = Image.fromarray(rgbArray,"RGB")
pillImage.save('{prefix}_{method}_merged.jpg'.format(prefix=prefix, method=method), "JPEG")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment