Merge pull request #1 from griembauer/initial_keyhole_workflow

grass-keyhole: add addons,scripts,readmes

README.md

@@ -0,0 +1,42 @@
+# Overview
+## Summary
+This repository contains GRASS GIS Addons, bash scripts, and a detailed step-by-step documentation of a GRASS GIS based georectification process of satellite images from the KEYHOLE reconnaissance missions KH-4A/B (CORONA), KH-7 (GAMBIT), and KH-9 (HEXAGON). Further background can be found [here](https://www.mundialis.de/en/georeferenzierung-von-corona-spionagesatellitendaten/).
+## Repo Content
+ - **documentation**: contains a detailed step-by-step guide to rectify KEYHOLE reconnaissance missions based on GRASS GIS. See [documentation/README.md](documentation/README.md) to get started.
+ - **grass_addons**: contains two addons for GRASS GIS that are required for the rectification process
+ - **scripts**: contains shell scripts that are required for the rectification process
+
+
+ ## KEYHOLE missions overview:
+
+ ##### KH-4A, KH-4B - CORONA
+
+ - Panorama camera
+     - -> orthorectification with panorama correction
+ - Film size: 70 mm x 756.9 mm
+ - flight direction seems to be south
+ - in the scanned film, north is approximately up
+
+ ##### KH-7 - GAMBIT
+
+ - strip camera, acquired imagery in continuous lengthwise sweeps of the terrain.
+     - -> orthorectification does not really work, bad geolocation
+     - -> panorama correction not applicable
+     - -> standard rectification
+ - Film size: 9 inch (228.6 mm) x variable length
+ - flight direction seems to be south
+ - in the scanned film, north is approximately right
+
+ ##### KH-9 - HEXAGON
+
+ - Panorama camera
+     -> orthorectification with panorama correction
+ - flight direction seems to be south
+ - in the scanned film, north is approximately up
+ - focal_length: 152.4 cm
+
+ ##### Further KEYHOLE missions
+
+ - KH-6 Lanyard: Lanyard was a failed mission - no data with acceptable quality available
+ - KH-5 Argon: ground resolution too low - irrelevant
+ - KH-8 Gambit 3: no data available

documentation/0_preparation.md

@@ -0,0 +1,70 @@
+# Preparation
+## Operating System
+The described process was developed and tested using Linux Mint 20. It is in general optimized for use with Linux. No tests on other operating systems were performed.
+
+## Required Software
+- [GRASS GIS >= v7.9](https://grasswiki.osgeo.org/wiki/Installation_Guide)
+- [QGIS](https://www.qgis.org/en/site/forusers/alldownloads.html#)
+- (optional): [Hugin](https://hugin.sourceforge.io/) for the automatic stitching process:
+    - for Ubuntu: use `apt-get install`
+            --> `apt-get install hugin`
+    - for Fedora: use `dnf install`
+            --> `dnf install hugin`
+
+- GRASS GIS Addons:
+    - Start GRASS GIS and install the GRASS modules [i.ortho.position](../grass_addons/i.ortho.position) and [i.ortho.corona](../grass_addons/i.ortho.corona) from this repository:
+    ```
+    g.extension extension=i.ortho.position url=/path/to/this/repo/grass_addons/i.ortho.position
+    g.extension extension=i.ortho.corona url=/path/to/this/repo/grass_addons/i.ortho.corona
+    ```
+    - install the [r.in.nasadem](https://grass.osgeo.org/grass82/manuals/addons/r.in.nasadem.html) addon from the official GRASS Addons repository:
+    ```
+    g.extension extension=r.in.nasadem
+    ```
+    To use this addon, you will need a registered user account at [NASA EARTHDATA](https://urs.earthdata.nasa.gov/)
+
+
+## Data preparation
+**0)** Download the raw Keyhole Scene you would like to process and the vector file containing its rough extent and metadata from the USGS Earth Explorer.
+- if the vector dataset is not available yet:
+  - Browse https://earthexplorer.usgs.gov/ and navigate to the scene of interest using the ENTITY ID (in Step 2: "Data Sets" select Declass I for CORONA, Declass II for GAMBIT or Declass III for HEXAGON; enter the ENTITY ID in Step 3: "Additional Criteria"). Click at "Click here to export your results" (above the search results) to export a .zip-file with the vector dataset of the scene
+  - Note: For some HEXAGON scenes there is no specific vector dataset that can be downloaded. Instead the .zip-file will contain footprints for all HEXAGON scenes around the world. To cut out only the desired footprint you need to open the attribute table of this dataset in e.g. QGIS, search for the desired footprint using the ENTITY ID of the scene and create a new layer that only contains this footprint. Save this layer as a vector dataset and continue at the step above.  
+
+**1)** Launch GRASS GIS and create a GRASS xy location, e.g. corona_xy (you can name it e.g. gambit_xy for GAMBIT scenes - Note: in the walkthrough below the location is always called corona_xy)
+- xy means no CRS, no georeferencing information as the original USGS scans do not have any CRS information
+
+**2)** create a GRASS location in the projection of the target UTM zone
+- example for zone 32N:
+    `grass -c epsg:32632 ~/grassdata/corona_utm32n`
+
+**3)** download/import the DEM (Note: only applies to KH-4A/KH-4B CORONA and KH-9 HEXAGON scenes)
+- if the DEM is already available, import it into the UTM location using `r.import input=/path/to/the/DEM.tif output=dem_<country>`
+- if no DEM is available yet:
+   - import the vector dataset of the rough scene extent using `v.import`: `v.import input=/path/to/scene_extent.shp output=<scene>_extent`
+   - set the region to this vector and add some km of buffer around it (due to inaccuracies of the footprints):
+    `g.region vector=<scene>_extent n=n+60000 s=s-60000 w=w-60000 e=e+60000 res=30 -pa`
+   - import the nasadem:
+    `r.in.nasadem user=<USER> password=<PASSWORD> output=dem_<COUNTRY> memory=6000 resolution=30 method=bilinear`
+Note: The NASA only provides DEMs that cover an area of the Earth that consists at least proportionally of land mass. All parts of the earth that lie in water are not covered by the DEMs. Due to the enlargement of the region by 60km, it may happen that `r.in.nasadem` wants to download DEMs that are completely in water, especially when the area of interest is a coastal region - these DEMs are not available! Consequently, `r.in.nasadem` will raise an error. In order to get the DEMs anyway, either the region in the corresponding cardinal direction must be set smaller than 60km, or the DEMs must be downloaded manually for the corresponding area under https://search.earthdata.nasa.gov/.
+   - set DEM NULLs to 0 to avoid nodata in oceans:
+    `r.null dem_<COUNTRY> null=0`
+   - round the DEM to integer values to save disk space:
+    ` r.mapcalc expression="dem_<country> = round(dem_<country>)" --o`
+
+Note: In the walkthrough below, it is assumed that the target UTM Zone is 32N.
+For other zones, this has to be adapted in steps 2_gcp_collection, 3_georectification and 3_orthorectification
+
+**4)** collect scene specifications
+
+A scene specification document is not required as input, but it is helpful to keep track of the processing in it.
+1. Get the four approximate corner coordinates.
+    - Use the scene extent vector dataset
+    - OR, check the USGS earth explorer (https://earthexplorer.usgs.gov/):
+        - select Dataset -> Declassifed -> declassI (for CORONA), declassII (for GAMBIT) or declassIII (for HEXAGON)
+        - additional criteria: scene id
+        - get metadata + browse overlay
+2. Check if the image needs to be flipped --> use the USGS earth explorer thumbnails to find out
+3. For CORONA and HEXAGON note the look direction of the camera:
+    - It's indicated in the scene name (F for forward, A for aft)
+
+**Continue** with the stitching of the individual scene scans in [1_stitching.md](1_stitching.md) for manual stitching or [1_automatic_stitching.md](1_automatic_stitching.md) for automatic stitching using the Hugin panorama software

documentation/1_automatic_stitching.md

@@ -0,0 +1,33 @@
+# Automatic Stitching
+The Keyhole scenes are delivered in several pieces: CORONA scenes usually in 2-4 pieces, GAMBIT in 2 pieces and HEXAGON in 4-14 pieces. These pieces are not georeferenced and must be stitched together.
+
+Automatic stitching is successfully tested for CORONA (KH-4A/KH-4B) scenes, GAMBIT 1 (KH-7) scenes and HEXAGON (KH-9) scenes with a maximum of 4 pieces. Scenes with a very high proportion of undifferentiable land cover (>80%), e.g. deserts or ocean, or scenes with poor image quality may not be successfully stitched.
+Automatic stitching is particularly suitable for scenes that can also be easily stitched manually and scenes that do not have clearly identifiable infrastructure but do have clearly recognisable land structures such as rock formations.
+
+The automatic stitching is processed with [Hugin](https://hugin.sourceforge.io/).
+
+
+## Requirements
+
+- 30++ GB RAM for CORONA scenes
+- 55++ GB RAM for HEXAGON scenes
+- 65++ GB RAM for GAMBIT scenes
+
+Note: some HEXAGON scenes are delivered in more than 4 pieces. The automatic stitching process for these scenes would possibly need more than 60 GB RAM, which has not been successfully tested yet. Splitting scenes with e.g. 7 pieces in two stitching groups has not been successful too due to the same RAM problem. For those scenes it is possibly practical to enlarge the physical RAM with a swap file to a total of > 80 GB.
+
+## Stitching process
+
+**1)** Create a folder in your local file system that only contains the [hugin_automatic_stitching.sh](../scripts/hugin_automatic_stitching.sh) script and the pieces of the scene that you would like to stitch together. Note: The script has to be adapted for HEXAGON scenes with more or fewer pieces (currently for scene pieces `_a` to `_f`) --> lines 35, 45 in the script
+
+**2)** Navigate to the created folder in your terminal and run the `hugin_automatic_stitching.sh` script passing the name of the satellite mission (`CORONA`, `GAMBIT` or `HEXAGON`) as an argument:
+
+e.g. `bash hugin_automatic_stitching.sh CORONA`
+
+It outputs `stitched_D*.tif` which is the result of the stitching process.
+
+**3)** The stitching process produces a `.tif` that contains 2 bands. Only the first band is relevant for further processing steps. Therefore only import the first band to the GRASS xy location using the `band=1`parameter.
+
+`r.in.gdal input=stitched_D*.tif output=stitched_D* band=1`
+Note: Make sure the output name does not include `-`, otherwise `r.mapcalc` in the next step  will not work. Use the `g.rename` command if required.
+
+**Continue** with the cropping of the remaining physical film border in [2_cropping.md](2_cropping.md).

