1.11. Phase Change Analysis


Saeid Aminjafari

GeoComput. & ML

2 June 2021

1.11.1. Project description

In this project I downloaded 42 Sentinel-1 images (42 dates i.e. every six days from 24 March 2019 to 25 November 2019) from the Alaska Satellite Facility. Then I calculated the phase change between paired images successively by SARScape software. We want to test if the accumulated phase changes of pixels inside a lake correlate with the water level changes in that lake. In cases with good correlation, this method can be replaced by conventional methods to estimate water level changes in lakes. This method has high spatial and temporal resolution, is cost-efficient and fast.

1.11.1.1. Geocomputain datasets

For the course Geo-computation, I have two sets of data: 1. The phase difference between successive dates (41 maps of phase change!). 2. The quality images of the phase change maps.

1.11.1.2. Data processing and computation procedure:

  1. For each quality image, we mask the pixels with quality less than 0.25 (41 mask.tif images).

  2. Multiply all mask images to get the pixels that in all maps have high quality (maskAll.tif).

  3. Multiply the maskAll.tiff by each phase-change map to get the phase change of high-quality pixels.

  4. Crop the masked phase-change maps around the lake area with a shape file.

  5. Get the values of selected pixels and create a time series of phase change for each pixel (42 dates).

  6. Accumulate the time series of the phase change for each pixel.

  7. Calculate the correlation between the accumulated phase change and the water level of the gauging station.

The quality images are stored in a file named ccAll which has 41 bands i.e. the number of images. The code below shows the inforamtion of this file.

[1]:
!gdalinfo Project/ccAll
Driver: ENVI/ENVI .hdr Labelled
Files: Project/ccAll
       Project/ccAll.hdr
Size is 9219, 4226
Coordinate System is:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137.0,298.257223563]],
    PRIMEM["Greenwich",0.0],
    UNIT["degree",0.0174532925199433]]
