GEE: extarct and plot data from sentinel 2A_B imagery

 

First, Let us pull out a toa image from the S2_TOA collection by using consecutive filter methods on the collection. Before that, we will be pulling out a sentinel image from Sioux Falls, South Dakota. The tile number is T14TPP.

Lets, look for the day of 06/02/2019 because I know for sure that we have an image on that day

I also have a local copy of that image to tinker within QGIS. So, we can verify the results locally.

First, define some date ranges for our collection filter,

var start_date='2019-06-01';

var end_date='2019-06-03';

Then, apply the filters to the collection , then sort in ascending order and then just select the very first one from the sorted list

var s2_im=ee.Image(S2_TOA

                  .filterDate(start_date,end_date)

                  .filterBounds(point)

                  .sort('CLOUD_COVER')

                  .first());

 

Pack the visulaiztion parameters in a dictionary object

var trueColor={bands:['B4','B3','B2'],min:0,max:6000};

 

Use the addlayer method of Map class to add the image in the map

Map.addLayer(s2_im,trueColor,'True Color')

 

 print(s2_im,'all bands of the image')

//Map.addLayer(s2_im)

Clicking on the console tab we can explore bands available for the s2_im variable.

The expansion of the Image on the console yields the following. It has all the bands from Sentinel. 

Now, we will filter only the VNIR bands (bands which are Red (B4), Green (B3), Blue(B2) and NIR(B8) )

To select only VNIR bands we need to crate a list of the band IDs . After creating the list we need to pass the bandID  list to the ‘select’ method of the Image ‘s2_im’ object.

 

Now, to select only VNIR bands, make a variable that corresponds to the VNIR bands

var vnir_list=['B2','B3','B4','B8'];

var s2_im_vnir=s2_im.select(vnir_list);

var wavelengths=[0.490, 0.560, 0.665, 0.842];//optional for legend only

Check the s2_im  by printing to see if only vnir is selected

print(s2_im_vnir,'vnir filterd')

The output should be like the following after expansion in the console tab.

So, we have succesfully selected the only VNIR bands. The following part will do the tricks for extraction and plotting of the data--

Extarcting and plotting the VNIR data from the point geometry:

We need to configure the visualization option parameter first. To do that, we make a dictionary/object with the options.

var options={

title: 'VNIR raw data extarction from S2 using GEE (point Geometry)',

hAxis:{title:'Bands'},

vAxis:{title:'Reflectance(Raw DN)'},

lineWidth:2,

pointSize:7

};

 

 

This will create a chart in the console window when we print the chart object.

var chart=ui.Chart.image.regions(

  s2_im_vnir,point,null,null,null,vnir_list

  ).setOptions(options)

 print(chart,'Extract from point geometry')

Extarcting and plotting the VNIR data from the polygon geometry:

Similar method is used for reducing a polygon geometry. Also, we need to configure the visualization option parameter again here.

var options={

title: 'VNIR raw data extarction from S2 using GEE (Polygon Geometry)',

hAxis:{title:'Bands'},

vAxis:{title:'Reflectance(Raw DN)'},

lineWidth:2,

pointSize:7

};

var chart=ui.Chart.image.regions(

  s2_im_vnir,poly,ee.Reducer.mean(),null,null,vnir_list

  ).setOptions(options)

 

print(chart,'Extract from rectangle geometry')

The results are shown below:

Validation part: 

So the pixel values received from the google earth engine is plotted on the above picture.

To do the validation we will be using QGIS and a hand downloaded scene (same scene)

.Before that we want to export the point and polygon as kml/shp file. The process of doing that is fairly simple.

This snippet of comments are for extracting the geometry to kml or shape file. It will be exported to the google drive and can be  accessed through the account

First, make a collection of points.

var features = ee.FeatureCollection([

  ee.Feature(poly, {name: 'Point_of_interest_S2'})

]);

Export the FeatureCollection to a KML/SHP file.

Export.table.toDrive({

  collection: features,

  description:'rect_angle',

  fileFormat: 'SHP'

});

We can see these export code in the Tasks tab. It will not be exported to the google drive until we run the code from the Tasks list.

Opening and verification in QGIS:

I have already exported the point and polygon in the google drive and loaded in the local computer. So, I am doing a screenshot of  the QGIS window. and the QGIS window output is the following.