documentation/1_stitching.md

@@ -0,0 +1,48 @@
+# Stitching
+The Keyhole scenes are delivered in several pieces: CORONA scenes usually in 2-4 pieces, GAMBIT in 2 pieces and HEXAGON in 4-14 pieces. These pieces are not georeferenced and must be stitched together:
+
+**1)** Start GRASS GIS in the xy location that was created in the first step, select the mapset `PERMANENT`
+
+**2)** import the westernmost part using `r.in.gdal` (usually *_d for CORONA, but check with the USGS earth explorer thumbnail):
+   `r.in.gdal in=DS1043-2201DF006_d.tif out=DS1043-2201DF006_d`
+
+**3)** Move the remaining parts to the correct place, while the westermost part stays the reference. For each part after the westernmost, do the following:
+
+ * create a new mapset named after the part, e.g. *_c:
+    `g.mapset -c DS1043-2201DF006_c`
+ * import the part to stitch, e.g. *_c:
+     `r.in.gdal in=DS1043-2201DF006_c.tif out=DS1043-2201DF006_c`
+ * create an image group named after and consisting of the raster, e.g. *_c:
+     `i.group group=DS1043-2201DF006_c in=DS1043-2201DF006_c`
+ * In the GUI: Launch the "Georectify" dialog. source location: corona_xy,
+     source mapset: *_c, e.g. DS1043-2201DF006_c. click next, keep the settings
+     in the next dialog, click next, source-map: the map to be stitched, e.g. *_c in
+     the corresponding mapset. target map: The previous reference map, e.g. *_d in the
+     PERMANENT mapset.
+ * Collect 4-8 GCPs that you find in the overlap area of the two scene pieces, e.g. four at the edge, and two within the images to compensate for rotation. Use distinct pixels.
+ * Under "Georectifier Settings" --> "Georectification", make sure 1.st order is ticked
+ * Tick all GCPs, and click on "Recalculate RMS error". Since we don't have a CRS,
+     the unit of the errors is pixels. It should be < 1.0.
+ * click on "Georectify"
+ * Sometimes the writing to the PERMANENT mapset does not work, if so: copy the "rectified"
+     raster to the PERMANENT mapset and if necessary rename it, e.g.:
+     `g.mapset PERMANENT`
+     `g.copy DS1043-2201DF006_c_georect17884@DS1043-2201DF006_c,DS1043-2201DF006_c`
+ * Load the result in the GUI check if it fits. Therefore, right-click on the reference map, e.g. *_d and choose "Display layer". Then navigate in the GUI to "Layers" (at the bottom) and click "Add Raster Map" and choose the georectified raster.
+ * Delete the temporary mapset in your local file system inside the grassdb
+ * Repeat step **3)** for the next pieces of the scene (e.g. *_b), until all pieces have been moved to the correct position.
+
+**4)** Run the script [raster_patch_xy.sh](../scripts/raster_patch_xy.sh). It patches the pieces together and outputs a raster in the GRASS location/mapset, as well as a GeoTiff and a .vrt with a rough georeferencing.
+* modify in the script the parameters `ORDEREDLIST` (line 21) and `FLIP` (line 28) according to the respective scene. `ORDEREDLIST` has to begin with the part that was stitched last (e.g. *_a, see the comments in the script)
+* run the script by passing the parameters of the scene ID and the corner coordinates:
+
+   `bash raster_patch_xy.sh SCENEID NWLON NWLAT NELON NELAT SWLON SWLAT SELON SELAT`
+
+   e.g.:
+
+   `bash raster_patch_xy.sh DS1043-2201DF006 7.639 36.330 10.720 36.792 7.642 36.116 10.816 36.591`
+
+   (latitude and longitude coordinates can be found in the attribute table of the vector map)
+   Note: KH-7 Gambit scenes may be flipped by rougly 90 degrees so the coordinates passed to the script may not be correct. This however only affects the output .vrt, which is just a visual aid for GCP collection and not crucial for the process.
+
+**Continue** with the cropping of the remaining physical film border in [2_cropping.md](2_cropping.md).

