Buscar este blog

viernes, 6 de octubre de 2017

Some steps for processing Satellite Images with Pipeline

Workshop VDAP:

Go into the sabancaya folder (where all the satellite images) and type the following command (serves also to verify if program is running) :-)

$ wv_correct

To load an image file:

$ wv_correct 16sep10151208-p1bs_r1c1-500885382010_01_p001.ntf 16sep10151208-p1bs_r1c1-500885382010_01_p001.xml 16sep10151208-p1bs_r1c1-500885382010_01_p001.tif 

Wee need to do this 4 times, for each NTF file and its corresponding XML file. This will generate new TIF files

16sep10151208-p1bs_r1c1-500885382010_01_p001.ntf
16sep10151208-p1bs_r2c1-500885382010_01_p001.ntf
16sep10151321-p1bs_r1c1-500885382010_01_p001.ntf
16sep10151321-p1bs_r2c1-500885382010_01_p001.ntf

16sep10151208-p1bs_r1c1-500885382010_01_p001.xml
16sep10151208-p1bs_r2c1-500885382010_01_p001.xml
16sep10151321-p1bs_r1c1-500885382010_01_p001.xml
16sep10151321-p1bs_r2c1-500885382010_01_p001.xml


Then we have to run the following command for the newTIF files generated:

$ dg_mosaic

$ dg_mosaic 16sep10151208*tif --output-prefix 16sep10151208

This will generate a single TIF file which is composed from the previous 2 TIF files (16sep10151208*tif). Then we do the same for the other couple of files (16sep10151321*tif) as follows:

$ dg_mosaic 16sep10151321*tif --output-prefix 16sep10151321

The next step is absolutely necessary!!!! some manuals says it is not true!!! ;-)

We need to download the DEM global coverage from earth explorer.gov for our area of work.
SRTM version 3 is the preferred DEM dataset. We can use also other source such as ASTER or LIDAR.

Be careful with the next step, If you have installed GMT so the command map project which also exist on GMT and Pipeline do not the same stuff. So here I have

$ /usr/local/geodesy/stereopipeline/StereoPipeline-2.5.3/bin/mapproject -t rpc --t_srs EPSG:32718 --tr 0.5 Sabancaya_utm.tif 16sep10151208.r100.tif 16sep10151208.r100.xml 16sep10151208_map.tif 


From the software directory copy the file which is into the example directory called "stereo.default.example", this file should be copied within the local directory, in this case is the "Sabancaya2016Sep10" directory. Here rename the file as follows:

mv stereo.default.example stereo.default

Open this file and edit as your prefer:
*****
# -*- mode: sh -*-

# Pre-Processing / stereo_pprc
################################################################

# Pre-alignment options
#
# Available choices are (however not all are supported by all sessions):
#    NONE           (Recommended for anything map projected)
#    EPIPOLAR       (Recommended for Pinhole Sessions)
#    HOMOGRAPHY     (Recommended for ISIS wide-angle shots)
#    AFFINEEPIPOLAR (Recommended for ISIS narrow-angle and DG sessions)
alignment-method none

# Intensity Normalization
force-use-entire-range       # Use entire input range

# Select a preprocessing filter:
#
# 0 - None
# 1 - Subtracted Mean
# 2 - Laplacian of Gaussian (recommended)
prefilter-mode 2

# Kernel size (1-sigma) for pre-processing
#
# Recommend 1.4 px for Laplacian of Gaussian
# Recommend 25 px for Subtracted Mean
prefilter-kernel-width 1.4

# Integer Correlation / stereo_corr
################################################################

# Select a cost function to use for initialization:
#
# 0 - absolute difference (fast)
# 1 - squared difference  (faster .. but usually bad)
# 2 - normalized cross correlation (recommended)
cost-mode 2

# Initialization step: correlation kernel size
corr-kernel 21 21

# Initializaion step: correlation window size
# corr-search -80 -2 20 2

# Subpixel Refinement / stereo_rfne
################################################################

# Subpixel step: subpixel modes
#
# 0 - disable subpixel correlation (fastest)
# 1 - parabola fitting (draft mode - not as accurate)
# 2 - affine adaptive window, bayes EM weighting (slower, but much more accurate)
# 3 - affine window, (intermediate speed, results similar to bayes EM)
subpixel-mode 2  (THIS HAS BEEN CHANGED)!!!!!

