[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.

Screenshot from 2020-09-18 18-48-17.png

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.

Screenshot from 2020-09-18 17-27-26.png

About the observation image during the initial function confirmation period of the Advanced Land Observing Satellite-2 "DAICHI-2" (ALOS-2) (Part 5 / Ship Edition)

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.

Screenshot from 2020-08-16 11-17-18.png

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. Screenshot from 2020-09-18 17-36-23.png

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

Screenshot from 2020-09-18 17-38-23.png

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()

Screenshot from 2020-09-18 18-22-25.png

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')

Screenshot from 2020-09-18 18-24-36.png

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()

Screenshot from 2020-09-18 18-25-53.png

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.

Cell count

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')

Screenshot from 2020-09-18 18-30-40.png

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()

Screenshot from 2020-09-18 18-33-24.png

Enlarge the part where many ships were imaged.

plt.imshow(distmap[0:200,0:200])
plt.colorbar()
plt.show()

Screenshot from 2020-09-18 18-34-31.png

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)

Screenshot from 2020-09-18 18-41-32.png

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()

Screenshot from 2020-09-18 18-48-17.png

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

Cell count

Recommended Posts

I tried to find the trend of the number of ships in Tokyo Bay from satellite images.
I tried fitting the exponential function and logistics function to the number of COVID-19 positive patients in Tokyo
I tried to predict the number of people infected with coronavirus in consideration of the effect of refraining from going out
I tried to find the entropy of the image with python
How to find the optimal number of clusters in k-means
I tried to find the average of the sequence with TensorFlow
Implementation of recommendation system ~ I tried to find the similarity from the outline of the movie using TF-IDF ~
I tried to sort out the objects from the image of the steak set meal-② Overlap number sorting
(Python) I tried to analyze 1 million hands ~ I tried to estimate the number of AA ~
I tried to display the altitude value of DTM in a graph
[Completed version] Try to find out the number of residents in the town from the address list with Python
I tried to refactor the template code posted in "Getting images from Flickr API with Python" (Part 2)
I tried to find 100 million digits of pi
I tried to touch the API of ebay
Find the number of days in a month
I tried to correct the keystone of the image
I tried to predict the price of ETF
I tried to vectorize the lyrics of Hinatazaka46!
I tried to predict the number of people infected with coronavirus in Japan by the method of the latest paper in China
python beginners tried to predict the number of criminals
I want to detect images of cats from Instagram
How to get the number of digits in Python
I tried to detect the iris from the camera image
I tried to summarize the basic form of GPLVM
I tried to visualize the spacha information of VTuber
I tried to erase the negative part of Meros
I tried to classify the voices of voice actors
I tried to summarize the string operations of Python
[Linux] I learned LPIC lv1 in 10 days and tried to understand the mechanism of Linux.
I tried to sort out the objects from the image of the steak set meal-④ Clustering
I used Python to find out about the role choices of the 51 "Yachts" in the world.
Find a guideline for the number of processes / threads to set in the application server
I tried to find out the outline about Big Gorilla
[Horse Racing] I tried to quantify the strength of racehorses
Maya | Find out the number of polygons in the selected object
I tried to get the location information of Odakyu Bus
I tried the accuracy of three Stirling's approximations in python
I tried to create API list.csv in Python from swagger.yaml
Examine the margin of error in the number of deaths from pneumonia
I tried to tabulate the number of deaths per capita of COVID-19 (new coronavirus) by country
[For beginners] I want to explain the number of learning times in an easy-to-understand manner.
[Deep Learning from scratch] I tried to explain the gradient confirmation in an easy-to-understand manner.
I tried to get the number of days of the month holidays (Saturdays, Sundays, and holidays) with python
Python --Find out number of groups in the regex expression
[Python] Programming to find the number of a in a character string that repeats a specified number of times.
I tried to summarize the code often used in Pandas
I tried to illustrate the time and time in C language
How to increase the number of machine learning dataset images
[Python] I tried to visualize the follow relationship of Twitter
I tried to create a Python script to get the value of a cell in Microsoft Excel
I tried to summarize the commands often used in business
I tried to implement the mail sending function in Python
[Machine learning] I tried to summarize the theory of Adaboost
I tried to find the affine matrix in image alignment (feature point matching) using affine transformation.
I tried to predict the genre of music from the song title on the Recurrent Neural Network
I tried to fight the Local Minimum of Goldstein-Price Function
I tried changing the python script from 2.7.11 to 3.6.0 on windows10
I tried to sort out the objects from the image of the steak set meal-① Object detection
I tried to put HULFT IoT (Edge Streaming) in the gateway Rooster of Sun Electronics
[Natural language processing] I tried to visualize the remarks of each member in the Slack community
I tried to get various information from the codeforces API