Skip to content

Instantly share code, notes, and snippets.

@ashnair1
Last active November 1, 2022 06:29
Show Gist options
  • Save ashnair1/316ba2289b463b013f7304632767aa00 to your computer and use it in GitHub Desktop.
Save ashnair1/316ba2289b463b013f7304632767aa00 to your computer and use it in GitHub Desktop.

Geospatial

Contents

Software Installation Notes

SNAP

  1. Create a conda environment as follows:
conda create -n sar -c conda-forge pyrosar python=3.6
  1. Install the SNAP toolbox from here.

For Ubuntu download the shell script, and simply execute it. When you are asked to configure your python interpreter, specify the python from the conda environment created above.

  1. Navigate to your snap installation dir
$ cd <snap-install-dir>/bin
# Run the following command to configure your python to use the SNAP Python interface
# ./snappy-conf <python-exe> <snappy-dir>

# Mine was as follows:
$ ./snappy-conf /home/ashwin/anaconda3/envs/sar/bin/python /home/ashwin/anaconda3/envs/sar/lib/python3.6/site-packages/
  1. Test imports
$ conda activate sar
$ python
Python 3.6.12 | packaged by conda-forge | (default, Dec  9 2020, 00:36:02) 
[GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import snappy
INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: Incompatible GDAL 3.2.1 found on system. Internal GDAL 3.0.0 from distribution will be used.
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: Internal GDAL 3.0.0 set to be used by SNAP.
INFO: org.esa.snap.core.util.EngineVersionCheckActivator: Please check regularly for new updates for the best SNAP experience.

For MacOS

Installation of SNAP is easy i.e. download the .dmg file and install it. However, configuring snap to work with python is a bit tricky.

  1. Install the SNAP toolbox. It will be located in /Applications/snap
  2. Check whether you have Java installed and see if $JAVA_HOME is configured. If not, install this version from Oracle's website.
  3. Install jpy
# install jpy
git clone https://github.com/bcdev/jpy.git
cd jpy/
python setup.py bdist_wheel
cp dist/*.whl "/Users/ashwinnair/anaconda3/envs/sar/lib/python3.6/site-packages/snappy"
  1. Try running
./snappy-conf /Users/ashwinnair/anaconda3/envs/sar/bin/python /Users/ashwinnair/anaconda3/envs/sar/lib/python3.6/site-packages

This step might work but most likely you'll run into something like this:

Configuring SNAP-Python interface...
java.io.IOException: Python configuration failed.
Command [/Users/ashwinnair/anaconda3/envs/sar/bin/python ./snappyutil.py --snap_home /Applications/snap --java_module /Applications/snap/snap/modules/org-esa-snap-snap-python.jar --force --log_file ./snappyutil.log --jvm_max_mem 11G --java_home /Applications/snap/.install4j/jre.bundle/Contents/Home --req_arch x86_64]
failed with return code 97.
Please check the log file '/Users/ashwinnair/anaconda3/envs/sar/lib/python3.6/site-packages/snappy/snappyutil.log'.
	at org.esa.snap.python.PyBridge.configureJpy(PyBridge.java:232)
	at org.esa.snap.python.PyBridge.installPythonModule(PyBridge.java:149)
	at org.esa.snap.rcp.cli.SnapArgsProcessor.processPython(SnapArgsProcessor.java:103)
	at org.esa.snap.rcp.cli.SnapArgsProcessor.process(SnapArgsProcessor.java:49)
	at org.netbeans.modules.sendopts.DefaultProcessor.process(DefaultProcessor.java:202)
	at org.netbeans.spi.sendopts.Option$1.process(Option.java:387)
	at org.netbeans.api.sendopts.CommandLine.process(CommandLine.java:317)
	at org.netbeans.modules.sendopts.HandlerImpl.execute(HandlerImpl.java:62)
	at org.netbeans.modules.sendopts.Handler.cli(Handler.java:69)
	at org.netbeans.CLIHandler.notifyHandlers(CLIHandler.java:234)
	at org.netbeans.core.startup.CLICoreBridge.cli(CLICoreBridge.java:82)
	at org.netbeans.CLIHandler.notifyHandlers(CLIHandler.java:234)
	at org.netbeans.CLIHandler$1.exec(CLIHandler.java:268)
	at org.netbeans.CLIHandler.finishInitialization(CLIHandler.java:447)
	at org.netbeans.MainImpl.finishInitialization(MainImpl.java:256)
	at org.netbeans.Main.finishInitialization(Main.java:92)
	at org.netbeans.core.startup.Main.start(Main.java:316)
	at org.netbeans.core.startup.TopThreadGroup.run(TopThreadGroup.java:123)
	at java.lang.Thread.run(Thread.java:748)
Python configuration error: Python configuration failed.
Command [/Users/ashwinnair/anaconda3/envs/sar/bin/python ./snappyutil.py --snap_home /Applications/snap --java_module /Applications/snap/snap/modules/org-esa-snap-snap-python.jar --force --log_file ./snappyutil.log --jvm_max_mem 11G --java_home /Applications/snap/.install4j/jre.bundle/Contents/Home --req_arch x86_64]
failed with return code 97.
Please check the log file '/Users/ashwinnair/anaconda3/envs/sar/lib/python3.6/site-packages/snappy/snappyutil.log'.

You'll probably be prompted by macOS itself saying that it needs JRE v6. If you're using a version below Catalina, you could just install it. But the easier solution is to make the following change:

sudo vim /Library/Java/JavaVirtualMachines/jdk1.8.0_281.jdk/Contents/Info.plist

Replace

<key>JVMCapabilities</key>
 <array>
  <string>CommandLine</string>
 </array>

with

<key>JVMCapabilities</key>
 <array>
  <string>JNI</string>
  <string>BundledApp</string>
  <string>CommandLine</string>
 </array>

Refer here for why this change is required.

Running step 4 now should successfully configure the SNAP-Python interface.

Configuring SNAP-Python interface...
Done. The SNAP-Python interface is located in '/Users/ashwinnair/anaconda3/envs/sar/lib/python3.6/site-packages/snappy'
When using SNAP from Python, either do: sys.path.append('/Users/ashwinnair/anaconda3/envs/sar/lib/python3.6/site-packages')
or copy the 'snappy' module into your Python's 'site-packages' directory.

Make sure to add the SNAP binary to path so that pyroSAR detects it.

References:

  1. https://senbox.atlassian.net/wiki/spaces/SNAP/pages/50855941/Configure+Python+to+use+the+SNAP-Python+snappy+interface
  2. https://forum.step.esa.int/t/problem-with-snappy-and-error/13606/5

Orfeo Toolbox

One major requirement,

Your libpython so file corresponding to your python version needs to be located here /usr/lib/x86_64-linux-gnu/ Adding to LD_LIBRARY_PATH doesnt work despite it saying it searched those directories.

So if you're using python3.7, libpython3.7m.so.1.0 needs to be inside /usr/lib/x86_64-linux-gnu/. You might need to recompile python bindings but I think it's not necessary.

It's a similar issue for Macos, except you can't put it in /usr/lib due to this. So pick another local directory from the list

Use locate <.so file> (Linux) or mdfind -name <.so file> (MacOS) to find the library

Once installed, go to QGIS -> Plugins -> Enable OTB. Then configure it by going to Processing -> Options -> Providers -> OTB

Geoserver

Refer this article

Editing GeoJSON files

QGIS allows editing GeoJSON files, with good topology control to snap vertices and add/delete/move nodes, and export directly to GeoJSON. (Oddly, QGIS allows import but not live editing of GeoJSON files. You need to save as in Shapefile format before editing is enabled.)

Image Registration

  • Using Orfeo Toolbox

Steps:

How does FineRegistration work? 3

This application computes a disparity map between image 1 and image 2, i.e. the shift (in pixels) between the pixels of the two images. For each pixel p1 of image 1 the algorithm will search in image 2 which pixel p2 in the neighborhood of p1 is the "closest" to p1. The size of this (rectangular) neighoborhood is determined by the erx and ery parameters (respectively in x and y). To determine how "close" two pixels are, a metric is computed on rectangular patches around the two pixels (for example cross correlation, this is a parameter of the application), the size of these patches is determined by the parameters mrx and mry (size in x and y). This metric is computed for all pixels in the neigborhood, and the corresponding shift is then deduced by finding the extremum of the metric.

Also take a look at 2 for more info

Instalation -> conda install -c conda-forge arosics

Sample command:

$ arosics_cli.py global 2016_test.tif 2014_test.tif -fmt_out GTIFF -max_shift 10
Calculating actual data corner coordinates for reference image...
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Bounding box of calculated footprint for reference image:
	(778094.0, 2567286.0, 787984.0, 2580890.0)
Calculating actual data corner coordinates for image to be shifted...
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Bounding box of calculated footprint for image to be shifted:
	(780670.0, 2567572.0, 788862.0, 2583956.0)
Matching window position (X,Y): 784296.2454573751/2574318.8850265993
Detected integer shifts (X/Y):                            10/4
Detected subpixel shifts (X/Y):                           -0.40409536782327765/-0.4598884162158167
Calculated total shifts in fft pixel units (X/Y):         9.595904632176723/3.5401115837841832
Calculated total shifts in reference pixel units (X/Y):   9.595904632176723/3.5401115837841832
Calculated total shifts in target pixel units (X/Y):      9.595904632176723/3.5401115837841832
Calculated map shifts (X,Y):				  19.191809264360927/-7.080223167780787
Calculated absolute shift vector length in map units:     20.456175178786225
Calculated angle of shift vector in degrees from North:   290.2499761632992
Original map info: ['UTM', 1.0, 1.0, 780670.0, 2583956.0, 2.0, 2.0, 39, 'North', 'WGS-84']
Updated map info:  ['UTM', 1.0, 1.0, '780689.1918092644', '2583948.919776832', 2.0, 2.0, 39, 'North', 'WGS-84']
Image similarity within the matching window (SSIM before/after correction): 0.4406 => 0.6084
Estimated reliability of the calculated shifts:  54.7 %
Correcting geometric shifts...
Warping progress     |==================================================| 100.0% Complete  => 0:00:03
Writing GeoArray of size (8193, 4097, 3) to /home/ashwin/Downloads/2014_test__shifted_to__2016_test.tif.

total processing time: 8.34s

Much easier to use, more specialised

Transform geo to pixel coordinates

Two equivalent ways. f is a DatasetReader

coords = [rowcol(f.transform, x, y) for x, y in geom]
coords = [f.index(x, y) for x, y in geom]

References:

  1. https://forum.orfeo-toolbox.org/t/gaofen-1-image-pansharpening/566
  2. https://gis.stackexchange.com/questions/389097/image-registration-params-using-orfeo-toolbox/389102#389102
  3. http://otb-users.37221.n3.nabble.com/Fine-registration-understanding-parameters-td4030701.html
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment