feat: add sentinel 2 script #335
Conversation
Python script and Topo50 index sheets shapefile for processing Sentinel 2 imagery.
| rgb_files = glob.glob("R10m/*B0[2-4]_10m.jp2") | ||
|
|
||
| blue = rgb_files[0] | ||
| green = rgb_files[1] | ||
| red = rgb_files[2] |
There was a problem hiding this comment.
If rgb_files is supposed to have exactly three files you can unpack all of them in one go:
| rgb_files = glob.glob("R10m/*B0[2-4]_10m.jp2") | |
| blue = rgb_files[0] | |
| green = rgb_files[1] | |
| red = rgb_files[2] | |
| blue, green, red = glob.glob("R10m/*B0[2-4]_10m.jp2") |
| import subprocess | ||
|
|
||
|
|
||
| # user definated variables |
| TILE_INDEX = "nz-150k-tile-index.shp" | ||
| MIN = "800" | ||
| MAX = "3000" | ||
| CWD = "C:/Temp/image-processing/test" |
There was a problem hiding this comment.
"CWD" ("current working directory") has a specific meaning in operating systems, being the directory from which a program was executed (not the directory where the script lives). You can get this directory using os.getcwd().
| MAX = "3000" | ||
| CWD = "C:/Temp/image-processing/test" | ||
| SENTINEL_IMAGE_DIR = "R10m" | ||
| GDAL_DIR = "C:/Program Files/QGIS 3.34.4/apps/Python39/Scripts" |
There was a problem hiding this comment.
This is very brittle, because the directory path will be different for anyone using a different OS or version of QGIS. This can be made more flexible in a couple of ways:
- Pass the directory as an option, like
--gdal-dir="C:/Program Files/QGIS 3.34.4/apps/Python39/Scripts". - If files from that directory are discoverable by the surrounding shell, that is, if you can run for example
gdal_merge.pywithout specifying its directory, then you can useshutil.which("gdal_merge.py")to get its absolute path.
| if os.path.exists(path): | ||
| shutil.rmtree(path) |
There was a problem hiding this comment.
This is very risky! It would be pretty easy for a programming mistake to delete a bunch of files which shouldn't be deleted. I'd recommend either:
- raising an exception telling the user what to do, or
- using
tempfile.mkdtempto create a new random directory every run and tell the user which directory to look for their data in.
| """ | ||
| Create directors for the various processed images to be saved into. | ||
| """ | ||
| path = os.path.join(CWD, dir_name) |
There was a problem hiding this comment.
If you from os.path import join you can shorten these expressions. But this is a stylistic choice, so I'll leave it to you whether you want to do that.
| print("Merging image tiles into one RGB image") | ||
| subprocess.run( | ||
| [ | ||
| "python.exe", |
There was a problem hiding this comment.
If you run just "python" instead, it should work on other platforms without changes. But this line shouldn't even be necessary, since gdal_merge.py should be executable.
| ogr_info = subprocess.Popen( | ||
| ["ogrinfo", "-ro", "-al", TILE_INDEX], stdout=subprocess.PIPE | ||
| ) | ||
|
|
||
| # tile into Topo50 sheet extents | ||
| ogr_output = ogr_info.communicate()[0].decode("utf-8") |
There was a problem hiding this comment.
It looks like you run this and then just wait for it to finish immediately. In that case I'd recommend ogr_output = subprocess.run([…], capture_output=True, check=True, encoding="utf-8").output.
| MIN = "800" | ||
| MAX = "3000" |
There was a problem hiding this comment.
These are passed to gdal_translate -scale; based on its documentation, would it be useful to name them something like MIN_SCALE and MAX_SCALE?
| "1", | ||
| "255", |
There was a problem hiding this comment.
According to the documentation, "If omitted the output range is 0 to 255". Do we actually want to use the output range 1 to 255 instead?
| import os | ||
| import glob | ||
| import shutil | ||
| from osgeo import gdal |
There was a problem hiding this comment.
It doesn't look like you're using this import anywhere. Did you mean to change the subprocess uses below to use this library? That would be preferable.
Python script and Topo50 index sheets shapefile for processing Sentinel 2 imagery.