{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Phase Change Analysis\n", "\n", "-----------------------\n", "***Saeid Aminjafari***\n", "\n", "GeoComput. & ML\n", " \n", "2 June 2021" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## **Project description**\n", "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](https://search.asf.alaska.edu/). 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.\n", "\n", "### **Geocomputain datasets**\n", "For the course Geo-computation, I have two sets of data:\n", "1. The phase difference between successive dates (41 maps of phase change!).\n", "2. The quality images of the phase change maps.\n", "\n", "### **Data processing and computation procedure:**\n", "1. For each quality image, we mask the pixels with quality less than 0.25 (41 mask.tif images).\n", "2. Multiply all mask images to get the pixels that in all maps have high quality (maskAll.tif).\n", "3. Multiply the maskAll.tiff by each phase-change map to get the phase change of high-quality pixels.\n", "4. Crop the masked phase-change maps around the lake area with a shape file.\n", "5. Get the values of selected pixels and create a time series of phase change for each pixel (42 dates).\n", "6. Accumulate the time series of the phase change for each pixel.\n", "7. Calculate the correlation between the accumulated phase change and the water level of the gauging station.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**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.**" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Driver: ENVI/ENVI .hdr Labelled\n", "Files: Project/ccAll\n", " Project/ccAll.hdr\n", "Size is 9219, 4226\n", "Coordinate System is:\n", "GEOGCS[\"GCS_WGS_1984\",\n", " DATUM[\"WGS_1984\",\n", " SPHEROID[\"WGS_84\",6378137.0,298.257223563]],\n", " PRIMEM[\"Greenwich\",0.0],\n", " UNIT[\"degree\",0.0174532925199433]]\n", "Origin = (11.809735920000000,58.110175499999997)\n", "Pixel Size = (0.000208320000000,-0.000208320000000)\n", "Metadata:\n", " Band_1=Layer (Band 1:IS_20191119_m_40_20191125_s_41_cc_geo)\n", " Band_10=Layer (Band 1:IS_20190926_m_31_20191002_s_32_cc_geo)\n", " Band_11=Layer (Band 1:IS_20190920_m_30_20190926_s_31_cc_geo)\n", " Band_12=Layer (Band 1:IS_20190914_m_29_20190920_s_30_cc_geo)\n", " Band_13=Layer (Band 1:IS_20190908_m_28_20190914_s_29_cc_geo)\n", " Band_14=Layer (Band 1:IS_20190902_m_27_20190908_s_28_cc_geo)\n", " Band_15=Layer (Band 1:IS_20190827_m_26_20190902_s_27_cc_geo)\n", " Band_16=Layer (Band 1:IS_20190821_m_25_20190827_s_26_cc_geo)\n", " Band_17=Layer (Band 1:IS_20190815_m_24_20190821_s_25_cc_geo)\n", " Band_18=Layer (Band 1:IS_20190809_m_23_20190815_s_24_cc_geo)\n", " Band_19=Layer (Band 1:IS_20190803_m_22_20190809_s_23_cc_geo)\n", " Band_2=Layer (Band 1:IS_20191113_m_39_20191119_s_40_cc_geo)\n", " Band_20=Layer (Band 1:IS_20190728_m_21_20190803_s_22_cc_geo)\n", " Band_21=Layer (Band 1:IS_20190722_m_20_20190728_s_21_cc_geo)\n", " Band_22=Layer (Band 1:IS_20190716_m_19_20190722_s_20_cc_geo)\n", " Band_23=Layer (Band 1:IS_20190710_m_18_20190716_s_19_cc_geo)\n", " Band_24=Layer (Band 1:IS_20190704_m_17_20190710_s_18_cc_geo)\n", " Band_25=Layer (Band 1:IS_20190628_m_16_20190704_s_17_cc_geo)\n", " Band_26=Layer (Band 1:IS_20190622_m_15_20190628_s_16_cc_geo)\n", " Band_27=Layer (Band 1:IS_20190616_m_14_20190622_s_15_cc_geo)\n", " Band_28=Layer (Band 1:IS_20190610_m_13_20190616_s_14_cc_geo)\n", " Band_29=Layer (Band 1:IS_20190604_m_12_20190610_s_13_cc_geo)\n", " Band_3=Layer (Band 1:IS_20191107_m_38_20191113_s_39_cc_geo)\n", " Band_30=Layer (Band 1:IS_20190529_m_11_20190604_s_12_cc_geo)\n", " Band_31=Layer (Band 1:IS_20190523_m_10_20190529_s_11_cc_geo)\n", " Band_32=Layer (Band 1:IS_20190517_m_9_20190523_s_10_cc_geo)\n", " Band_33=Layer (Band 1:IS_20190511_m_8_20190517_s_9_cc_geo)\n", " Band_34=Layer (Band 1:IS_20190505_m_7_20190511_s_8_cc_geo)\n", " Band_35=Layer (Band 1:IS_20190429_m_6_20190505_s_7_cc_geo)\n", " Band_36=Layer (Band 1:IS_20190423_m_5_20190429_s_6_cc_geo)\n", " Band_37=Layer (Band 1:IS_20190417_m_4_20190423_s_5_cc_geo)\n", " Band_38=Layer (Band 1:IS_20190411_m_3_20190417_s_4_cc_geo)\n", " Band_39=Layer (Band 1:IS_20190405_m_2_20190411_s_3_cc_geo)\n", " Band_4=Layer (Band 1:IS_20191101_m_37_20191107_s_38_cc_geo)\n", " Band_40=Layer (Band 1:IS_20190330_m_1_20190405_s_2_cc_geo)\n", " Band_41=Layer (Band 1:IS_20190324_m_0_20190330_s_1_cc_geo)\n", " Band_5=Layer (Band 1:IS_20191026_m_36_20191101_s_37_cc_geo)\n", " Band_6=Layer (Band 1:IS_20191020_m_35_20191026_s_36_cc_geo)\n", " Band_7=Layer (Band 1:IS_20191014_m_34_20191020_s_35_cc_geo)\n", " Band_8=Layer (Band 1:IS_20191008_m_33_20191014_s_34_cc_geo)\n", " Band_9=Layer (Band 1:IS_20191002_m_32_20191008_s_33_cc_geo)\n", "Image Structure Metadata:\n", " INTERLEAVE=BAND\n", "Corner Coordinates:\n", "Upper Left ( 11.8097359, 58.1101755) ( 11d48'35.05\"E, 58d 6'36.63\"N)\n", "Lower Left ( 11.8097359, 57.2298152) ( 11d48'35.05\"E, 57d13'47.33\"N)\n", "Upper Right ( 13.7302380, 58.1101755) ( 13d43'48.86\"E, 58d 6'36.63\"N)\n", "Lower Right ( 13.7302380, 57.2298152) ( 13d43'48.86\"E, 57d13'47.33\"N)\n", "Center ( 12.7699870, 57.6699953) ( 12d46'11.95\"E, 57d40'11.98\"N)\n", "Band 1 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191119_m_40_20191125_s_41_cc_geo)\n", "Band 2 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191113_m_39_20191119_s_40_cc_geo)\n", "Band 3 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191107_m_38_20191113_s_39_cc_geo)\n", "Band 4 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191101_m_37_20191107_s_38_cc_geo)\n", "Band 5 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191026_m_36_20191101_s_37_cc_geo)\n", "Band 6 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191020_m_35_20191026_s_36_cc_geo)\n", "Band 7 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191014_m_34_20191020_s_35_cc_geo)\n", "Band 8 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191008_m_33_20191014_s_34_cc_geo)\n", "Band 9 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20191002_m_32_20191008_s_33_cc_geo)\n", "Band 10 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190926_m_31_20191002_s_32_cc_geo)\n", "Band 11 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190920_m_30_20190926_s_31_cc_geo)\n", "Band 12 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190914_m_29_20190920_s_30_cc_geo)\n", "Band 13 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190908_m_28_20190914_s_29_cc_geo)\n", "Band 14 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190902_m_27_20190908_s_28_cc_geo)\n", "Band 15 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190827_m_26_20190902_s_27_cc_geo)\n", "Band 16 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190821_m_25_20190827_s_26_cc_geo)\n", "Band 17 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190815_m_24_20190821_s_25_cc_geo)\n", "Band 18 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190809_m_23_20190815_s_24_cc_geo)\n", "Band 19 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190803_m_22_20190809_s_23_cc_geo)\n", "Band 20 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190728_m_21_20190803_s_22_cc_geo)\n", "Band 21 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190722_m_20_20190728_s_21_cc_geo)\n", "Band 22 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190716_m_19_20190722_s_20_cc_geo)\n", "Band 23 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190710_m_18_20190716_s_19_cc_geo)\n", "Band 24 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190704_m_17_20190710_s_18_cc_geo)\n", "Band 25 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190628_m_16_20190704_s_17_cc_geo)\n", "Band 26 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190622_m_15_20190628_s_16_cc_geo)\n", "Band 27 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190616_m_14_20190622_s_15_cc_geo)\n", "Band 28 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190610_m_13_20190616_s_14_cc_geo)\n", "Band 29 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190604_m_12_20190610_s_13_cc_geo)\n", "Band 30 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190529_m_11_20190604_s_12_cc_geo)\n", "Band 31 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190523_m_10_20190529_s_11_cc_geo)\n", "Band 32 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190517_m_9_20190523_s_10_cc_geo)\n", "Band 33 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190511_m_8_20190517_s_9_cc_geo)\n", "Band 34 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190505_m_7_20190511_s_8_cc_geo)\n", "Band 35 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190429_m_6_20190505_s_7_cc_geo)\n", "Band 36 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190423_m_5_20190429_s_6_cc_geo)\n", "Band 37 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190417_m_4_20190423_s_5_cc_geo)\n", "Band 38 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190411_m_3_20190417_s_4_cc_geo)\n", "Band 39 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190405_m_2_20190411_s_3_cc_geo)\n", "Band 40 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190330_m_1_20190405_s_2_cc_geo)\n", "Band 41 Block=9219x1 Type=Float32, ColorInterp=Undefined\n", " Description = Layer (Band 1:IS_20190324_m_0_20190330_s_1_cc_geo)\n" ] } ], "source": [ "!gdalinfo Project/ccAll" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**The code below stores the number of bands in a variable called *nBands*. I use this command thoughout my code several times.**" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "41\n" ] } ], "source": [ "%%bash\n", "nBands=$(gdalinfo Project/ccAll | grep Band_ | wc -l)\n", "echo $nBands" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**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.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%bash\n", "nBands=$(gdalinfo Project/ccAll | grep Band_ | wc -l)\n", "for band in `seq 1 $nBands`; do # a for loop through all bands (images)\n", " echo \"Band $band:\" # show the band number in progress\n", " gdal_translate -of XYZ -b $band Project/ccAll Project/\"B$band\".txt #convert the quality image to a text file (B*.txt)\n", " 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)\n", " 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\n", " rm -f Project/\"B$band\".txt # remove the created text files in previous steps\n", " rm -f Project/\"mask$band\".txt # remove the created text files in previous steps\n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**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.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%bash\n", "for i in {A..Z}; do\n", "echo $i >> letters.txt # type letters A to Z in 26 lins of letters.txt\n", "done" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A\n", "B\n", "C\n", "D\n", "E\n", "F\n", "G\n", "H\n", "I\n", "J\n", "K\n", "L\n", "M\n", "N\n", "O\n", "P\n", "Q\n", "R\n", "S\n", "T\n", "U\n", "V\n", "W\n", "X\n", "Y\n", "Z\n" ] } ], "source": [ "!cat Project/letters.txt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign each letter with a dash to a mask*.tif file\n", "%%bash\n", "nBands=$(gdalinfo ccAll | grep Band_ | wc -l)\n", "for count in `seq 1 26`; do\n", "awk -v i=$count 'NR==i{ print \"-\"$1\" mask\"i\".tif\" }' letters.txt \n", "done > list.txt" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-A mask1.tif\n", "-B mask2.tif\n", "-C mask3.tif\n", "-D mask4.tif\n", "-E mask5.tif\n", "-F mask6.tif\n", "-G mask7.tif\n", "-H mask8.tif\n", "-I mask9.tif\n", "-J mask10.tif\n", "-K mask11.tif\n", "-L mask12.tif\n", "-M mask13.tif\n", "-N mask14.tif\n", "-O mask15.tif\n", "-P mask16.tif\n", "-Q mask17.tif\n", "-R mask18.tif\n", "-S mask19.tif\n", "-T mask20.tif\n", "-U mask21.tif\n", "-V mask22.tif\n", "-W mask23.tif\n", "-X mask24.tif\n", "-Y mask25.tif\n", "-Z mask26.tif\n" ] } ], "source": [ "!cat Project/list.txt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# change list.txt file to a single line file\n", "!grep - list.txt | tr '\\n' ' ' > ls.txt \n", "!var=$(cat ls.txt)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-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 " ] } ], "source": [ "!echo $var" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# put an asterisk in front of each letter for gdal_calc.py sytanx\n", "%%bash\n", "for i in `seq 1 26`; do\n", "awk -v i=$i '{ if (NR<26) print $1,\"*\"; else print $1; }' letters.txt > multiply.txt \n", "done" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A *\n", "B *\n", "C *\n", "D *\n", "E *\n", "F *\n", "G *\n", "H *\n", "I *\n", "J *\n", "K *\n", "L *\n", "M *\n", "N *\n", "O *\n", "P *\n", "Q *\n", "R *\n", "S *\n", "T *\n", "U *\n", "V *\n", "W *\n", "X *\n", "Y *\n", "Z\n" ] } ], "source": [ "!cat Project/multiply.txt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# change multiply.txt file to a single line file\n", "!grep . multiply.txt | tr '\\n' ' ' > multiply2.txt " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "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 \n" ] } ], "source": [ "!cat Project/multiply2.txt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# multiply the first 26 mask*.tif images\n", "!gdal_calc.py --type=Byte --overwrite $var --co=COMPRESS=DEFLATE --co=ZLEVEL=9 --calc=\"( $(cat multiply2.txt) )\" --outfile=mask1_26.tif" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# the same procedure is done for the mask images mask27.tif to mask41.tif\n", "%%bash\n", "for count in `seq 27 41`; do\n", "awk -v i=$count 'NR==i-26{ print \"-\"$1\" mask\"i\".tif\" }' letters.txt\n", "done > list2.txt\n", "############################################################\n", "grep - list2.txt | tr '\\n' ' ' > ls2.txt\n", "var=$(cat ls2.txt)\n", "echo $var\n", "############################################################\n", "for i in `seq 27 41`; do\n", "awk -v i=$i '{ if (NR<42-27) print $1,\"*\"; else if (NR==42-27) print $1; else }' letters.txt > multiply3.txt\n", "done\n", "############################################################\n", "grep . multiply3.txt | tr '\\n' ' ' > multiply4.txt\n", "cat multiply4.txt\n", "############################################################\n", "gdal_calc.py --type=Byte --overwrite $var --co=COMPRESS=DEFLATE --co=ZLEVEL=9 --calc=\"( $(cat multiply4.txt) )\" --outfile=mask27_41.tif" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Multiply the last two masks we got and create the final maskAll.tif file\n", "!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" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Multiply each phase difference map (stored in bands of ifgAll file) by the maskAll.file to get the phase change of high-quality pixels\n", "%%bash\n", "for j in `seq 1 41`; do\n", "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\n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Crop the lake area**" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Crop the lake area for each phase change map with a shape file\n", "%%bash\n", "for j in `seq 1 41`; do\n", "gdalwarp -cutline shape/tvarsjon60.shp -crop_to_cutline -dstnodata 0 \"ifg$j\".tif \"ifg_crop$j\".tif\n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**The codes below create the time series of phase change and the accumulated phase change!**" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "tags": [] }, "outputs": [], "source": [ "# get the phase change of each pixel and store the time series of each pixel in each field of file timesS.txt\n", "%%bash\n", "nBands=$(gdalinfo ccAll | grep Band_ | wc -l)\n", "for l in `seq 1 $nBands`; do\n", " 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!\n", " 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!\n", "done" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Accumulate each field of timesS.txt\n", "!awk '{ for (i=1; i<=NF; ++i) {sum[i]+=$i; $i=sum[i] }; print $0}' timeS.txt > timesAcc.txt" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "482\n", "0\n" ] } ], "source": [ "# the number of high-quality pixels\n", "!awk '{print NF}' Project/timesAcc.txt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!awk 'BEGIN{ORS=\" \";} {if ($3!=0) {print \"0\" } }' G1.txt > zero.txt #create a line of zeros\n", "!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)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " ] } ], "source": [ "!cat Project/zero.txt" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n", "483\n" ] } ], "source": [ "# number of pixels\n", "!awk '{print NF}' Project/tsAcc.txt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n", "2\n", "3\n", "4\n", "5\n", "6\n", "7\n", "8\n", "9\n", "10\n", "11\n", "12\n", "13\n", "14\n", "15\n", "16\n", "17\n", "18\n", "19\n", "20\n", "21\n", "22\n", "23\n", "24\n", "25\n", "26\n", "27\n", "28\n", "29\n", "30\n", "31\n", "32\n", "33\n", "34\n", "35\n", "36\n", "37\n", "38\n", "39\n", "40\n", "41\n", "42\n" ] } ], "source": [ "# number of dates\n", "!awk '{print NR}' Project/tsAcc.txt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# add the water level of the gauging station to the last column of the tsAcc.txt file\n", "! awk '{getline to_add < \"field.txt\"; print $0,to_add}' tsAcc.txt > tsAcc.txt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# 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\n", "%%bash\n", "for count in `seq 1 483`; do\n", "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\n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Here, we see 10 pixels with highest correletation with field data" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.665665\n", "0.681602\n", "0.697547\n", "0.706916\n", "0.714161\n", "0.747608\n", "0.77914\n", "0.786415\n", "0.806715\n", "1\n" ] } ], "source": [ "!sort -k1 -n corr.txt | tail -10" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[NbConvertApp] Converting notebook myProject.ipynb to html\n", "[NbConvertApp] Writing 614090 bytes to myProject.html\n" ] } ], "source": [ "!jupyter nbconvert myProject.ipynb --to html" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[NbConvertApp] Converting notebook myProject.ipynb to html\n", "[NbConvertApp] Writing 614318 bytes to myProject.html\n" ] } ], "source": [ "!jupyter nbconvert SaeidAminjafariProject.ipynb --to html" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 4 }