# Subpixel step: correlation kernel size
subpixel-kernel 21 21

# Post Filtering / stereo_fltr
################################################################

# Fill in holes up to 100,000 pixels in size with an inpainting method
# disable-fill-holes

# Automatic "erode" low confidence pixels
filter-mode 1
rm-half-kernel 5 5
max-mean-diff 3
rm-min-matches 60
rm-threshold 3
rm-cleanup-passes 1

# Triangulation / stereo_tri
################################################################

# Size max of the universe in meters and altitude off the ground.
# Setting both values to zero turns this post-processing step off.
near-universe-radius 0.0
far-universe-radius 0.0

****

run the following command:
$ stereo -t dgmaprpc --subpixel-mode 2 --alignment-method none 16sep10151208_map.tif 16sep10151321_map.tif 16sep10151208.r100.xml 16sep10151321.r100.xml SABANCAYA/SABANCAYA Sabancaya_utm.tif 

On oct 9 2017

$ point2dem

$ point2dem --nodata-value -9999 --t_srs EPSG:32718 -s 2.0 --dem-hole-fill-len 100 SABANCAYA/SABANCAYA-PC.tif --orthoimage-hole-fill-len 100 SABANCAYA/SABANCAYA-L.tif --errorimage

If an error message appears (No ortho image was requested, yet texture files were passed as inputs) just take out the orthophoto option:
$point2dem --nodata-value -9999 --t_srs EPSG:32718 -s 2.0 --dem-hole-fill-len 100 SABANCAYA/SABANCAYA-PC.tif  --errorimage

If your are using ARGIS just open the TIF output file and convert it to a hillshade.

If not, to create a HILLSHADE using command line:
hillshade SABANCAYA/SABANCAYA-DEM.tif -o SABANCAYA/SABANCAYA-HLLSHD.tif -e 45 -a 315

Then, open it in QGIS 

PROCESSING IMAGED FROM 2014, in order to compare both epochs
====
Angie provide us with a satellite image of Sabancaya dating back to 2014.  The name of the Folder is "2014Sabancaya28July" inside the folder there are several files as follows: 
450M Aug  8  2014 14JUL28150826-P1BS_R03C1-500088628040_04_P001.NTF
707K Jul 28  2014 14JUL28150826-P1BS_R03C1-500088628040_04_P001.XML
452M Oct  9  2017 14JUL28150826-P1BS_R04C1-500088628040_04_P001.NTF
707K Jul 28  2014 14JUL28150826-P1BS_R04C1-500088628040_04_P001.XML
373M Oct  9  2017 14JUL28150912-P1BS_R03C1-500088628040_04_P001.NTF
662K Jul 28  2014 14JUL28150912-P1BS_R03C1-500088628040_04_P001.XML
374M Oct  9  2017 14JUL28150912-P1BS_R04C1-500088628040_04_P001.NTF
662K Jul 28  2014 14JUL28150912-P1BS_R04C1-500088628040_04_P001.XML

Then copy from previous/older directory the following files:

450M Aug  8  2014 14JUL28150826-P1BS_R03C1-500088628040_04_P001.NTF*
707K Jul 28  2014 14JUL28150826-P1BS_R03C1-500088628040_04_P001.XML*
452M Oct  9  2017 14JUL28150826-P1BS_R04C1-500088628040_04_P001.NTF*
707K Jul 28  2014 14JUL28150826-P1BS_R04C1-500088628040_04_P001.XML*
373M Oct  9  2017 14JUL28150912-P1BS_R03C1-500088628040_04_P001.NTF*
662K Jul 28  2014 14JUL28150912-P1BS_R03C1-500088628040_04_P001.XML*
374M Oct  9  2017 14JUL28150912-P1BS_R04C1-500088628040_04_P001.NTF*
662K Jul 28  2014 14JUL28150912-P1BS_R04C1-500088628040_04_P001.XML*
 92B Oct  9 14:41 Sabancaya_utm.tfw*
 15M Oct  9 14:41 Sabancaya_utm.tif*
