-
Notifications
You must be signed in to change notification settings - Fork 0
/
sentinel-test30-timelapse.py
386 lines (297 loc) · 13.5 KB
/
sentinel-test30-timelapse.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
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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 IPython.display import HTML
HTML(r'<iframe width="960" height="480" src="https://www.keene.edu/campus/maps/tool/" frameborder="0"></iframe>')
#上記にて取得した地理情報をコピペする
AREA = [
[
-220.291841,
35.6593884
],
[
-220.2932143,
35.4817801
],
[
-220.1380324,
35.4817801
],
[
-220.1421523,
35.6493456
],
[
-220.291841,
35.6593884
]
]
print(AREA)
print(len(AREA))
print(AREA[0])
print(AREA[0][0])
# Convert coordinate to 360degree format
for i in range(len(AREA)):
AREA[i][0] = AREA[i][0] +360
# Check the result before the conversion
print(AREA)
print(len(AREA))
print(AREA[0])
print(AREA[0][0])
# Convert coordinate to Polygon format
m=Polygon([AREA])
print(m)
#ファイル名を定義. 好きな名称に設定してください.
object_name = 'sentinel-test30-timelapse'
# Convert the Polygon data to GeoJSON format
with open(str(object_name) +'.geojson', 'w') as f:
json.dump(m, f)
# Convert a GeoJSON object to Well-Known Text. Intended for use with OpenSearch queries.
footprint_geojson = geojson_to_wkt(read_geojson(str(object_name) +'.geojson'))
print(footprint_geojson)
# Sentinel API credentials
user = 'username'
password = 'password'
api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus')
# Create a base map, simply by passing the starting coordinates to Folium
# In this case get the coordinate from the first and last x,y divided by 2
x_map = (AREA[0][1]+AREA[len(AREA)-1][1])/2
y_map = (AREA[0][0]+AREA[len(AREA)-1][0])/2
xyzoom_start = 11
m = folium.Map([x_map, y_map], zoom_start=xyzoom_start)
# Add the GeoJSON data (coordinate info) to folium
folium.GeoJson(str(object_name) +'.geojson').add_to(m)
# Display the Folium map on OpenStreetMap
m
# Dont forget to download the font file from https://osdn.net/dl/mplus-fonts/mplus-TESTFLIGHT-063a.tar.xz
# wget https://osdn.net/dl/mplus-fonts/mplus-TESTFLIGHT-063a.tar.xz
# ls -l mplus-TESTFLIGHT-063a.tar.xz
# -rw-r--r-- 1 jovyan users 10371708 Apr 23 2019 mplus-TESTFLIGHT-063a.tar.xz
!xz -dc mplus-TESTFLIGHT-*.tar.xz | tar xf -
cwd = os.getcwd()
suffix = '.ttf'
base_filename = 'mplus-1c-bold'
fontfile = os.path.join(cwd,'mplus-TESTFLIGHT-063a',base_filename + suffix)
def Sentinel2_get_init(i):
products = api.query(footprint_geojson,
date = (Begin_date, End_date1), #画像取得月の画像取得
platformname = 'Sentinel-2',
processinglevel = 'Level-2A', #Leve-1C
cloudcoverpercentage = (0,100)) #被雲率（0%〜50%）
products_gdf = api.to_geodataframe(products)
#画像取得月の被雲率の小さい画像からソートする．
products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
title = products_gdf_sorted.iloc[i]["title"]
print(str(title))
uuid = products_gdf_sorted.iloc[i]["uuid"]
product_title = products_gdf_sorted.iloc[i]["title"]
#画像の観測日を確認．
date = products_gdf_sorted.iloc[i]["ingestiondate"].strftime('%Y-%m-%d')
print(date)
#Sentinel-2 data download
api.download(uuid)
file_name = str(product_title) +'.zip'
#ダウロードファイルの解凍．
with zipfile.ZipFile(file_name) as zf:
zf.extractall()
#フォルダの画像ファイルのアドレスを取得．
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]) +'/R10m'
files3 = os.listdir(pathB)
path_b2 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B02_10m.jp2')
path_b3 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B03_10m.jp2')
path_b4 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B04_10m.jp2')
# Band４，３，２をR,G,Bとして開く．
b4 = rio.open(path_b4)
b3 = rio.open(path_b3)
b2 = rio.open(path_b2)
#RGBカラー合成（GeoTiffファイルの出力）
with rio.open(str(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()
#ポリゴン情報の取得．
nReserve_geo = gpd.read_file(str(object_name) +'.geojson')
#取得画像のEPSGを取得
epsg = b4.crs
nReserve_proj = nReserve_geo.to_crs({'init': str(epsg)})
#マスクTiffファイルの一時置き場
os.makedirs('./Image_tiff', exist_ok=True)
#カラー合成画像より関心域を抽出．
with rio.open(str(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_' +str(object_name) +'.tif', "w", **out_meta) as dest:
dest.write(out_image)
#抽出画像のjpeg処理．
scale = '-scale 0 250 0 30'
options_list = [
'-ot Byte',
'-of JPEG',
scale
]
options_string = " ".join(options_list)
#ディレクトリの作成
os.makedirs('./Image_jpeg_'+str(object_name), exist_ok=True)
#jpeg画像の保存．
gdal.Translate('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(object_name) +'.jpg',
'./Image_tiff/Masked_' +str(object_name) +'.tif',
options=options_string)
#画像への撮像日を記載
img = Image.open('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(object_name) +'.jpg')
#print(img.size)
#print(img.size[0])
x = img.size[0]/100 #日付の記載位置の設定
y = img.size[1]/100 #日付の記載位置の設定
fs = img.size[0]/50 #日付のフォントサイズの設定
fs1 = int(fs)
#print(fs1)
#print(type(fs1))
obj_draw = ImageDraw.Draw(img)
obj_font = ImageFont.truetype(fontfile, fs1)
obj_draw.text((x, y), str(date), fill=(255, 255, 255), font=obj_font)
obj_draw.text((img.size[0]/2, img.size[1]-y - img.size[1]/20 ), 'produced from ESA remote sensing data', fill=(255, 255, 255), font=obj_font)
img.save('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(object_name) +'.jpg')
#大容量のダウンロードファイル，および解凍フォルダの削除
shutil.rmtree( str(product_title) + '.SAFE')
os.remove(str(product_title) +'.zip')
return
def Sentinel2_get():
products = api.query(footprint_geojson,
date = (Begin_date, End_date), #取得希望期間の入力
platformname = 'Sentinel-2',
processinglevel = 'Level-2A', #Leve-1C
cloudcoverpercentage = (0,100)) #被雲率（０％〜５０％）
products_gdf = api.to_geodataframe(products)
products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
#同一シーンの画像を取得するため，placenumberを固定する．
products_gdf_sorted = products_gdf_sorted[products_gdf_sorted["title"].str.contains(placenumber)]
title = products_gdf_sorted.iloc[0]["title"]
print(str(title))
uuid = products_gdf_sorted.iloc[0]["uuid"]
product_title = products_gdf_sorted.iloc[0]["title"]
date = products_gdf_sorted.iloc[0]["ingestiondate"].strftime('%Y-%m-%d')
print(date)
api.download(uuid)
file_name = str(product_title) +'.zip'
with zipfile.ZipFile(file_name) as zf:
zf.extractall()
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]) +'/R10m'
files3 = os.listdir(pathB)
path_b2 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B02_10m.jp2')
path_b3 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B03_10m.jp2')
path_b4 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B04_10m.jp2')
b4 = rio.open(path_b4)
b3 = rio.open(path_b3)
b2 = rio.open(path_b2)
with rio.open(str(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()
os.makedirs('./Image_tiff', exist_ok=True)
nReserve_geo = gpd.read_file(str(object_name) +'.geojson')
epsg = b4.crs
nReserve_proj = nReserve_geo.to_crs({'init': str(epsg)})
with rio.open(str(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_' +str(object_name) +'.tif', "w", **out_meta) as dest:
dest.write(out_image)
from osgeo import gdal
scale = '-scale 0 250 0 30'
options_list = [
'-ot Byte',
'-of JPEG',
scale
]
options_string = " ".join(options_list)
os.makedirs('./Image_jpeg_'+str(object_name), exist_ok=True)
gdal.Translate('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(object_name) +'.jpg',
'./Image_tiff/Masked_' +str(object_name) +'.tif',
options=options_string)
#画像への撮像日の記載
img = Image.open('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(object_name) +'.jpg')
#print(img.size)
#print(img.size[0])
x = img.size[0]/100 #日付の記載位置の設定
y = img.size[1]/100 #日付の記載位置の設定
fs = img.size[0]/50 #日付のフォントサイズの設定
fs1 = int(fs)
#print(fs1)
#print(type(fs1))
obj_draw = ImageDraw.Draw(img)
obj_font = ImageFont.truetype(fontfile, fs1)
obj_draw.text((x, y), str(date), fill=(255, 255, 255), font=obj_font)
obj_draw.text((img.size[0]/2, img.size[1]-y - img.size[1]/20 ), 'produced from ESA remote sensing data', fill=(255, 255, 255), font=obj_font)
img.save('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(object_name) +'.jpg')
shutil.rmtree( str(product_title) + '.SAFE')
os.remove(str(product_title) +'.zip')
return
#L2Aデータ取得のため，２０２０年１月以降を開始日とすること．
Begin_date = '20200501'
End_date1 = Begin_date[:6] +'30'#月の終わりを２８としているのは，２月を考慮
Sentinel2_get_init(0)
placenumber = '_T54SUE_' #上記で取得したplacenumberを貼り付けます
#L2Aデータ取得の場合は、2019年1月以降を開始日とすること
Begin_date = '20190701'
End_date = '20200731'
#年をまたぐ画像を取得する場合のため、以下のコードにて月数をカウントする
d_year = int(End_date[2:4]) - int(Begin_date[2:4])
d_year = d_year*12 + int(End_date[4:6]) - int(Begin_date[4:6])
# 各月の衛星画像の取得し、関心領域のjpeg画像を作成
for i in range(d_year):
if i < 1:
m = int(Begin_date[4:6])
else:
m = int(Begin_date[4:6]) +1
if m <10:
Begin_date = Begin_date[:4] +'0'+ str(m) + Begin_date[6:]
elif m <13:
Begin_date = Begin_date[:4] + str(m) + Begin_date[6:]
else:
y = int(Begin_date[2:4]) +1
Begin_date = Begin_date[:2] + str(y) + '01' + Begin_date[6:]
End_date = Begin_date[:6] +'28'
Sentinel2_get()
# 作成したjpeg画像を、古い観測画像よりGIFアニメーションにする
images =[]
files = sorted(glob.glob('./Image_jpeg_'+str(object_name) +'/*.jpg'))
images = list(map(lambda file: Image.open(file), files))
images[0].save('./Image_jpeg_'+str(object_name) +'/' + str(object_name) + '.gif', save_all=True, append_images=images[1:], duration=1000, loop=0)