Origin = (11.809735920000000,58.110175499999997)
Pixel Size = (0.000208320000000,-0.000208320000000)
Metadata:
  Band_1=Layer (Band 1:IS_20191119_m_40_20191125_s_41_cc_geo)
  Band_10=Layer (Band 1:IS_20190926_m_31_20191002_s_32_cc_geo)
  Band_11=Layer (Band 1:IS_20190920_m_30_20190926_s_31_cc_geo)
  Band_12=Layer (Band 1:IS_20190914_m_29_20190920_s_30_cc_geo)
  Band_13=Layer (Band 1:IS_20190908_m_28_20190914_s_29_cc_geo)
  Band_14=Layer (Band 1:IS_20190902_m_27_20190908_s_28_cc_geo)
  Band_15=Layer (Band 1:IS_20190827_m_26_20190902_s_27_cc_geo)
  Band_16=Layer (Band 1:IS_20190821_m_25_20190827_s_26_cc_geo)
  Band_17=Layer (Band 1:IS_20190815_m_24_20190821_s_25_cc_geo)
  Band_18=Layer (Band 1:IS_20190809_m_23_20190815_s_24_cc_geo)
  Band_19=Layer (Band 1:IS_20190803_m_22_20190809_s_23_cc_geo)
  Band_2=Layer (Band 1:IS_20191113_m_39_20191119_s_40_cc_geo)
  Band_20=Layer (Band 1:IS_20190728_m_21_20190803_s_22_cc_geo)
  Band_21=Layer (Band 1:IS_20190722_m_20_20190728_s_21_cc_geo)
  Band_22=Layer (Band 1:IS_20190716_m_19_20190722_s_20_cc_geo)
  Band_23=Layer (Band 1:IS_20190710_m_18_20190716_s_19_cc_geo)
  Band_24=Layer (Band 1:IS_20190704_m_17_20190710_s_18_cc_geo)
  Band_25=Layer (Band 1:IS_20190628_m_16_20190704_s_17_cc_geo)
  Band_26=Layer (Band 1:IS_20190622_m_15_20190628_s_16_cc_geo)
  Band_27=Layer (Band 1:IS_20190616_m_14_20190622_s_15_cc_geo)
  Band_28=Layer (Band 1:IS_20190610_m_13_20190616_s_14_cc_geo)
  Band_29=Layer (Band 1:IS_20190604_m_12_20190610_s_13_cc_geo)
  Band_3=Layer (Band 1:IS_20191107_m_38_20191113_s_39_cc_geo)
  Band_30=Layer (Band 1:IS_20190529_m_11_20190604_s_12_cc_geo)
  Band_31=Layer (Band 1:IS_20190523_m_10_20190529_s_11_cc_geo)
  Band_32=Layer (Band 1:IS_20190517_m_9_20190523_s_10_cc_geo)
  Band_33=Layer (Band 1:IS_20190511_m_8_20190517_s_9_cc_geo)
  Band_34=Layer (Band 1:IS_20190505_m_7_20190511_s_8_cc_geo)
  Band_35=Layer (Band 1:IS_20190429_m_6_20190505_s_7_cc_geo)
  Band_36=Layer (Band 1:IS_20190423_m_5_20190429_s_6_cc_geo)
  Band_37=Layer (Band 1:IS_20190417_m_4_20190423_s_5_cc_geo)
  Band_38=Layer (Band 1:IS_20190411_m_3_20190417_s_4_cc_geo)
  Band_39=Layer (Band 1:IS_20190405_m_2_20190411_s_3_cc_geo)
  Band_4=Layer (Band 1:IS_20191101_m_37_20191107_s_38_cc_geo)
  Band_40=Layer (Band 1:IS_20190330_m_1_20190405_s_2_cc_geo)
  Band_41=Layer (Band 1:IS_20190324_m_0_20190330_s_1_cc_geo)
  Band_5=Layer (Band 1:IS_20191026_m_36_20191101_s_37_cc_geo)
  Band_6=Layer (Band 1:IS_20191020_m_35_20191026_s_36_cc_geo)
  Band_7=Layer (Band 1:IS_20191014_m_34_20191020_s_35_cc_geo)
  Band_8=Layer (Band 1:IS_20191008_m_33_20191014_s_34_cc_geo)
  Band_9=Layer (Band 1:IS_20191002_m_32_20191008_s_33_cc_geo)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  11.8097359,  58.1101755) ( 11d48'35.05"E, 58d 6'36.63"N)
Lower Left  (  11.8097359,  57.2298152) ( 11d48'35.05"E, 57d13'47.33"N)
Upper Right (  13.7302380,  58.1101755) ( 13d43'48.86"E, 58d 6'36.63"N)
Lower Right (  13.7302380,  57.2298152) ( 13d43'48.86"E, 57d13'47.33"N)
Center      (  12.7699870,  57.6699953) ( 12d46'11.95"E, 57d40'11.98"N)
Band 1 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191119_m_40_20191125_s_41_cc_geo)
Band 2 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191113_m_39_20191119_s_40_cc_geo)
Band 3 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191107_m_38_20191113_s_39_cc_geo)
Band 4 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191101_m_37_20191107_s_38_cc_geo)
Band 5 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191026_m_36_20191101_s_37_cc_geo)
Band 6 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191020_m_35_20191026_s_36_cc_geo)
Band 7 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191014_m_34_20191020_s_35_cc_geo)
Band 8 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191008_m_33_20191014_s_34_cc_geo)
Band 9 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20191002_m_32_20191008_s_33_cc_geo)
Band 10 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190926_m_31_20191002_s_32_cc_geo)
Band 11 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190920_m_30_20190926_s_31_cc_geo)
Band 12 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190914_m_29_20190920_s_30_cc_geo)
Band 13 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190908_m_28_20190914_s_29_cc_geo)
Band 14 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190902_m_27_20190908_s_28_cc_geo)
Band 15 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190827_m_26_20190902_s_27_cc_geo)
Band 16 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190821_m_25_20190827_s_26_cc_geo)
Band 17 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190815_m_24_20190821_s_25_cc_geo)
Band 18 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190809_m_23_20190815_s_24_cc_geo)
Band 19 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190803_m_22_20190809_s_23_cc_geo)
Band 20 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190728_m_21_20190803_s_22_cc_geo)
Band 21 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190722_m_20_20190728_s_21_cc_geo)
Band 22 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190716_m_19_20190722_s_20_cc_geo)
Band 23 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190710_m_18_20190716_s_19_cc_geo)
Band 24 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190704_m_17_20190710_s_18_cc_geo)
Band 25 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190628_m_16_20190704_s_17_cc_geo)
Band 26 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190622_m_15_20190628_s_16_cc_geo)
Band 27 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190616_m_14_20190622_s_15_cc_geo)
Band 28 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190610_m_13_20190616_s_14_cc_geo)
Band 29 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190604_m_12_20190610_s_13_cc_geo)
Band 30 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190529_m_11_20190604_s_12_cc_geo)
Band 31 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190523_m_10_20190529_s_11_cc_geo)
Band 32 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190517_m_9_20190523_s_10_cc_geo)
Band 33 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190511_m_8_20190517_s_9_cc_geo)
Band 34 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190505_m_7_20190511_s_8_cc_geo)
Band 35 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190429_m_6_20190505_s_7_cc_geo)
Band 36 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190423_m_5_20190429_s_6_cc_geo)
Band 37 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190417_m_4_20190423_s_5_cc_geo)
Band 38 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190411_m_3_20190417_s_4_cc_geo)
Band 39 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190405_m_2_20190411_s_3_cc_geo)
Band 40 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190330_m_1_20190405_s_2_cc_geo)
Band 41 Block=9219x1 Type=Float32, ColorInterp=Undefined
  Description = Layer (Band 1:IS_20190324_m_0_20190330_s_1_cc_geo)