documentation/2_cropping.md

@@ -0,0 +1,32 @@
+# Cropping
+The scenes are delivered with black borders that should be removed for asthetic and functional reasons.
+
+**0)** If you manually stitched the raster: import the stitched raster using `r.import`, e.g.:
+    `r.import in=DS1036-2155DF050_stitched.tif out=DS1036_stitched_raster` Note: Make sure the output name does not include `-`, otherwise `r.mapcalc` in the next step  will not work. Use the `g.rename` command if required.
+
+**1)** Manually set the region in the GUI to the extent of the actual image (only the relevant pixels, no black border etc.):
+
+Therefore right-click on the stitched layer --> "Display Layer" --> "Different Zoom Settings" --> "Set computational region extent interactively" --> Draw a rectangle only around the image data (leave out the black film border)
+
+OR
+
+Use "Display Layer" --> "Analyze Map" --> "Measure Distance" to measure the distance of the four edges of the raster to the start of actual image data, and then set the region accordingly, e.g.: `g.region n=n-4593 s=s+710 e=e-3075 w=w+12865`
+
+**2)** Use `r.mapcalc` to recalculate the stitched raster map. Black areas (with a pixel value of 0) are recalculated to 1 to use the 0 as NoData later on: `r.mapcalc expression="DS1036_stitched_raster_calc = if(DS1036_stitched_raster == 0,1,DS1036_stitched_raster)"`
+
+**3)** `r.mapcalc` changed the colors of the raster map - set it back to grey scale using `r.colors`.
+    `r.colors map=DS1036_stitched_raster_calc color=grey255`
+
+**4)** Run `r.info` to get the basic information about the raster that are used in the next step (rows & columns of the image):
+    `r.info DS1036_stitched_raster_calc`
+
+**5)** Use the values of the rows and columns from the previous step to reset the image coordinates of the raster map so that the image coordinate with the value 0,0 is in the lower left corner of the image. Use the rows value for `north` and the columns value for `east`:
+    `r.region map=DS1036_stitched_raster_calc s=0 n=<ROWS_VALUE> w=0 e=<COLUMNS_VALUE>`
+
+**6)** Set the region to the raster:
+    `g.region raster=DS1036_stitched_raster_calc`
+
+**7)** Export the cropped raster map as a GeoTIFF:
+    `r.out.gdal -cm overviews=5 createopt="COMPRESS=LZW,PREDICTOR=2,TILED=YES" in=DS1036_stitched_raster_calc out=DS1036_stitched_raster_calc.tif`
+
+**Continue** with the selection of GCPs in QGIS in [3_gcp_collection.md](3_gcp_collection.md)

