[PYTHON] [Sequel] Artificial satellite image analysis by Google Earth Engine and Google Colab-Satellite image analysis starting for free (practice)-

Introduction

-Artificial satellite image analysis by Google Earth Engine and Google Colab-Satellite image analysis starting for free (Introduction)-, with Google Earth Engine (GEE) Introduced how to acquire satellite images using Google Colab ――In this article, we will introduce the analysis method of typical ** artificial satellite data (optical image, vegetation index, ground surface temperature) ** that can be used in GEE. --Various other satellite datasets are stored in GEE's Data Catalog, and the variables in the code below (satellite name and satellite name) If you change the band name), you can perform the same analysis.

Code to get satellite data

--Use the code below (see Previous article for a detailed explanation of the code) --"Satellite name", "Band name", "Data period", "Target area", and "Save folder name" are variables, and when the function is executed, the satellite image is saved in the specified folder (Google Drive).

Function definition

#Earth Engine Python API Import
import ee

#GEE authentication / initialization
ee.Authenticate()
ee.Initialize()

#GEE data load
def load_data(snippet, from_date, to_date, geometry, band):
    #Extract data according to parameter conditions
    dataset = ee.ImageCollection(snippet).filter(
    ee.Filter.date(from_date, to_date)).filter(
    ee.Filter.geometry(geometry)).select(band)
    #Convert to list type
    data_list = dataset.toList(dataset.size().getInfo())
    #Output the number of target data and data list
    return dataset.size().getInfo(), data_list

#Save satellite imagery to Google Drive
def save_on_gdrive(image, geometry, dir_name, file_name, scale):
    task = ee.batch.Export.image.toDrive(**{
        'image': image,#Satellite information to load
        'description': file_name,#File name to save
        'folder':dir_name,#Save destination folder name
        'scale': scale,#resolution
        'region': geometry.getInfo()['coordinates'],#Target area
        'crs': 'EPSG:4326'
    })
    # Run exporting
    task.start()
    print('Done.')

Variable (parameter) setting

##Parameter specification
#Designate satellite
snippet = 'NOAA/DMSP-OLS/NIGHTTIME_LIGHTS'
#Specify the band name
band = 'avg_vis'
#Specify the period
from_date='2010-01-01'
to_date='2012-12-31'
#Specify area(Specify latitude / longitude in Japan area)
geometry = ee.Geometry.Rectangle([128.60, 29.97, 148.43, 46.12])
#Folder name to save
dir_name = 'GEE_download'
#file name
file_name = 'file_name'
#resolution
scale = 1000

Execution of processing

##Execution of processing----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

##Save all (File name uses satellite ID)
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

Satellite data analysis ① ~ Optical image ~

Overview

--Color image expressed by combining visible light bands (RGB) --Image of a commonly used "camera" photo --Each sensor data of RGB can be acquired from GEE and visualized by RGB composition.

data set

--Satellite = "Sentinel 2" --Band = TCI_R, TCI_G, TCI_B --Resolution: 10m --Area: Tokyo (mainly the Imperial Palace)

code

--Change the parameters and try to get TCI_R (red band of visible light)

##Parameter specification
#Designate satellite
snippet = 'COPERNICUS/S2_SR'
#Specify the band name
band = 'TCI_R'

#Specify the period
from_date='2020-04-01'
to_date='2020-04-15'
#Specify area(Specify latitude / longitude in Japan area)
geometry = ee.Geometry.Rectangle([139.686, 35.655, 139.796, 35.719])

#Folder name to save
dir_name = 'GEE_Sentinel_Red'
#file name
file_name = 'file_name'
#resolution
scale = 10

##Execution of processing----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

##Save all (File name uses satellite ID)
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

-Next, let's visualize the downloaded satellite image.

#Package installation&import
!pip install rasterio
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import rasterio
import glob

#Data reading
with rasterio.open('/content/drive/My Drive/GEE_Sentinel_Red/COPERNICUS_S2_SR_20200402T012651_20200402T012651_T54SUE.tif') as src:
    arr = src.read()

#Visualization
plt.imshow(arr[0], cmap='Reds')

Optical image of Tokyo (around the Imperial Palace) (red band)

――You can see the area that looks like the Imperial Palace in the center, and Yoyogi Park and Shinjuku Gyoen on the west side. ――Next, I will try to get other green and blue bands as well.

#Green band data acquisition
##Parameter specification
#Specify the band name
band = 'TCI_G'
#Folder name to save
dir_name = 'GEE_Sentinel_Green'

##Execution of processing----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

##Save all (File name uses satellite ID)
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

#Blue band data acquisition
##Parameter specification
#Specify the band name
band = 'TCI_B'
#Folder name to save
dir_name = 'GEE_Sentinel_Blue'

##Execution of processing----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

##Save all (File name uses satellite ID)
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

――I was able to get the green band and blue band as well as the red band. --Similarly, read the data and visualize it.

#Data reading
# Red
with rasterio.open('/content/drive/My Drive/GEE_Sentinel_Red/COPERNICUS_S2_SR_20200402T012651_20200402T012651_T54SUE.tif') as src:
    red = src.read()
# Green
with rasterio.open('/content/drive/My Drive/GEE_Sentinel_Green/COPERNICUS_S2_SR_20200402T012651_20200402T012651_T54SUE.tif') as src:
    green = src.read()
# Blue
with rasterio.open('/content/drive/My Drive/GEE_Sentinel_Blue/COPERNICUS_S2_SR_20200402T012651_20200402T012651_T54SUE.tif') as src:
    blue = src.read()

#Visualization
plt.figure(figsize=(15, 10))
plt.subplot(131); plt.imshow(red[0], cmap='Reds'); plt.title('Red')
plt.subplot(132); plt.imshow(green[0], cmap='Greens'); plt.title('Green')
plt.subplot(133); plt.imshow(blue[0], cmap='Blues'); plt.title('Blue')

Optical image of Tokyo (around the Imperial Palace) (red band, green band, blue band)

--Finally, these three color images are RGB-combined and converted into a color image. --RGB composition is just combined with numpy dstack

#RGB composition
## np.Connect with dstack (Red, Green,Note the order of Blue)
rgb = np.dstack((red[0], green[0], blue[0]))

#Visualization of RGB composited image
plt.imshow(rgb); plt.title('RGB Image')

Optical image of Tokyo (around the Imperial Palace) (RGB composition)

――By RGB composition, it looks like a satellite image. ――This time, we visualized only at a specific area timing, but data can be acquired at any timing within the period when the artificial satellite is in operation. ――It seems that you can see various things by changing the seasons or comparing before and after the disaster (optical images are easily affected by clouds, so the images may be unclear in the rainy season or on bad weather. There is)

Satellite data analysis ② ~ Vegetation index ~

Overview

――Indicator showing the distribution and activity of plants --Measured using a sensor with a wavelength that easily reacts to plants --Vegetation indexes have been proposed, but the NDVI (Normalized Difference Vegetation Index) is famous. --Used for desertification, forest fire monitoring, urban use classification, etc.

data set

--Satellite = "Terra (MODIS)" --Band = NDVI -* Normally, NDVI requires calculation by combining multiple bands, but GEE stores the calculated data (convenient ... !!) --Resolution: 250m --Area: Tokyo (mainly the Imperial Palace)

code

--Try to get NDVI (vegetation index) by changing various parameters

##Parameter specification
#Designate satellite
snippet = 'MODIS/006/MOD13Q1'
#Specify the band name
band = 'NDVI'

#Specify the period
from_date='2005-01-01'
to_date='2005-12-31'
#Specify area(Specify latitude / longitude in Japan area)
geometry = ee.Geometry.Rectangle([139.686, 35.655, 139.796, 35.719])

#Folder name to save
dir_name = 'GEE_NDVI'
#resolution
scale = 250

##Execution of processing----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

##Save all (File name uses satellite ID)
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

```python
## Parameter specification
# Designate satellite
snippet = 'MODIS/006/MOD13Q1'
# Specify the band name
band = 'NDVI'
# Specify the period
from_date='2005-01-01'
to_date='2005-12-31'
# Specify area (specify latitude / longitude for Japan area)
geometry = ee.Geometry.Rectangle([139.686, 35.655, 139.796, 35.719])

# Folder name to save
dir_name = 'GEE_NDVI'
# resolution
scale = 250

## Execution of processing ------------------------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## Save all (File name uses satellite ID)
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

-I will read the acquired data in the same way and visualize it

# Data reading
with rasterio.open('/content/drive/My Drive/GEE_NDVI/MODIS_006_MOD13Q1_2005_01_01.tif') as src:
    arr = src.read()

# Visualization
plt.figure(figsize=(20, 8))
plt.subplot(121); plt.imshow(rgb); plt.title('RGB Image'); plt.title('Optical Image')
plt.subplot(122); plt.imshow(arr[0], cmap='YlGn'); plt.title('NDVI')

Optical RGB image and vegetation index of Tokyo (around the Imperial Palace)

-The Imperial Palace and parks in Tokyo have high vegetation indexes. -Next, let's visualize all the data of the acquired period (1 year in 2005). -Vegetation is generally early summer~It activates in the fall and declines in the winter, so let's check how it looks.

# Visualization in chronological order
files = glob.glob('/content/drive/My Drive/GEE_NDVI/*tif')
files.sort()

plt.figure(figsize=(15, 10))
j=0
for i in range(len(files)):
 # Get images one scene at a time and visualize
  with rasterio.open(files[i]) as src:
      arr = src.read()
 j + = 1 # Shift and place the plot position of the image
  plt.subplot(5,6,j)
  plt.imshow(arr[0], cmap='YlGn', vmin=-2000, vmax=10000)
 plt.title (files [i] [-14:]) # Get date from filename
  plt.tight_layout()

Changes in the vegetation index in Tokyo (around the Imperial Palace) (January)~December)

• 5~You can see that the vegetation index is activating around September. -It's a little difficult to compare, so let's calculate the total value of the vegetation index and see the time series trend.

# Get the total value of NDVI in the area
sum_NDVI = []
label_NDVI = []
for i in range(len(files)):
 # Get images one scene at a time and visualize
  with rasterio.open(files[i]) as src:
      arr = src.read()
  sum_NDVI.append(arr.sum())
  label_NDVI.append(files[i][-14:-4])

# Visualization
fig,ax = plt.subplots(figsize=(15,3))
plt.plot(sum_NDVI, marker='o')
ax.set_xticks(np.arange(0,23))
ax.set_xticklabels(label_NDVI, rotation=90)
plt.show()

Changes in vegetation index (January)~December)

-You can see how the vegetation index peaks over the summer -It is possible that June is a major phenomenon because the satellite cannot acquire data correctly due to the influence of the rainy season (clouds). -Optical satellites have the weakness of being easily affected by clouds, but have the advantage of being able to acquire a variety of information (on the other hand, there are satellites of the type called SAR that are not affected by clouds). -Since it is difficult to understand the change in seasonality in a single year, I tried to acquire and visualize the data for 15 years. -If you also display the moving average in layers, you can see the seasonal changes in the vegetation index more prominently.

Changes in vegetation index (15 years)

##Satellite data analysis ③~Ground surface temperature~ ###Overview

• "temperature"Is the temperature in the air, whereas"Ground surface temperature"Represents the temperature of the surface of the earth
• Land Surface Temperature(LST)Called -Artificial structures (buildings and roads) tend to have higher surface temperature than forests and soil. -It is used for monitoring the heat island phenomenon, etc.

###data set -satellite= "Terra(MODIS)" -band= LST (Land Surface Temerature) -resolution: 1000m -area:Tokyo (mainly the Imperial Palace)

###code -Try to get LST by changing various parameters -Visualize the acquired data in the same way

## Parameter specification
# Designate satellite
snippet = 'MODIS/006/MOD11A2'
# Specify the band name
band = 'LST_Day_1km'

# Specify the period
from_date='2005-01-01'
to_date='2005-12-31'
# Specify area (specify latitude / longitude for Japan area)
geometry = ee.Geometry.Rectangle([139.686, 35.655, 139.796, 35.719])

# Folder name to save
dir_name = 'GEE_LST'
# resolution
scale = 1000

## Execution of processing ------------------------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## Save all (File name uses satellite ID)
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

# Data reading
with rasterio.open('/content/drive/My Drive/GEE_LST/MODIS_006_MOD11A2_2005_08_05.tif') as src:
    arr = src.read()

# Visualization
plt.figure(figsize=(20, 8))
plt.subplot(121); plt.imshow(rgb); plt.title('RGB Image'); plt.title('Optical Image')
plt.subplot(122); plt.imshow(arr[0], cmap='inferno'); plt.title('Land Surface Temperature')

Optical RGB image and ground surface temperature in Tokyo (around the Imperial Palace) -It is a little difficult to understand because the resolution is low compared to other data. -Looking at all the applicable periods, it seems that the temperature tends to be higher in the west side. -One of the features of the ground surface temperature is that while the time resolution is high (data is acquired frequently), there are many images for which data cannot be acquired correctly due to the influence of clouds.

Changes in surface temperature in Tokyo (around the Imperial Palace) (January)~December)

Changes in surface temperature in Tokyo (around the Imperial Palace) (for 15 years)

-I also acquired long-term data as well as the vegetation index and visualized it including the moving average. -Like vegetation, it seems to be seasonal (higher in summer and lower in winter) -Since the resolution is low, it seems better to calculate the surface temperature over a wider area and analyze the global trend.

##at the end -Introduced typical satellite image analysis using Google Earth Engine and Google Colab -By using GEE and Colab as described above,**Various satellite data can be analyzed simply by changing the variable name of the satellite name or band name.I think you can see -AlsoAfter acquiring satellite data, you can use Python's convenient analysis library in the same environment.**Is very convenient -In the actual analysis, detailed preprocessing is added according to the requirements and a model is created to remove the bias, but by using these services, the initial introduction step is greatly streamlined. I think I can do it -Also,Analysis of uploading business data at hand to Google Drive and merging it with satellite data on ColabSeems easy to do -I hope this article will help you to utilize artificial satellite data in various fields.