-
Notifications
You must be signed in to change notification settings - Fork 0
/
sentinel2.py
265 lines (213 loc) · 10.3 KB
/
sentinel2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
# IMPORTANT: Make sure you install below packages in advance
#!pip install sentinelsat
#!pip install rasterio
#!apt install gdal-bin python-gdal python3-gdal
#!apt install python3-rtree
#!pip install git+git://github.com/geopandas/geopandas.git
#!pip install descartes
#!pip install shapely
#!pip install six
#!pip install pyproj
#!pip install descartes
#!pip install geopandas
#!pip install folium
import os
import numpy as np
from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
from shapely.geometry import MultiPolygon, Polygon
from rasterio.plot import show
from geojson import Polygon
from osgeo import gdal
from PIL import Image, ImageDraw, ImageFont
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import json
import shutil
import glob
import rasterio as rio
import rasterio.mask
import fiona
import folium
import zipfile
from notebook import notebookapp
import urllib
import ipykernel
#Generate coordinate of Area of Interest by using tools from Keene University below
#from IPython.display import HTML
#HTML(r'<iframe width="960" height="480" src="https://www.keene.edu/campus/maps/tool/" frameborder="0"></iframe>')
# Sentinel API credentials
user = 'username'
password = 'password'
api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus')
def sentinel2_hello():
print("Hello from Sentinel2")
return
# (JupyterNotebook only) Get notebookpath
def notebook_path():
"""Returns the absolute path of the Notebook or None if it cannot be determined
NOTE: works only when the security is token-based or there is also no password
"""
connection_file = os.path.basename(ipykernel.get_connection_file())
kernel_id = connection_file.split('-', 1)[1].split('.')[0]
for srv in notebookapp.list_running_servers():
try:
if srv['token']=='' and not srv['password']: # No token and no password, ahem...
req = urllib.request.urlopen(srv['url']+'api/sessions')
else:
req = urllib.request.urlopen(srv['url']+'api/sessions?token='+srv['token'])
sessions = json.load(req)
for sess in sessions:
if sess['kernel']['id'] == kernel_id:
return os.path.join(srv['notebook_dir'],sess['notebook']['path'])
except:
pass # There may be stale entries in the runtime directory
return None
# Convert the Polygon data to GeoJSON format
def Sentinel2_convert_polygon_to_json(object_name, polygon_object):
print("Converting polygon to GeoJSON...")
with open(object_name +'.geojson', 'w') as f:
json.dump(polygon_object, f)
print(object_name +".geojson created")
return
## Get Sentinel Geodataframe
def Sentinel2_get_productsgdf(footprint_geojson, start_date, end_date):
print("Start Date: ", start_date)
print("End Date: ", end_date)
print("Querying Sentinel-2...")
products = api.query(footprint_geojson,
date = (start_date, end_date), #Desired date for the beginning and ending time of Sentinel-2 image
platformname = 'Sentinel-2',
processinglevel = 'Level-2A',
cloudcoverpercentage = (0,70))
products_gdf = api.to_geodataframe(products)
products_gdf_sorted_by_date = products_gdf.sort_values(['ingestiondate'], ascending=[True])
return products_gdf_sorted_by_date
## Get Sentinel all scenes
def Sentinel2_get_all_data(products_gdf, index, fontfile, object_name, start_date, resolution):
print("\n---START---")
products_gdf = products_gdf
title = products_gdf.iloc[index]["title"]
print(" Sentinel-2 title: \t", str(title))
uuid = products_gdf.iloc[index]["uuid"]
product_title = products_gdf.iloc[index]["title"]
ingestion_date = products_gdf.iloc[index]["ingestiondate"].strftime('%Y%m%d-%H:%M')
print(" IngestionDate: \t", ingestion_date)
cloudcoverpercentage = products_gdf.iloc[index]["cloudcoverpercentage"]
print(" Cloud percentage: \t", cloudcoverpercentage)
datasize = products_gdf.iloc[index]["size"]
print(" Data size: \t", datasize)
#Download Sentinel-2 data
print(" Download data now...")
api.download(uuid)
file_name = str(product_title) +'.zip'
print(" Successfully downloaded data as: \t", file_name)
#Extract downloaded Sentinel-2 data
print("\nExtracting data.zip file...")
with zipfile.ZipFile(file_name) as zf:
zf.extractall()
#Get image's folder path
print("\nFile identification...")
print(" Resolution: \t", resolution)
path = str(product_title) + '.SAFE/GRANULE'
files = os.listdir(path)
pathA = str(product_title) + '.SAFE/GRANULE/' + str(files[0])
files2 = os.listdir(pathA)
pathB = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R' + resolution + 'm'
files3 = os.listdir(pathB)
path_b2 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R' + resolution + 'm/' +str(files3[0][0:23] +'B02_'+resolution+'m.jp2')
path_b3 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R' + resolution + 'm/' +str(files3[0][0:23] +'B03_'+resolution+'m.jp2')
path_b4 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R' + resolution + 'm/' +str(files3[0][0:23] +'B04_'+resolution+'m.jp2')
print(" ", path_b2)
print(" ", path_b3)
print(" ", path_b4)
#Open Band4(Blue), 3(Green) and 2(Red)
b4 = rio.open(path_b4)
b3 = rio.open(path_b3)
b2 = rio.open(path_b2)
#RGB color compose (output file as GeoTiff: public domain metadata standard which allows georeferencing information to be embedded within a TIFF file)
print("\nRGB Color composing...")
with rio.open(object_name +'.tiff','w',driver='Gtiff', width=b4.width, height=b4.height,
count=3,crs=b4.crs,transform=b4.transform, dtype=b4.dtypes[0]) as rgb:
rgb.write(b4.read(1),1)
rgb.write(b3.read(1),2)
rgb.write(b2.read(1),3)
rgb.close()
#Read polygon from .geojson
nReserve_geo = gpd.read_file(object_name +'.geojson')
# Reference https://rasterio.readthedocs.io/en/latest/api/rasterio.crs.html#module-rasterio.crs
# EPSGの326544というのは、地図投影する際によく利用される、WGS84のUTM座標系のことを表します。
epsg = b4.crs
# Since this RGB image is large and huge you save both computing power and time to clip and use only the area of interest.
# We will clip the Natural reserve area from the RGB image.
nReserve_proj = nReserve_geo.to_crs({'init': str(epsg)})
#Create directory for temporary Masked Tiff
print("Creating directory Image_tiff...")
os.makedirs('./Image_tiff', exist_ok=True)
#Extract image for Area of Interest from the composed color image
print("\nCalculating area of interest...")
with rio.open(object_name +'.tiff') as src:
out_image, out_transform = rio.mask.mask(src, nReserve_proj.geometry,crop=True)
out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
with rasterio.open('./Image_tiff/' +'Masked_' + object_name +'.tif', "w", **out_meta) as dest:
dest.write(out_image)
#jpeg processing from the extracted images
scale = '-scale 0 250 0 30'
options_list = [
'-ot Byte',
'-of JPEG',
scale
]
options_string = " ".join(options_list)
#Create directory for jpeg processed image
print("Creating Image_jpeg...")
os.makedirs('./Image_jpeg_'+object_name, exist_ok=True)
#Save jpeg image
#https://pypi.org/project/GDAL/
#GDAL Geospatial Data Abstraction Library
gdal.Translate('./Image_jpeg_'+object_name +'/' + str(ingestion_date) + 'Masked_' +object_name +'.jpg',
'./Image_tiff/Masked_' +object_name +'.tif',
options=options_string)
#Print the ingestion_date on the image
img = Image.open('./Image_jpeg_'+object_name +'/' + str(ingestion_date) + 'Masked_' +object_name +'.jpg')
#print(img.size)
#print(img.size[0])
x = img.size[0]/100 #x coordinate for the print position
y = img.size[1]/100 #y coordinate for the print position
fs = img.size[0]/50 #font size
fs1 = int(fs)
#print(fs1)
#print(type(fs1))
obj_draw = ImageDraw.Draw(img)
obj_font = ImageFont.truetype(fontfile, fs1)
obj_draw.text((x, y), '-- 75 Tahun INDONESIA MERDEKA! --\n' + str(ingestion_date) + 'produced from European Space Agency (ESA) remote sensing data', fill=(255, 255, 255), font=obj_font)
#obj_draw.text((img.size[0]/2, img.size[1]-y - img.size[1]/20 ), 'produced from European Space Agency (ESA) remote sensing data', fill=(255, 255, 255), font=obj_font)
img.save('./Image_jpeg_'+str(object_name) +'/' + str(ingestion_date) + 'Masked_' +object_name +'.jpg')
#Remove downloaded files
print("\nRemoving files...")
shutil.rmtree( str(product_title) + '.SAFE')
os.remove(str(product_title) +'.zip')
print("---DONE---\n")
return
## Get Sentinel satellite scene
def Sentinel2_get_sorted_data(index, fontfile, object_name, footprint_geojson, start_date, end_date, resolution):
print("Querying Sentinel-2...")
products = api.query(footprint_geojson,
date = (start_date, end_date), #Desired date for the beginning and ending time of Sentinel-2 image
platformname = 'Sentinel-2',
processinglevel = 'Level-2A',
cloudcoverpercentage = (0,100))
products_gdf = api.to_geodataframe(products)
#Sort the value of cloud coverage percentage from the small one
print("Sorting data from the lowest Cloud Coverage Percentage...")
products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
print(" Data sorted! Information as follow:")
#Get Sentinel-2 data
Sentinel2_get_all_data(products_gdf_sorted, index, fontfile, object_name, start_date, resolution)
print("Get sorted data finished.")
return