The code below stores the number of bands in a variable called nBands. I use this command thoughout my code several times.

[8]:
%%bash
nBands=$(gdalinfo Project/ccAll | grep Band_ | wc -l)
echo $nBands
41

In this snippet of code I create a mask file for each band (image) masks all the pixels with low quality ( < 0.25). Each line is commented with its functionality.

[ ]:
%%bash
nBands=$(gdalinfo Project/ccAll | grep Band_ | wc -l)
for band in `seq 1 $nBands`; do   # a for loop through all bands (images)
    echo "Band $band:"            # show the band number in progress
    gdal_translate -of  XYZ -b $band Project/ccAll Project/"B$band".txt  #convert the quality image to a text file (B*.txt)
    awk '{if ($3<0.25) {print $1,$2,0 } else {print $1,$2,1 }}' Project/"B$band".txt > Project/"mask$band".txt # The pixel values with quality<0.25 are set to 0 and otherwise are set to 1 in file B*.txt and write the new values to a new file (mask*.txt)
    gdal_translate -co COMPRESS=DEFLATE -co ZLEVEL=9  -ot Byte Project/"mask$band".txt Project/"mask$band".tif # Convert mask*.txt files to mask*.tif files
    rm -f Project/"B$band".txt # remove the created text files in previous steps
    rm -f Project/"mask$band".txt # remove the created text files in previous steps
done

In order to get a single mask file to ensure that a pixel has high quality in all images, we need to multiply all the mask images (with 0 and 1 values). However, for multiplication with gdal_calc.py we need to create a letter for each mask file, but we have only 26 letters because gdal version 2 only supports capital letters. Thus, we first multiply the first 26 mask*.tif images, then we multiply the rest, then we multiply the last two files and finally we multiply the maskAll.tif by all phase change maps.