2.3K Oct  9 14:41 Sabancaya_utm.tif.aux.xml*
5.0M Oct  9 14:41 Sabancaya_utm.tif.ovr*
5.6K Oct  9 14:41 Sabancaya_utm.tif.xml*
2.3K Oct  9 14:41 stereo.default
Then run the WV_CORRECT program in order to generate the tif files. Do this for each of the four couple of NTF and XML files.


$ wv_correct 14JUL28150826-P1BS_R03C1-500088628040_04_P001.NTF 14JUL28150826-P1BS_R03C1-500088628040_04_P001.XML 14JUL28150826-P1BS_R03C1-500088628040_04_P001.tif

$ wv_correct 14JUL28150826-P1BS_R04C1-500088628040_04_P001.NTF 14JUL28150826-P1BS_R04C1-500088628040_04_P001.XML 14JUL28150826-P1BS_R04C1-500088628040_04_P001.tif

$ wv_correct 14JUL28150912-P1BS_R03C1-500088628040_04_P001.NTF 14JUL28150912-P1BS_R03C1-500088628040_04_P001.XML 14JUL28150912-P1BS_R03C1-500088628040_04_P001.tif


$ wv_correct 14JUL28150912-P1BS_R04C1-500088628040_04_P001.NTF 14JUL28150912-P1BS_R04C1-500088628040_04_P001.XML 14JUL28150912-P1BS_R04C1-500088628040_04_P001.tif

This is an ls to the directory:
2.9K Oct  9 14:44 14JUL28150826-P1BS_R03C1-500088628040_04_P001.IMD
450M Aug  8  2014 14JUL28150826-P1BS_R03C1-500088628040_04_P001.NTF*
707K Jul 28  2014 14JUL28150826-P1BS_R03C1-500088628040_04_P001.XML*
3.4G Oct  9 14:47 14JUL28150826-P1BS_R03C1-500088628040_04_P001.tif
2.9K Oct  9 14:48 14JUL28150826-P1BS_R04C1-500088628040_04_P001.IMD
452M Oct  9  2017 14JUL28150826-P1BS_R04C1-500088628040_04_P001.NTF*
707K Jul 28  2014 14JUL28150826-P1BS_R04C1-500088628040_04_P001.XML*
3.4G Oct  9 14:51 14JUL28150826-P1BS_R04C1-500088628040_04_P001.tif
2.9K Oct  9 14:52 14JUL28150912-P1BS_R03C1-500088628040_04_P001.IMD
373M Oct  9  2017 14JUL28150912-P1BS_R03C1-500088628040_04_P001.NTF*
662K Jul 28  2014 14JUL28150912-P1BS_R03C1-500088628040_04_P001.XML*
2.9G Oct  9 14:54 14JUL28150912-P1BS_R03C1-500088628040_04_P001.tif
2.9K Oct  9 14:54 14JUL28150912-P1BS_R04C1-500088628040_04_P001.IMD
374M Oct  9  2017 14JUL28150912-P1BS_R04C1-500088628040_04_P001.NTF*
662K Jul 28  2014 14JUL28150912-P1BS_R04C1-500088628040_04_P001.XML*
2.9G Oct  9 14:56 14JUL28150912-P1BS_R04C1-500088628040_04_P001.tif
 92B Oct  9 14:41 Sabancaya_utm.tfw*
 15M Oct  9 14:41 Sabancaya_utm.tif*
2.3K Oct  9 14:41 Sabancaya_utm.tif.aux.xml*
5.0M Oct  9 14:41 Sabancaya_utm.tif.ovr*
5.6K Oct  9 14:41 Sabancaya_utm.tif.xml*
2.3K Oct  9 14:41 stereo.default


Then run the dg_mosaic for each of the TIF timestamp 14JUL281508  & 14JUL281509
Will create a new mosaic image with its corresponding camera model/parameters.

$dg_mosaic 14JUL28150826-P1BS_R0*.tif --output-prefix 14JUL28150826 

$ dg_mosaic 14JUL28150912-P1BS_R0*.tif --output-prefix 14JUL28150912

These are the files that have been generated:

