[PYTHON] I tried to find the trend of the number of ships in Tokyo Bay from satellite images.
Overview
From the satellite image obtained from the previous Google Earth Engine, we calculated the transition of the number of cars in the parking lot of Tokyo Disneyland.
I checked the usage status of the parking lot from satellite images.
This time, as in the previous time, we will obtain satellite images from Google Earth Engine and find the changes in the number of vessels in Tokyo Bay.
This is a count of ships in the images taken by artificial satellites near the entrance to Tokyo Bay, and the trends are calculated. I will introduce it in this article including the explanation of this graph. 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
Last time, we calculated the number of cars in the parking lot of Tokyo Disneyland using the radio wave sensor (SAR) image in the satellite image.
I checked the usage status of the parking lot from satellite images.
It is difficult to identify the type of vehicle such as a car or truck in the SAR image, but it can be qualitatively evaluated if a large number of cars are parked in the parking lot, and it is not affected by clouds. Therefore, it can be evaluated stably. Therefore, we will use the SAR image again to find the number of vessels from the satellite image near the entrance of Tokyo Bay. This is because I thought that it would be possible to evaluate how the number of vessels has changed due to the influence of the new coronavirus, whether it has not changed, and the extent to which it has changed. I think this site will be helpful for how ships can be seen in SAR images.
This captured image of DAICHI-2 is an image in observation mode with a resolution of 3 m, and if you are a large ship, you can see its shape. The important point here is that the ship appears bright and the surrounding sea appears dark. SAR satellites irradiate the earth with radio waves and image the intensity of the reflected signal. At this time, ** radio waves are radiated from an oblique direction with respect to the ground surface **, so if the ground surface is flat, many things are reflected in the incident direction of the radio waves, that is, in the direction opposite to the satellite that receives the radiated radio waves. To do. Therefore, ** If the target ground surface is flat, the received signal is weak and the image will be dark. The surface of water such as lakes and the sea is extremely flat, so it looks darker **. On the other hand, when there are structures such as buildings, buildings, cars and ships, the objects are uneven, so more radio waves are reflected in the same direction as the incident direction, so the image looks bright. Let's count the number of ships in the satellite image, using this ** brightly reflected part ** as a ship.
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 obtain it using the following site created to easily find and obtain the latitude and 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 this area near the entrance to Tokyo Bay.
Then paste the copied map information below and enter it.
A = {"type":"FeatureCollection","features":[{"properties":{"note":"","distance":"35937.06 m","drawtype":"rectangle","area":"12713.15 ha"},"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[139.71197873353958,35.33475055335104],[139.71197873353958,35.452996930987666],[139.81840878725052,35.452996930987666],[139.81840878725052,35.33475055335104],[139.71197873353958,35.33475055335104]]]}}]}
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 = 'tokyobay3'
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=15)
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-1 SAR GRD: C-band Synthetic Aperture Radar Ground Range Detected, log scaling
From this page, you can use the data of ** Sentinel-1 observation image from October 3, 2014 **.
Now, get the Sentinel-1 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.
#Specify the period
from_date='2019-01-01'
to_date='2020-08-31'
#Folder name to save
dir_name_s1 = 'GEE_Sentinel1_' + object_name
Now, let's set the image conditions for Sentinel-1.
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()
Sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD').filterBounds(region).filterDate(parse(from_date),parse(to_date)).filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING')).select(['VV'])
imageList_s1 = Sentinel1.toList(300)
Here, since the observation image of Sentinel-1 uses only the observation image of satellite B at 6:00 pm,'ASCENDING' is selected in'orbitProperties_pass'. When set to "Descending", the observation image will be at 6:00 am.
Now, under the above satellite image acquisition conditions, the image of the Sentinel-1 area of interest is acquired.
for i in range(imageList_s1.size().getInfo()):
image = ee.Image(imageList_s1.get(i))
ImageExport(image.reproject(crs='EPSG:4326',scale=10),image.get('system:index').getInfo(),dir_name_s1,region['coordinates'][0],10)
2.4 Display of satellite images and evaluation of signal strength
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
s1_path = '/content/drive/My Drive/' + dir_name_s1 + '/'
files =os.listdir(s1_path)
files.sort()
plt.figure(figsize=(20, 40))
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(s1_path + files[i]) as src:
arr = src.read()
j+=1#Shift and place the plot position of the image
plt.subplot(v,5,j)
plt.imshow(arr[0], cmap='gray')
plt.title(files[i][33:41])#Get the date from the file name
plt.tight_layout()
The observed image of the SAR satellite is a black and white image, and the bright part in this image seems to be a ship. Next, count the number of bright parts in the image as ships.
#Data reading
n = 2
with rasterio.open(s1_path + files[n]) as src:
arr = src.read()
print(files[n][33:41])
#Visualization
plt.imshow(arr[0], cmap='gray')
When I magnified the satellite image, it was confirmed that there were many bright spots. Now, let's find the intensity distribution of this image.
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.hist(arr[0], bins=50)
fig.show()
The signal strength is often 0 or less, but some signals of 0 or more have also been observed. From this condition, the binarization process is performed first with 0 or more and less than 0.
2.5 Calculate the number of ships from satellite images.
For the counting method of objects in the image, I referred to the following site.
This method binarizes with an arbitrary signal strength, finds the center position if there is a spread, checks whether the center position is the same object or not in an arbitrary window, and determines the number of objects. Count. Then, I will introduce the method concretely. First, the sample image was the latest captured image in the acquired image.
#Data reading
n = len(files) -1
with rasterio.open(s1_path + files[n]) as src:
arr = src.read()
print(files[n][33:41])
#Visualization
plt.imshow(arr[0], cmap='gray')
You can see that there are many bright spots. Next, binarization processing is performed. Here, the binarization process is performed with 0 or more and less than 0 with reference to the above signal strength distribution. In addition, a map of the center position is created from the obtained binarized image.
#Binarization process
import cv2
binimg = (arr[0]>0) #Let 0 be the threshold.
plt.imshow(binimg)
plt.colorbar()
plt.show()
#Distance map(distance map)To calculate
binimg = binimg.astype(np.uint8)
distmap = cv2.distanceTransform(binimg,1,3)
plt.imshow(distmap)
plt.colorbar()
plt.show()
Enlarge the part where many ships were imaged.
plt.imshow(distmap[0:200,0:200])
plt.colorbar()
plt.show()
You can see that the signal that seems to be a ship extends in a certain direction. Next, the distribution with the maximum signal strength is obtained from this image. At this time, the size of the ship was set to 30 pixels. In other words, we define that there is only one ship in 30 pixels.
#Calculate the location with the maximum value of the distance map
out = distmap*0
ksize=30 #Maximum value within 30 pixels
for x in range(ksize,distmap.shape[0]-ksize*2):
for y in range(ksize,distmap.shape[1]-ksize*2):
if distmap[x,y]>0 and distmap[x,y]==np.max(distmap[x-ksize:x+ksize,y-ksize:y+ksize]):
out[x,y]=1
plt.imshow(distmap[0:200,0:200])
plt.colorbar()
plt.show()
Next, if the signal strength is high at a pixel near the position that seems to be a ship, it may be falsely detected, so if it is within a certain range, it will be ignored.
#Inflate, contour detect and count
out = cv2.dilate(out,(10,10)) #10*Inflate to 10 and do not count if within the same range
contours, _ = cv2.findContours(out.astype(np.uint8), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
arr=[]
for i in contours:
x_=0
y_=0
for j in i:
x_ += j[0][0]
y_ += j[0][1]
arr.append([x_/len(i), y_/len(i)])
arr = np.array(arr)
plt.imshow(distmap[0:200,0:200])
plt.colorbar()
plt.show()
print("Number of ships: ", len(arr))
print(arr)
Here, the number of ships and their positions are output. Using this method, the number of ships is calculated from the acquired satellite image, and the transition is graphed.
#Time series graph of the estimated number of ships in the area
#Visualization in chronological order
s1_path = '/content/drive/My Drive/' + dir_name_s1 + '/'
files =os.listdir(s1_path)
files.sort()
sum_ship = []
label_signal = []
for i in range(len(files)):
label_signal.append(files[i][33:41])
#Acquire and visualize images one scene at a time
with rasterio.open(s1_path + files[i]) as src:
arr = src.read()
#Binarization process
binimg = (arr[0]>0)
#Distance map(distance map)To calculate
binimg = binimg.astype(np.uint8)
distmap = cv2.distanceTransform(binimg,1,3)
out = distmap*0
ksize=30 #Maximum value within 30 pixels
for x in range(ksize,distmap.shape[0]-ksize*2):
for y in range(ksize,distmap.shape[1]-ksize*2):
if distmap[x,y]>0 and distmap[x,y]==np.max(distmap[x-ksize:x+ksize,y-ksize:y+ksize]):
out[x,y]=1
out = cv2.dilate(out,(10,10)) #10*Inflate to 10 and do not count if within the same range
contours, _ = cv2.findContours(out.astype(np.uint8), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
num=[]
for i in contours:
x_=0
y_=0
for j in i:
x_ += j[0][0]
y_ += j[0][1]
num.append([x_/len(i), y_/len(i)])
num = np.array(num)
sum_ship.append(len(num))
#Visualization
fig,ax = plt.subplots(figsize=(15,6))
plt.plot(sum_ship, marker='o')
ax.set_xticks(np.arange(0,len(files)))
ax.set_xticklabels(label_signal, rotation=90)
plt.show()
The horizontal axis is the observation date, and the vertical axis is the number of vessels calculated from the observation image. It will reach its maximum value around March 2019, but after that it will average around 30 vessels. I wondered if the movement by ship was restricted or attenuated in the corona disaster after March 2020, but as far as I can see, it seems that the number of parking lots in Tokyo Tisneyland is extremely decreasing. not. In fact, I searched the Internet to see how many ships in Tokyo Bay were sailing, and found an article about cruise ships refraining from operating, but many other ships are declining. There was nothing. I think that everyone was operating while paying attention to safety in order to stabilize the distribution even in the corona virus. Thank you. I tried to investigate the places where many ships are navigating other than Tokyo Bay by the same method, but there is no tendency that the number of coronaviruses decreases drastically like this time, and the people who are operating ships Through my efforts, I realized that stable logistics are being maintained.
3. Finally
Using the Google Earth Engine provided by Google, we introduced how to calculate the number of vessels in Tokyo Bay as an example of satellite image acquisition and analysis. I hope this will be an opportunity for more people to be interested in satellite imagery. 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.
[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-1 SAR GRD: C-band Synthetic Aperture Radar Ground Range Detected, log scaling