documentation/3.1_gcps_overlapping_scenes.md

@@ -0,0 +1,16 @@
+## GCPS for overlapping scenes
+If you want to rectify scenes from the same overflight that have overlap areas to each other, there is a way to quicker rectifying them than doing it for each scene individually.
+
+Example: You want to rectify the scenes `DS1036-2155DF050`, `DS1036-2155DF051`, `DS1036-2155DF052` and `DS1036-2155DF053`. Instead of collecting GCPs for each scene individually, you can use the overlap area between two scenes to process them quicker. In the following workflow it is assumed that you want to rectify the scenes in the direction from top to bottom (`DS1036-2155DF050` to `DS1036-2155DF053`) - of course it is also possible the other way round.
+
+**1)** Collect GCPs for the first scene (e.g. `DS1036-2155DF050`). Set 10-15 GCPs at the top of the scene and 10-15 at the bottom. The GCPs at the bottom needs to be in an area of the image that overlaps with the next scene (in this example: `DS1036-2155DF051`). Due to little differences in the width of the scenes, which is a result of the flight direction of the satellite, and due to the possibility of poor image quality in the corners of the image you should make sure that you leave a free space at the left and the right border of the scene where you do not set any GCPs (~2,000 pixels at each side of the image).
+
+**2)** Use the collected GCPs to continue the process at the GCP conversion in `3_gcp_collection.md` for the first scene (`DS1036-2155DF050`).
+
+**3)** After finishing the rectification process of the first scene (`DS1036-2155DF050`), you can start off the GCP collection process for the next scene (in this example: `DS1036-2155DF051`) with already half a set of GCPs: Open the file with the collected GCPs from the first scene (the `.points` file) and delete all points that were set in the top area of the image. Save this file and open it in the QGIS-georectify window of the new scene (here: `DS1036-2155DF051`). You now already have the GCPs that belong to the top area of this scene. Those GCPs already have the correct reference coordinates. The x-image coordinates are also almost correct and give a good orientation about the actual position of the points in the image. Only the y-image coordinates are completely wrong. Using the "move" function, the GCPs in the image can now be moved up to their actual position.
+
+**4)** Afterwards, collect further 10-15 GCPs in the lower image area, which are located in the overlap area to the next scene (in this example: `DS1036-2155DF052`). Use the collected GCPs to process the geo- and orthorectification of the second scene (`DS1036-2155DF051`).
+For the third scene the same is done as described in steps 1 to 4.
+Repeat steps 1 to 4 until you have rectified all the required scenes of the same flyover.
+
+By reusing the points from the previous scene, the processing time of the following scene is significantly reduced. If several scenes of the same flyover are to be rectified, this can lead to great time savings with the same quality of results compared to the individual GCP collection for each scene.