2.9K Oct  9 14:44 14JUL28150826-P1BS_R03C1-500088628040_04_P001.IMD
450M Aug  8  2014 14JUL28150826-P1BS_R03C1-500088628040_04_P001.NTF*
707K Jul 28  2014 14JUL28150826-P1BS_R03C1-500088628040_04_P001.XML*
3.4G Oct  9 14:47 14JUL28150826-P1BS_R03C1-500088628040_04_P001.tif
2.9K Oct  9 14:48 14JUL28150826-P1BS_R04C1-500088628040_04_P001.IMD
452M Oct  9  2017 14JUL28150826-P1BS_R04C1-500088628040_04_P001.NTF*
707K Jul 28  2014 14JUL28150826-P1BS_R04C1-500088628040_04_P001.XML*
3.4G Oct  9 14:51 14JUL28150826-P1BS_R04C1-500088628040_04_P001.tif
5.2G Oct  9 15:04 14JUL28150826.r100.tif
703K Oct  9 15:00 14JUL28150826.r100.xml
2.9K Oct  9 14:52 14JUL28150912-P1BS_R03C1-500088628040_04_P001.IMD
373M Oct  9  2017 14JUL28150912-P1BS_R03C1-500088628040_04_P001.NTF*
662K Jul 28  2014 14JUL28150912-P1BS_R03C1-500088628040_04_P001.XML*
2.9G Oct  9 14:54 14JUL28150912-P1BS_R03C1-500088628040_04_P001.tif
2.9K Oct  9 14:54 14JUL28150912-P1BS_R04C1-500088628040_04_P001.IMD
374M Oct  9  2017 14JUL28150912-P1BS_R04C1-500088628040_04_P001.NTF*
662K Jul 28  2014 14JUL28150912-P1BS_R04C1-500088628040_04_P001.XML*
2.9G Oct  9 14:56 14JUL28150912-P1BS_R04C1-500088628040_04_P001.tif
4.5G Oct  9 15:09 14JUL28150912.r100.tif
659K Oct  9 15:06 14JUL28150912.r100.xml
 92B Oct  9 14:41 Sabancaya_utm.tfw*
 15M Oct  9 14:41 Sabancaya_utm.tif*
2.3K Oct  9 14:41 Sabancaya_utm.tif.aux.xml*
5.0M Oct  9 14:41 Sabancaya_utm.tif.ovr*
5.6K Oct  9 14:41 Sabancaya_utm.tif.xml*
2.3K Oct  9 14:41 stereo.default

Now we generate a map projected file using as initial reference the SRTM 30 meters  file for that area. 
$ /usr/local/geodesy/stereopipeline/StereoPipeline-2.5.3/bin/mapproject -t rpc --t_srs EPSG:32718 --tr 0.5 Sabancaya_utm.tif 14JUL28150826.r100.tif 14JUL28150826.r100.xml 14JUL28150826_map.tif 

$ /usr/local/geodesy/stereopipeline/StereoPipeline-2.5.3/bin/mapproject -t rpc --t_srs EPSG:32718 --tr 0.5 Sabancaya_utm.tif 14JUL28150912.r100.tif 14JUL28150912.r100.xml 14JUL28150912_map.tif



Using GUI  to zoom in and out 
$ stereo_gui -t dgmaprpc --subpixel-mode 2 --alignment-method none 14JUL28150826_map.tif 14JUL28150912_map.tif 14JUL28150826.r100.xml 14JUL28150912.r100.xml 2014SABANCAYA/2014SABANCAYA Sabancaya_utm.tif

Once this process ends, a window showing the two images will appear. Here we need to choose the interest area, a small area we want to study.

Then the command 

$point2dem


and finally the HILLSHADE
hillshade 

lastly open it on QGIS


########  10 october 2017  #####
Go into the folder  /stereopipeline/Sabancaya2016Sep10

Georeference

pc_align --max-displacement 10 --alignment-method point-to-point --csv-proj4 EPSG:32718 --csv-format '1:easting 2:northing 3:height_above_datum' --datum WGS_1984 SABANCAYA/SABANCAYA-DEM.tif igp_utm.csv --save-inv-transformed-reference-points --output-prefix SABANCAYA_ref/SABANCAYA_ref

pc_align --max-displacement 10 --alignment-method point-to-point --csv-proj4 EPSG:32718 --csv-format '1:easting 2:northing 3:height_above_datum' --datum WGS_1984 igp_utm.csv SABANCAYA/SABANCAYA-DEM.tif --save-transformed-source-points --output-prefix SABANCAYA_ref/SABANCAYA_ref

No hay comentarios:

Publicar un comentario