[ ]:
%%bash
for i in {A..Z}; do
echo $i >> letters.txt # type letters A to Z in 26 lins of letters.txt
done
[2]:
!cat Project/letters.txt
A
B
C
D
E
F
G
H
I
J
K
L
M
N
O
P
Q
R
S
T
U
V
W
X
Y
Z
[ ]:
# Assign each letter with a dash to a mask*.tif file
%%bash
nBands=$(gdalinfo ccAll | grep Band_ | wc -l)
for count in `seq 1 26`; do
awk -v i=$count 'NR==i{ print "-"$1" mask"i".tif" }' letters.txt
done > list.txt
[3]:
!cat Project/list.txt
-A mask1.tif
-B mask2.tif
-C mask3.tif
-D mask4.tif
-E mask5.tif
-F mask6.tif
-G mask7.tif
-H mask8.tif
-I mask9.tif
-J mask10.tif
-K mask11.tif
-L mask12.tif
-M mask13.tif
-N mask14.tif
-O mask15.tif
-P mask16.tif
-Q mask17.tif
-R mask18.tif
-S mask19.tif
-T mask20.tif
-U mask21.tif
-V mask22.tif
-W mask23.tif
-X mask24.tif
-Y mask25.tif
-Z mask26.tif
[ ]:
# change list.txt file to a single line file
!grep - list.txt | tr '\n' ' ' > ls.txt
!var=$(cat ls.txt)
[5]:
!echo $var
-A mask1.tif -B mask2.tif -C mask3.tif -D mask4.tif -E mask5.tif -F mask6.tif -G mask7.tif -H mask8.tif -I mask9.tif -J mask10.tif -K mask11.tif -L mask12.tif -M mask13.tif -N mask14.tif -O mask15.tif -P mask16.tif -Q mask17.tif -R mask18.tif -S mask19.tif -T mask20.tif -U mask21.tif -V mask22.tif -W mask23.tif -X mask24.tif -Y mask25.tif -Z mask26.tif
[ ]:
# put an asterisk in front of each letter for gdal_calc.py sytanx
%%bash
for i in `seq 1 26`; do
awk -v i=$i '{ if (NR<26) print $1,"*"; else print $1; }' letters.txt > multiply.txt
done
[6]:
!cat Project/multiply.txt
A *
B *
C *
D *
E *
F *
G *
H *
I *
J *
K *
L *
M *
N *
O *
P *
Q *
R *
S *
T *
U *
V *
W *
X *
Y *
Z
[ ]:
# change multiply.txt file to a single line file
!grep . multiply.txt | tr '\n' ' ' > multiply2.txt
[7]:
!cat Project/multiply2.txt
A * B * C * D * E * F * G * H * I * J * K * L * M * N * O * P * Q * R * S * T * U * V * W * X * Y * Z
[ ]:
# multiply the first 26 mask*.tif images
!gdal_calc.py --type=Byte --overwrite  $var --co=COMPRESS=DEFLATE --co=ZLEVEL=9   --calc="( $(cat multiply2.txt) )" --outfile=mask1_26.tif
[ ]:
# the same procedure is done for the mask images mask27.tif to mask41.tif
%%bash
for count in `seq 27 41`; do
awk -v i=$count 'NR==i-26{ print "-"$1" mask"i".tif" }' letters.txt
done > list2.txt
############################################################
grep - list2.txt | tr '\n' ' ' > ls2.txt
var=$(cat ls2.txt)
echo $var
############################################################
for i in `seq 27 41`; do
awk -v i=$i '{ if (NR<42-27) print $1,"*"; else if (NR==42-27) print $1; else }' letters.txt > multiply3.txt
done
############################################################
grep . multiply3.txt | tr '\n' ' ' > multiply4.txt
cat multiply4.txt
############################################################
gdal_calc.py --type=Byte --overwrite  $var --co=COMPRESS=DEFLATE --co=ZLEVEL=9   --calc="( $(cat multiply4.txt) )" --outfile=mask27_41.tif
[ ]:
# Multiply the last two masks we got and create the final maskAll.tif file
!gdal_calc.py --type=Byte --overwrite  -A mask1_26.tif -B mask27_41.tif --co=COMPRESS=DEFLATE --co=ZLEVEL=9   --calc="( A * B )" --outfile=maskAll.tif
[ ]:
# Multiply each phase difference map (stored in bands of ifgAll file) by the maskAll.file to get the phase change of high-quality pixels
%%bash
for j in `seq 1 41`; do
gdal_calc.py --type=Float32 --overwrite  -A ifgAll --A_band=$j -B maskAll.tif --co=COMPRESS=DEFLATE --co=ZLEVEL=9   --calc="( A * B.astype(float) )" --outfile="ifg$j".tif
done