If we zoom in we can see the point and the polygon.

Now to validate the point geometry we need to zoom in to the point and get to the pixel level and get the DN value of the desired layer. To do that we need the 'value tool' for qgis. it can be installed as plugin if not available. We will need a 'zonal statistics'  plugin also for polygon verification. 

Doing a zoom to the pixel level and checking the DN value of that pixel through the value tool reveals that the Band2 or Blue Band value is 1643. Now, go back to the chart generated from GEE for point geometry,w here I intentionally put the cursor to on the B2(blue point) where it revealed 1643 also. So we extracted the right pixel.

Now, to verify the mean DN value of the polygon we need to run zonal statistics with the polygon shapefile and B2. To do that click on the zonal statistics toolbox from the 'processing toolbox section' and select the polygon layer and Blue band and process it.

Now, if we check the attribute table of the shape file we will see the mean value of all pixels inside the polygon for Blue band. 

The value we got in return is 1729.59810. Going back to the GEE returned calculation, we see the value is 1729.635 . The mean value is also matching 

The codes all together can also be found here.

Today's Python Learning: Enumerate function in Python with example

Enumerate is a built-in function in python which adds a counter to an iterable (iterable are objects that can be used with loop e.g., lists, dictionaries, etc.). So, often, we set a separate counter variable in the loop; enumerate function can save us from that task. Let's look at the following example.



In the above example, we created a list of strings and passed it to the enumerate method. The enumerate method returns an enumerator object. For printing convenience, we converted it to a list and printed the tuples.

Then we did the tuple unpacking using a for loop to see the indexes.
The enumerate()  function also accepts the initial index setting, which is by default set to zero. So, the default start index is always zero.



Bonus  learning tips for Jupyter notebook:

Cells can be used with the shortcut in two ways -
1. Command Mode ( When the cell borders are blue)
2.Edit mode (When the cell border is green)
Clicking outside or inside the cells  can toggle between the edit modes

# If we want to insert a new cell below, we can use the shortcut button "B" in command mode
# Similarly double-pressing the "D" key twice will delete the current cell
# Shift+Enter(Return) will  execute the cell

Today's python learning: Cascaded Assignments, Simultaneous assignment, Multiline statement, Multiline strings, split() method and eval method

Today I learned:

Cascaded Assignments:

We can assign the same value to multiple variables at the same time. For example, if we write

x=y=z=10

there will be 3 separate variables each of them holding 10 in their content. We can see from there content. This variable type and content explorer can be viewed from the Spyder IDE, which is cool and intuitive.

Simultaneous assignment:

We can also do simultaneous assignments in Python which is not commonly available to other languages. To do that we have to use a comma between expressions and an equal (assignment) operator at the right side. The following variable swapping problem is a perfect example of a simultaneous assignment.

We can see from the output that it works perfectly.

Multiline Statement:

This can be constructed using either escaping a newline character or using the parentheses.

Multiline String:

The multiline string can be done using enclosing strings into either a set of triple-double quotation marks or single quotation marks.

Follow the following examples for clarification.

Split method:

This is a very important method. It returns a list of strings after breaking the given string by the specified separator It accepts two parameters. Separator and the Max split. The following is an example of the split method, where the default separator is space.


eval() function:

This is a powerful tool as it evaluates the strings into python expressions. So we need to be careful about this tool. It will be an amazingly handy tool for building a calculator app or some graphing tools. The example can be found below.

Python Practice problem 4: Check for palindrome

Palindromes are the words that do not get changed when you put it in reverse order. So, "Mom", "Madam" are examples of the palindrome. This is a very simple problem where we use reversing a list using negative list slicing.

Sunday Morning bird searching with my new Nikon D7500

    I woke up at 9 am and had some cereal. Then, I took out my Nikon D7500 (Recently I have sold my old Nikon D5300+18-55 kit lens in the Facebook marketplace and bought this one ) and attached the 70-300mm telephoto lens on it. The objective was to go out and search for some birds in this area. It is winter here and the surroundings are mostly covered with snow with no leaves in the trees. It is pretty hard to find any birds and animals outside except some squirrels here and there occasionally. However, today's temperature was not negative. It was averaging at 1-2 degrees celsius. So, I thought that it will be fun walking in Elmwood Park.

    This park is my recent favorite place because of its location near to the Big Sioux River. A lot of birds come and rest there in the afternoon. For some reason, the water is rarely frozen there. However, you won't find any bird at the river in the daytime. So I wanted to just walk by the rive and there is a nice road beside the river which in my opinion, one of the best places for my relaxation in Sioux falls. After parking my car, I walked to the riverside road and walked for some time. It was sunny whole day. You can walk inside the park for some time and then get to the riverside road.

    While walking through the woods, I heard chirpings of the bird but, interestingly, I was not able to see any bird in the woods. Later, I discovered that the birds were too small to skip if you don't look carefully. To add more complexity the trees a the park were  huge in height, so it posses a tough challenge to locate those tiny birds. I felt that my 70-300 is not enough for me (I am thinking of the 200-500mm for a long time and I might get one soon). After spending about one hour I got to manage some shots. After coming back to home I processed some of the images. Oh one thing, I started shooting RAW and it made a night and day difference in the photo I take. I am very surprised by the image quality of this awesome camera (next target is to get a full-frame; Nikon D780 or Z6 might be the best option for me ...LoL). However, I am attaching some pictures from the Elmwood Park.

Python Practice problem 3: Determining the overlaps of two lists.

Problem:

This is the practice problem 5 from https://www.practicepython.org/exercise/2014/03/05/05-list-overlap.html. The problem is -


1. We need to make two lists, the items of the lists might be unique or have common items within them.

2. We need to write a program that returns a list containing only the elements that are common between the lists while not considering the duplicates.

For example, we are given two lists- list_A and list_B as the following

list_A=[1,3,5,6,7,8,9,10,13,15,20]

list_B=[2,4,6,8,9,10,12,14]

Our program should be able to find duplicates.


Solving Part:

The problem is pretty straight forward, but our objective is to explore different ways and learn python through problem-solving.

While solving problems, one problem I often fetch is indexing error which is caused by putting variable() instead of variable[]. The problem arises because I came from a Matlab background where the function calling and the variable indexing uses the same "()" which is confusing actually.

Anyway, let's go to coding. The first quick draft is shown at the following which is not taking care of duplicate values. The first thing comes to mind is the use of Brute Force. which is -

1. Take any of the two lists and iterate through it
2. While iterating through the items, check for the existence in the second list using "in" keyword.

The code is shown below-



You can obviously see the problem here. Look at the program output. [6, 8, 9, 10, 10]. We have duplicates in the result. It is expected. The last piece of code I  want to add to this program is, well I know it's not the best way to do this, is using a set operation on the generated list to remove the duplicate items i.e., to keep the unique items only in the list.

So the final version is here.

Learnings:

Few things are important to note from here

"in" operator:
The "in" operator in python is for checking membership of value against a sequence. For example,
If we declare a list and then if we want to check for the existence of an item in that list we can use this keyword "in". I am writing these expressions for clarification

There are two operators which are called membership operators. 
a. in operator
b. not in operator
About 'not in' operator, It basically does the opposite of  "in" operator. i.e., the expression containing it evaluates to true if it does not find the variable in the sequence.

We need to learn about Identity operators too. There are two identity operators like the membership operator.

1. 'is' operator: It checks if the object on the left side of  'is' is pointing to the same object on the right side. Personally, at first, I was confused about this operator.  But, when I think of the variables as a just human-readable version of memory tags,  it becomes a little bit clear. The following example and the output is a good way to understand the behavior of the 'is' operator. Here we took three variables p,q, and r. Then we initialized p and q with the same content. In the third line of code, we make p and r equal. As expected, although the content of the p and q are the same, they are pointing to the different locations in the memory. The code snippets will clarify the rest-


2.' is not' operator: just the opposite of the 'is' operator. 

Python Practice Problem 2: Print the Fibonacci series up to a input point

Let us go to practice problem 14 ( https://www.practicepython.org/exercise/2014/04/30/13-fibonacci.html) and solve the problem of the Fibonacci series. If you want more information on this special type of numbers you can google the term or you can go to https://en.wikipedia.org/wiki/Fibonacci_number, they have some practical examples also.
I am putting my solution here .