Getting acquainted with the Google Earth engine happened in March 2020, when the coronavirus pandemic hit Indonesia and various sectors were affected. Almost all office activities were “laid off”, while food did not stop. Then how do we estimate the harvest area and harvest report if we don’t see the field?
After browsing the internet and having additional conversations with friends who are competent in their fields, it was decided to try the Google Earth engine application. The Google Earth engine is actually not much different from other mapping applications such as ArcGis and others. However, the specialty of Google Earth is:
Open source
No need to think about licenses when using Google Earth Engine. This application is very suitable for people from beginner to advanced levels. in fact, the image output from the google earth engine has been used by various international journal articles. Simply register to the google earth engine https://earthengine.google.com/ with the condition that you already have a gmail account because this is one group with other gmail web applications. Later there will be a few questions about what this application will be used for and what your job is, and others. just fill it in, and then there will be an email that comes in as activation and there is a link to https://code.earthengine.google.com/ to start using it.
Complete dataset of various satellites
I personally have only ever used arcgis, but as a layman, I am still confused if I want to use sentinel satellites, landsat, Modis, etc., whereas in the Google Earth engine, all we have to do is write the coding, and everything is provided by Google, even last week’s or yesterday’s satellite imagery if the satellite updates daily.
Custom
I say very custom because this application uses coding to get the desired image. I personally think that if you use ArcGIS, you have to look for books or literature if you want to clip maps or calculate areas, legends, etc. because the menu has been determined. But in Google Earth, everything is up to us according to our purpose for using it. Indeed, this can be a drawback because a layman usually does not know what to type. However, it turns out that many guides and tutorials have been scattered across blogs and YouTube. I myself code from journals that usually provide source links. So no need to worry; there are already many Google Earth users.
In this blog, I will write down some of the experiences I have had using Google Earth. In accordance with the background above, I got acquainted with Google Earth to calculate how much the growth phase area is so that we know how much the estimated harvest area is.
Actually, there is an application that has been prepared by the ministry of agriculture, namely SC2, or sentinel 2, with an area in each growth phase at the sub-district level. My impression is that because it only displays the land area, it is difficult to verify in the field. Of course, it is difficult to convince field instructors with direct data on a table of so many hectares without being able to show where the rice fields are located. Using Google Earth, in addition to knowing the area of each phase of rice growth, we will also get maps at the sub-district level, even villages if we have the source of village administration and the source of the raw land map of rice fields.
But before that, I will make a disclaimer that what I write here is the experience of a layman who is just learning for the first time with mapping and the Google Earth engine.
Introduction to Google Earth Engine Components
In this article, I will write about the parts of Google Earth in my own language. If you have successfully created an account and entered the Google Earth engine, it is divided into 3 columns. The rightmost part I call the library (no. 1) because it contains stored scripts, Google documents or tutorials, and files that we upload to the Google Earth engine with a capacity of 10 GB. The files we save can be maps in the form of SHP or map images from satellites.
Then number 2 is where the code will be written to get maps or other information from the Google Earth engine. In my experience, I don’t write code manually but look for code from other sources by copying and pasting, commonly known as Observe, Imitate, and Modify.
Column number 3 is console, inspector, and task. The console will contain tables, graphs, or images that appear according to the requests we write in the coding section. Task is a Google Earth engine processing for uploading and exporting images. While we can use the inspector to see the value at the point we select in the map area, we cannot see the values listed at that point. For example, the ndvi value and the ndwi value depend on what processing we write in the code.
Adding geometry assets into Google Earth Engine
Before going further in creating scripts in Google Earth Engine, we should understand how to add geometry assets in Google Earth Engine so that when we play the script, we are already focused on the AOI (area of interest).
On the left panel of the Google Earth Engine display, there are script, Docs, and Asset tabs. Asset itself is a file that we can save in the form of images, SHP in table form, and also a collection of image collections from the images we save. When this article is written, each account will get a free capacity of 10 GB.
Actually, without importing assets, we can learn the Google Earth engine by using points, rectangles, or custom shapes that we click into the Google Earth engine map, as shown below. Click the geometry mark on the top right of the map, then select the area using the mouse or cursor. Then in the script will appear import (1 entry) var geometry, which indicates we have entered the AOI variable into the script. Furthermore, this geometry is a variable that can be used to cut images provided by collections from sentinal satellites, landsats, or others.
However, it would be nice if we had a direct AOI that focused on our specific goals using the Google Earth Engine. Specifically, making costum geometry in the form of shapes that are the output of mapping applications such as Arcgis, Map Info, and others
The first thing that needs to be done is to prepare the material to be uploaded. I used the Arcgis application to cut the shape of one district’s rice fields into one sub-district’s rice fields with the “clip” arch tool. Furthermore, the data in the table in the layer Sawah one sub-district is exported in the form of SHP. The results of the files that will be seen are seven files, which we then archive into one in the form of a zip.
The picture above is a map of rice fields in Sidamanik District that I will upload on Google Earth Engine for the process of working on the rice plant growth phase map. Rice field maps and sub-district administration can be downloaded at RBI (Rupa Bumi Indonesia). It’s free, we are only asked to create an account.
After the geometry data is complete, open Google Earth Engine and click the asset tab. Then click new – shape files (.shp. shx, … zip).
Click Select and navigate to the result file that we compressed earlier. Click Upload, and then the upload progress will appear on the right side of the task tab. Make sure the process is complete before you refresh the asset tab to see the geometry import results.
After a successful import, we can click, and the results are as shown below:
To use it or import it into the script in the next Google Earth engine process, we just click import at the bottom.
This image is an example of the results of a script that is clipped with AOI with the geometry of rice fields in Sidamanik sub-district, Simalungun Regency, North Sumatra Province.
One of the best ways to take clear images from satellites
After successfully importing geometry, or AOI, into the Google Earth engine, this article will discuss how to take pictures from satellites with the geometry restrictions that we have.
We open the Google Earth engine, import geometry or an AOI (area of interest) that we have, and then in the code editor section (middle section), we type the following code:
Var geometry = table
var dataset = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2020-01-01', '2020-01-30')
.filterBounds(geometry)
print (dataset)
In the code above, we describe two variables, namely geometry and dataset. The geometry variable is the same as the table that contains the AOI shp file that we have imported. Actually, this variable does not need to be listed, actually it does not matter, because the table itself has been described at the beginning when we import it with the name ‘var table’.
dataset contains a collection of images owned by sentinel 2 or Copernicus/S2 by providing a time limit (filter date function from January 1 to January 30). This time limit is important, because if it is not limited google earth will provide all sentinel image collections and it will cause overload.
Then the last but not least thing is to limit the area. Or the area that we will want to see is the geometry area. Giving this limitation is represented by .flterbounds(geometry).
Then on the last line “print(dataset)” serves to list the image images in the console space or the area to the right of the google earth engine. This means that all sentinel two images in the time range and region that match the filters we include will have their file names appear on the right.
After that, click run on the toolbar above. Wait until there will be information in the console. The result turned out to be 12 sentinel collection images in January 2020. We can click on the row or the arrow, then we see the metadata information one by one. Each has 16 bands.
The circular area in the image above is the name of the image. We can then process data specifically using this image by writing or copying and pasting the title.
Display sentinel image 2
Write the script below to display the image we have selected:
var image = ee.Image('COPERNICUS/S2/20200102T034141_20200102T034813_T47NMD')
var image2= image.clip(geometry);
var visParams = {bands: ['B8', 'B4', 'B3'], max: 3048, gamma: 1};
var visParams_ndvi = {min: -0.2, max: 0.8, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
'3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};
// Calculate NDVI
var image_ndvi = image2.normalizedDifference(['B8','B4']);
// Map results
Map.centerObject(geometry,13)
Map.addLayer(image2,visParams,'warna merah')
Map.addLayer(image_ndvi,visParams_ndvi,'warna NDVI')
Inside ee.image is the image name we want based on our previous selection. Just copy and paste the image title. Do not forget that we clip the AOI region that we want in the Google Earth engine with the image 2 variable.
The new variables visparams and visparam_ndvi are our explanation of how we want to display the map, this time I only use band 8 or B8, B4, and B3. The band description is in the description of the previously selected image. Or read the profile of the sentinel satellite Top of Atmosphere at the following link: sentinel 2 band description.
Then in the selected bands, we calculate the NDVI value that will be displayed in the new layer later. The NDVI discussion will be explained separately.
Then in the map result section we start to create a layer, with the center object we focus on the AOI with a focus size equal to 13. This focus value can be enlarged or reduced as needed. Then we name the layer “red color” using the description described in the visparam variable. And the ndvi color for the ndvi layer
Click run
The result looks as follows:
If we remove the clip geometry, the original image is actually square. As follows:
Script:
var image = ee.Image('COPERNICUS/S2/20200102T034141_20200102T034813_T47NMD')
var visParams = {bands: ['B8', 'B4', 'B3'], max: 3048, gamma: 1};
var visParams_ndvi = {min: -0.2, max: 0.8, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
'3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};
// Calculate NDVI
var image_ndvi = image.normalizedDifference(['B8','B4']);
// Map results
Map.centerObject(image,9)
Map.addLayer(image,visParams,'warna merah')
Map.addLayer(image_ndvi,visParams_ndvi,'warna NDVI')
Result:
Red color image:
It can be seen that there are clouds that cover a fairly large area of land. This is understandable because Indonesia has a tropical climate with cloud conditions that can be said to exist every day. It is likely that the eastern region of Indonesia will be cleaner from clouds because of the hotter weather.
One by one, images from the 12 images at the beginning can be tried to see which number of clouds is the least. But there is a more precise way to find out which image has the least number of clouds, namely by using the sort function on the number of clouds.
The image with the least number of clouds in a certain time range
Use the script below:
var geometry = table
var dataset = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2020-01-01', '2020-01-30')
.filterBounds(geometry)
var sorted = dataset.sort('CLOUD_COVER')
var image = sorted.first()
var visParams = {bands: ['B8', 'B4', 'B3'], max: 3048, gamma: 1};
var visParams_ndvi = {min: -0.2, max: 0.8, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
'3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};
// Calculate NDVI
var image_ndvi = image.normalizedDifference(['B8','B4']);
// Map results
Map.centerObject(geometry,12)
Map.addLayer(image.clip(geometry),visParams,'warna merah')
Map.addLayer(image_ndvi.clip(geometry),visParams_ndvi,'warna NDVI')
What is different from the above script is
var sorted = dataset.sort(‘CLOUD_COVER’)
var image = sorted.first()
So we immediately create a short in the cloud cover field and select the least cloud cover without having to specify or select the name of the image. So from the 12 images, without knowing the title of the image that will be selected, we immediately command to select the image with the least cloud cover in January 2020.
The result is as follows:
It can be seen that even though there is a slight cloud filter, the clouds still seem to be visible. Especially coastal areas or lake areas, where clouds are usually present every day.
Another process to remove clouds is by combining several images of the period to get a clearer layer called masking. So the pixels that are covered by clouds are replaced with other images that are not covered. This will probably be explained in the next article.
Calculating Area Based on NDVI Value
Making maps with NDVI descriptions is important, but it feels incomplete if we don’t add information on the area of each NDVI value range.
Based on my experience, I got a pdf map of plant growth phases from a well-known institution that has a special task of issuing such maps. As a result, I had to re-render the map to see how much land was in the vegetative phase and so on. That was quite troublesome.
The Google Earth engine provides a direct pixel calculation feature contained in the resulting map. The pixel calculation can then be multiplied by the scale used so that the area results come out directly, either in square meters or hectares.
In this article, I will divide the NDVI value into several ranges. There is a relationship between the plant growth phase and the NDVI value, so it can be grouped based on a certain range. This article does not discuss the phases of certain plants in a certain range but rather serves as an exercise on how to divide the NDVI range values and calculate the area.
Let’s open the Google Earth engine and write the following script:
Imports (2 entries)
var table : Table users/
var table2 : Table users/
var dataset = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2020-03-01', '2020-04-09')
.filterBounds(table);
var sorted = dataset.sort('CLOUD_COVER')
var scene = sorted.first().clip(table)
// Visualization parameters
var visParams = {bands: ['B8', 'B4', 'B3'], max: 3048, gamma: 1};
var visParams_ndvi = {min: -0.2, max: 0.8, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
'3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};
// Calculate NDVI
var image_ndvi = scene.normalizedDifference(['B8','B4']);
// Map results
Map.centerObject(table,12)
Map.addLayer(scene,visParams,'Sentinel-2 False Color Infrared')
Map.addLayer(table2,{},'kecamatan')
Map.addLayer(image_ndvi,visParams_ndvi,'Sentinel-2 NDVI')
var mergedAllFunction = function(image) {
var ndvi = scene.normalizedDifference(['B8','B4']).rename('NDVI');
var thres1 = ndvi.gte(0).rename('thres1')
var thres2=ndvi.gt(0.1921).and(ndvi.lt(0.739)).rename('thres2')
var thres3=ndvi.gt(0.74).and(ndvi.lt(0.778)).rename('thres3')
var thres4=ndvi.gt(0.192).and(ndvi.lte(1)).rename('thres4')
return image.addBands(ndvi).addBands([thres1,thres2,thres3,thres4]);
}
var median = dataset.median();
var merged = mergedAllFunction(median);
var areas = merged
.select(['thres1', 'thres2', 'thres3', 'thres4'])
.multiply(ee.Image.pixelArea())
.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: table, // a geometry
scale: 10, // scale = 10 for sentinel-2 'red' band
maxPixels: 1e9
});
print(areas,'square meters')
The first part is to import the assets of the rice field table (table 1) and the sub-district boundary (table 2). Then call the sentinel dataset collection by replacing the start and end boundary dates of the sentinel collection by applying the most minimal cloud filter.
The NDVI range division is located at:
var thres1 = ndvi.gte(0).rename(‘thres1’) var thres2=ndvi.gt(0.1921).and(ndvi.lt(0.739)).rename(‘thres2’) var thres3=ndvi.gt(0.74).and(ndvi.lt(0.778)).rename(‘thres3’) var thres4=ndvi.gt(0.192).and(ndvi.lte(1)).rename(‘thres4’)
the variable named thres1 has a range ndvi between 0 – 0.1920
thres2 has a range of 0.1921 – 0.738
thres3 has a range of 0.74 – 0.778
thres4 has a range of 0.192 – 1
The result of the script above is the image below. the result of the calculation of the area per class is in the console area with units of square meters. In addition, the map of rice fields with sub-district boundaries is also in the map section as explained in the previous article.
Leave a Reply