Crop the lake area

[18]:
# Crop the lake area for each phase change map with a shape file
%%bash
for j in `seq 1 41`; do
gdalwarp -cutline shape/tvarsjon60.shp -crop_to_cutline -dstnodata 0 "ifg$j".tif "ifg_crop$j".tif
done

The codes below create the time series of phase change and the accumulated phase change!

[19]:
# get the phase change of each pixel and store the time series of each pixel in each field of file timesS.txt
%%bash
nBands=$(gdalinfo ccAll | grep Band_ | wc -l)
for l in `seq 1 $nBands`; do
    gdal_translate -of  XYZ "ifg_crop$(($"42"-$l))".tif "G$l".txt # convert tif to text file, and also change the order of band beacause the oldest image in the original image is stored in band 41 and the newest image in band 1!
    awk 'BEGIN{ORS=" ";} {if ($3!=0) {print $3 }  else } END{print "\n"}' "G$l".txt >> timeS.txt # get only the pixels with non-zero values!
done
[ ]:
# Accumulate each field of timesS.txt
!awk '{ for (i=1; i<=NF; ++i) {sum[i]+=$i; $i=sum[i] }; print $0}' timeS.txt > timesAcc.txt
[10]:
# the number of high-quality pixels
!awk '{print NF}' Project/timesAcc.txt
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
482
0
[ ]:
!awk 'BEGIN{ORS=" ";} {if ($3!=0) {print "0" } }' G1.txt > zero.txt #create a line of zeros
!sed "1i\\$con" timesAcc.txt > tsAcc.txt #add zeros to the time series of the accumulated phase change (first line: we assume that the phase of the first date was zero)
[11]:
!cat Project/zero.txt
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[1]:
# number of pixels
!awk '{print NF}' Project/tsAcc.txt
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
483
[2]:
# number of dates
!awk '{print NR}' Project/tsAcc.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
[ ]:
# add the water level of the gauging station to the last column of the tsAcc.txt file
! awk '{getline to_add < "field.txt"; print $0,to_add}' tsAcc.txt > tsAcc.txt
[ ]:
# calculate the correlation between the last field (water level) and the accumulated phase change of all pixels (other fields) and store the values in corr.txt
%%bash
for count in `seq 1 483`; do
awk -v i=$count 'pass==1 {sx+=$i; sy+=$483; n+=1} pass==2 {mx=sx/(n); my=sy/(n); cov+=($i-mx)*($483-my); ssdx+=($i-mx)*($i-mx); ssdy+=($483-my)*($483-my);} END {print cov / ( sqrt(ssdx) * sqrt(ssdy) ) }' pass=1 tsAccT.txt pass=2 tsAccT.txt  >> corr.txt
done

1.11.2. Here, we see 10 pixels with highest correletation with field data

[6]:
!sort -k1 -n corr.txt | tail -10
0.665665
0.681602
0.697547
0.706916
0.714161
0.747608
0.77914
0.786415
0.806715
1
[4]:
!jupyter nbconvert myProject.ipynb --to html
[NbConvertApp] Converting notebook myProject.ipynb to html
[NbConvertApp] Writing 614090 bytes to myProject.html
[5]:
!jupyter nbconvert SaeidAminjafariProject.ipynb --to html
[NbConvertApp] Converting notebook myProject.ipynb to html
[NbConvertApp] Writing 614318 bytes to myProject.html
[ ]: