[PYTHON] Let's take a look at the forest fire on the west coast of the United States with satellite images.

(Update) Updated the method of creating composite images.

Overview

So far, we have introduced the acquisition of satellite images of the Sentinel series of European artificial satellites using Google Earth Engine and their analysis methods. This time, we will introduce the acquisition and analysis method of satellite images of Sentinel-2, which is an optical observation satellite familiar to everyone in the Sentinel series. This time, let's take a look at the forest fires on the west coast of the United States reported in overseas news. 　

This is a satellite image of a forest fire site near Portland on the west coast of the United States (taken on September 9, 2020).

The following is reported on overseas news sites. This site has a history of forest fires, and you can see where the fire broke out. 　Fire Maps: Glass and Zogg Wildfire Tracker

As for the situation of the site, the following sites may be real and squeezed. The New York Times has a lot of interesting articles that I read every day. 　Historic Wild Fires Range in Western States

Here, we will introduce an example of how to acquire the observed image of Sentinel-2 with Google Earth Engine and how to convert the image data into a color image. Also, the code introduced in the article is placed on Github, so if you want to try it yourself, please use this. ** Google colaboratory is used **, so you can use it regardless of PC as long as you have a network environment **.

1.First of all

Previously, using Sentinel's API, I acquired an image captured by an artificial satellite of an Australian forest fire in early 2020 and created a GIF animation of it.

I made a time-lapse of satellite images of Australian forest fires.

Here, the standard 100km square image of Sentinel-2 was acquired using the API, and the part of interest was cut out and color composition processing was performed. At this time, the scene will be different depending on the imaging timing, so the image has some notches. Even if they are different scenes, you can make an image of the entire area of interest by acquiring and merging each of them. However, handling multiple large size images is not very efficient considering the acquisition time. Therefore, this time, we will acquire the image cut out in the area of interest from Google Earth Engine (hereinafter referred to as GEE) and create a color composite image. 　 The Sentinel series is an earth observation series developed and operated by the European Space Agency, and has already been launched and operated from 1 to 5. This time, we will introduce image processing from the acquisition of satellite images of Sentinel-2, which captures optical observation images as well as digital camera images.

Sentinel-2 ©ESA 2000-2018.

For details on the observation plan and observation images of Sentinel-2, see How to obtain the latest satellite images for free. .

Sentinel-2 captures images of 12 different wavelengths.

(Credit: European Space Agency)

A general satellite image (True color) can be created by setting bands 4, 3, and 2 as R, G, and B.

This time, we will focus on forest fires and try to visualize the fires using infrared wavelength observation images. This time, I refer to the following articles.

How to use open source satellite data for your investigative reportingHow to use open source satellite data for your investigative reporting

From this article, we will look at forest fire sites using images of 12 and 11, which are long-wavelength bands that are sensitive to fire, and 8A, which is sensitive to forests.

2. Acquisition and evaluation of satellite images

2.1 Environment preparation (construction).

Please refer to the following article for details on how to acquire satellite images with GEE. Artificial satellite image analysis by Google Earth Engine and Google Colab-Satellite image analysis starting for free (Introduction)-

For the analysis method of satellite images obtained from GEE, please refer to the following article.

I checked the usage status of the parking lot from satellite images.

Although there are some overlaps, we will introduce the acquisition of satellite images and their processing.

import ee
import numpy as np
import matplotlib.pyplot as plt

ee.Authenticate()
ee.Initialize()

First, execute this to authenticate the GEE connection. When you execute it, the link will be returned. Click on it to perform the authentication procedure, copy the access code and enter it.

Next, authenticate the connection to Google Drive. Again, the flow is the same as GEE certification.

from google.colab import drive
drive.mount('/content/drive')

Next, we will perform work such as viewing the acquired satellite image and installing the modules necessary for digitizing and analyzing it.

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

import json
import os
import glob

import time
from datetime import datetime
from dateutil.parser import parse

Frequently used modules are already installed in Google Colaboratory, so no additional work is required, but this time we will use Geotiff, which is an image with map information added, so it is necessary for image processing ** Rasterio Install **.

Next, install a module called ** folium ** to check the set target area on the map.

!pip install folium

import folium

Now that the environment is ready, take a satellite image from GEE.

2.2 Setting the area of interest (target area)

In order to acquire satellite images from GEE, you need to enter the latitude / longitude information of the target area you are interested in. This time as well, we will use the following site created to easily find and obtain the ** latitude / longitude ** of the area of interest.

#Acquisition of polygon information of the area of interest.
from IPython.display import HTML
HTML(r'<iframe width="1000" height="580" src="https://gispolygon.herokuapp.com/" frameborder="0"></iframe>')

When you execute this, the following screen will be displayed.

After expanding the area of interest, select a square polygon from the icon on the left to display the polygon of the area of interest. Then click ** Show Features ** to see the polygon geographic information in the right window. Then click ** Copy ** at the bottom to copy this geographic information. This time, we will acquire satellite images of ** the eastern area of Portland ** in Oregon, USA.

Then paste the copied map information below and enter it.

A  = {"type":"FeatureCollection","features":[{"properties":{"note":"","distance":"210025.76 m","drawtype":"rectangle","area":"520197.34 ha"},"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[-122.60860204696657,44.80966973731369],[-122.60860204696657,45.400920907537866],[-121.60266637802125,45.400920907537866],[-121.60266637802125,44.80966973731369],[-122.60860204696657,44.80966973731369]]]}}]}

This geographic information is processed into the input format for GEE and for display in Folium.

#Set any file name to be used in the future. For example, the name of the area.
object_name = 'Portland'

with open(str(object_name) +'_2.geojson', 'w') as f:
    json.dump(A, f)

json_file = open(str(object_name) +'_2.geojson')
json_object = json.load(json_file)

#Extract only the latitude / longitude information of the area of interest from json.
AREA = json_object["features"][0]["geometry"]['coordinates'][0]

Now, let's check the set area of interest.

m = folium.Map([(AREA[0][1]+AREA[len(AREA)-2][1])/2,(AREA[0][0]+AREA[len(AREA)-3][0])/2], zoom_start=9)

folium.GeoJson(str(object_name) +'_2.geojson').add_to(m)
m

output

2.3 Acquisition of satellite images from GEE

A lot of satellite images and a lot of information already analyzed are set in GEE. For details, please refer to Data Catalog. Sentinel-1 and 2 are as follows.

Sentinel-2 MSI: MultiSpectral Instrument, Level-1C

From this page, you can use the data of ** Sentinel-2 observation image from June 23, 2015 **.

Now, get the Sentinel-2 image from GEE and save it in Google colaboratory.

First, prepare the format of the geographic information to be set in GEE.

region=ee.Geometry.Polygon(AREA)

Next, set the parameters of the information to be acquired. This time, the period of the acquired image and the save destination of the acquired image are specified. This time, we confirmed that we were observing the area on September 9, so we set the period to acquire only that image.

#Specify the period
from_date='2020-09-08'
to_date='2020-09-10'

#Folder name to save
dir_name_s2 = 'GEE_Sentinel2_' + object_name

Now, let's set the image conditions for Sentinel-2.

def ImageExport(image,description,folder,region,scale):
    task = ee.batch.Export.image.toDrive(image=image,description=description,folder=folder,region=region,scale=scale)
    task.start()

Sentinel2 = ee.ImageCollection('COPERNICUS/S2').filterBounds(region).filterDate(parse(from_date),parse(to_date)).filterMetadata('CLOUDY_PIXEL_PERCENTAGE','less_than', 80).select(['B12','B11','B8A'])

imageList_s2 = Sentinel2.toList(300) 

Here, we set to acquire the color composite images of ** B12 **, ** B11 ** and ** B8A ** that are used for color composition among the observation images of Sentinel-2. In addition, although it is possible to filter the proportion of clouds in the Sentinel-2 image, it is set to 80%, which is higher than usual, in order to avoid misrecognizing smoke from a fire as a cloud. Now, under the above satellite image acquisition conditions, the image of the Sentinel-2 area of interest is acquired.

for i in range(imageList_s2.size().getInfo()):
    image = ee.Image(imageList_s2.get(i))
    ImageExport(image.reproject(crs='EPSG:4326',scale=10),image.get('system:index').getInfo(),dir_name_s2,region['coordinates'][0],10)

2.4 Display of satellite images

Display and confirm the acquired satellite image. Satellite images are saved in the directory (folder) set in my drive of Google drive. Call it and display it.

#Visualization in chronological order
s2_path = '/content/drive/My Drive/' + dir_name_s2 + '/'
files =os.listdir(s2_path)
files.sort()

plt.figure(figsize=(25, 25))
j=0

v = len(files)//5 +1 
for i in range(len(files)):
  #Acquire and visualize images one scene at a time
  with rasterio.open(s2_path + files[i]) as src:
      arr = src.read()
  j+=1#Shift and place the plot position of the image
  plt.subplot(v,5,j)
  arrayImg = np.asarray(arr).transpose(1,2,0).astype(np.float16)/10000
  plt.imshow(arrayImg)
  plt.title(files[i][0:8])#Get the date from the file name
  plt.tight_layout()

We acquired four images of the area. The reason why it is divided into four is that the image of Sentinel-2 is divided and registered in 100km square, and this area of interest contains the four divided images. By synthesizing the first and second images, the area of interest is covered, so create a ** composite image ** of these two images. First, check each image.

#Check the image
s2_path = '/content/drive/My Drive/' + dir_name_s2 + '/'
files =os.listdir(s2_path)
files.sort()

#Data reading
n = 0

with rasterio.open(s2_path + files[n]) as src:
    arr = src.read()

print(files[n][0:8])
#Visualization
arrayImg0 = np.asarray(arr).transpose(1,2,0).astype(np.float16)/10000
plt.imshow(arrayImg0)

#Check the image
s2_path = '/content/drive/My Drive/' + dir_name_s2 + '/'
files =os.listdir(s2_path)
files.sort()

#Data reading
n = 1

with rasterio.open(s2_path + files[n]) as src:
    arr = src.read()

print(files[n][0:8])
#Visualization
arrayImg1 = np.asarray(arr).transpose(1,2,0).astype(np.float16)/10000
plt.imshow(arrayImg1)

Next, combine these two images. ** (Update) ** I think there are multiple methods for merging images, but I used to merge with numpy before, but I will combine them with the merge command of rasterio that was advised in the comment.

import rasterio.merge

dest, out_transform = rasterio.merge.merge([s2_path + files[0], s2_path + files[1]])

The composite image is converted by the transpose function because of the order of RGB channel, Vertical, and Horizontal.

image = dest/10000

img_transformed = image.transpose((1, 2, 0))
img_transformed.shape

output

(6596, 11199, 3)

Output the image.

plt.imshow(img_transformed)

You can synthesize it well. (Before the update, there was a black part at the joint.) Now, save this composite image as jpeg.

import cv2

new_image = ((img_transformed - img_transformed.min()) * (1/(img_transformed.max() - img_transformed.min()) * 255)).astype('uint8')


im_rgb = cv2.cvtColor(new_image, cv2.COLOR_BGR2RGB)
cv2.imwrite(str(object_name) +'.jpg', im_rgb )

img = Image.open(str(object_name) +'.jpg')

plt.figure(figsize=(15, 20))

plt.imshow(img)
plt.show()

Next, the date of imaging of the artificial satellite and the credit of the image are put on this image.

2.5 Enter the imaging date and credit on the satellite image.

Next, the shooting date and credit of the satellite image are described in this image. The fonts used here are those that are open to the public for free on the net. This is also your favorite.

#Font file download and settings
!wget https://osdn.net/dl/mplus-fonts/mplus-TESTFLIGHT-063a.tar.xz

!xz -dc mplus-TESTFLIGHT-*.tar.xz | tar xf -

fontfile = "./mplus-TESTFLIGHT-063a/mplus-1c-bold.ttf"

Now, write the characters on the image.

date = files[n][0:8]

img = Image.open(str(object_name) +'.jpg')

img = img.convert('RGB')

x = int(img.size[0]/1.3) #Setting the date description position
y = int(img.size[1]/20) #Setting the date description position
fs = int(img.size[0]/70) #Date font size setting

obj_draw = ImageDraw.Draw(img)
obj_font = ImageFont.truetype(fontfile, fs)
obj_draw.text((x, y), 'Observation Date: '+str(date), fill=(255, 255, 255), font=obj_font)
obj_draw.text((img.size[0]/1.6, img.size[1]-y*0.1 - img.size[1]/30 ), 'Contains modified Copernicus Sentinel data (2020)', fill=(255, 255, 255), font=obj_font)

img = img.resize((int(img.size[0] / 2) , int(img.size[1] / 2)))

img.save(str(object_name) +'.jpg')

plt.figure(figsize=(15, 20))

plt.imshow(img)
plt.show()

Now you have the necessary information in the image. Change the above code to adjust the character position, size, color, etc. to your liking. The created image file can be downloaded to your PC by right-clicking the target file name from the File list on the left and selecting ** Download **.

3. Finally

Using the ** Google Earth Engine ** provided by Google, we introduced the method of acquiring satellite images of ** Sentinel-2 ** and image processing of ** forest fires ** on the west coast of the United States as an analysis example. .. Forest fires are still ongoing and great damage has been reported. I pray that the activities of the fire brigade will be fruitful and that they will soon return to their normal lives. If you have any comments or questions, please feel free to comment. I'm happy.

Reference article

How to get the latest satellite imagery for free. I tried to automatically acquire the observation image of the artificial satellite (Sentinel-2) using API. Create an estimated distribution map of cars from artificial satellite images by PyTorch. I checked the usage status of the parking lot from satellite images. I made a time-lapse of satellite images of Australian forest fires.

[Sequel] Artificial satellite image analysis by Google Earth Engine and Google Colab-Satellite image analysis starting for free (practice)- Artificial satellite image analysis by Google Earth Engine and Google Colab-Satellite image analysis starting for free (Introduction)- [6. Local GEE execution by Python](https://sites.google.com/site/mizuochipublic/%E5%AE%9F%E8%B7%B5%E3%82%B3%E3%83%B3 % E3% 83% 86% E3% 83% B3% E3% 83% 84-remote-sensing-tutorials / google-earth-engine% E5% 85% A5% E9% 96% 80 / 6-python% E3% 81 % AB% E3% 82% 88% E3% 82% 8B% E3% 83% AD% E3% 83% BC% E3% 82% AB% E3% 83% AB% E3% 81% 8B% E3% 82% 89 % E3% 81% AEgee% E5% AE% 9F% E8% A1% 8C)

Kihon of satellite data-Understanding, type, frequency, resolution, application examples ~ Can you see people from artificial satellites? ~ By satellite, ground resolution / local time summary ~

Sentinel-2 MSI: MultiSpectral Instrument, Level-1C