LoginSignup
8
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@t-mat

JupyterLabで衛星画像を処理してみた。

Last updated at Posted at 2020-01-03

1.概要

@nigo1973さんの「無料で最新の衛星画像を入手する方法」というエントリを拝見し、欧州宇宙機関(ESA: Eurpean Space Agency)が地球観測衛星Sentinel-1とSentinel-2の観測画像を無料配布していることを知りました。
とても興味を持ったので、年末年始等の時間を利用してSentinel-2の観測画像をダウンロードし、画像処理をしてみました。
「SNAP」という衛星画像処理ソフトが欧州宇宙機関より配布されていますが、後々、「Tellus」でも色々試してみたいと思い、今回はJupyterLabを使用しました。

2.衛星画像のダウンロードと参考資料

衛星画像は「無料で最新の衛星画像を入手する方法」を参考に「Copernicus Open Data Hub」からダウンロードしました。
なお、今回は関東周辺の2019/10/30観測レベル1Cプロダクトをダウンロードしています。
Sentinel-2のプロダクトについては「Sentinel User Guide」の「product-types」を参考ください。
その他、Sentinel-2の情報は、(一財)リモート・センシング技術センターの「衛星総覧」を、GDAL/OGRについては、「 osgeo - Python Api」、「Python GDAL/OGR Cookbook」を参考にしました。

3.GDALでの衛星画像の読み込み

今回使用したモジュールは下記のとおりです。
なお、Jpeg2000画像を扱うので、OpenJPEGのインストールが必要です。

import numpy as np
from osgeo import gdal,gdalconst
from osgeo import osr 
from PIL import Image,ImageOps
from colour import Color
import matplotlib.pyplot as plt 
%matplotlib inline
Image.MAX_IMAGE_PIXELS = 1000000000

ダウンロードしたファイルを解凍し、パスを設定します。

path='S2B_MSIL1C_20191030T012759_N0208_R074_T54SUE_20191030T030910.SAFE\
      GRANULE\L1C_T54SUE_A013828_20191030T012755\IMG_DATA/T54SUE_20191030T012759_'

下記で解像度10ｍの赤・緑・青・近赤外バンド、解像度20ｍの赤・中間赤外バンドの画像を読み込みます。

#解像度10ｍのバンド
r_band = gdal.Open(path+"B04.jp2", gdalconst.GA_ReadOnly) 
g_band = gdal.Open(path+"B03.jp2", gdalconst.GA_ReadOnly)
b_band = gdal.Open(path+"B02.jp2", gdalconst.GA_ReadOnly)
nir_band=gdal.Open(path+"B08.jp2", gdalconst.GA_ReadOnly)
#解像度20ｍのバンド
r2_band = gdal.Open(path+"B05.jp2", gdalconst.GA_ReadOnly)
lir_band = gdal.Open(path+"B11.jp2", gdalconst.GA_ReadOnly)

読み込んだファイルの概要を確認してみます。
画像サイズは10980×10980ピクセル（解像度10ｍ）、座標系はWGS84/UTM 54N（EPGS：32654）であることがわかります。
また、平面座標系でのアフィン変換を確認します。

print(r_band.RasterXSize,r_band.RasterYSize,r_band.RasterCount)
> 10980 10980 1

print(r_band.GetProjection())
> PROJCS["WGS 84 / UTM zone 54N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,
> AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
> UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],
> PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",141],
> PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],
> UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32654"]]

print(r_band.GetGeoTransform())
> (300000.0, 10.0, 0.0, 4000020.0, 0.0, -10.0)

4.衛星画像をGeotifに出力する

1)カラーコンポジット画像の出力

トゥルーカラー、フォールスカラー、ナチュナルカラーのカラー合成画像をGeotiffで出力するため、以下の関数を作成しました。
なお、トゥルーカラー等についてはWikipedia「衛星画像」を参考ください。

def outputGeotifMulti(src,epgs,band1,band2,band3,path):
    tfw=src.GetGeoTransform()
    dtype = gdal.GDT_UInt16
    band = 3
    width=src.RasterXSize
    height=src.RasterYSize
    output = gdal.GetDriverByName('GTiff').Create(path, width, height, band, dtype, options = ['PHOTOMETRIC=RGB'])
    output.SetGeoTransform(tfw)
    crs = osr.SpatialReference()
    crs.ImportFromEPSG(epgs) # 32654 WGS84 UTM_54N
    output.SetProjection(crs.ExportToWkt())
    output.GetRasterBand(1).WriteArray(band1)
    output.GetRasterBand(2).WriteArray(band2)
    output.GetRasterBand(3).WriteArray(band3)
    output.FlushCache()
    output=None

下記で各画像のRGBにバンドを割り当て、Geotiffを出力します。

r  =r_band.GetRasterBand(1).ReadAsArray()
g  =g_band.GetRasterBand(1).ReadAsArray()
b  =b_band.GetRasterBand(1).ReadAsArray() 
nir=nir_band.GetRasterBand(1).ReadAsArray() 
outputGeotifMulti(r_band,32653,r,g,b,'img_true.tif')
outputGeotifMulti(r_band,32653,nir,r,g,'img_false.tif')
outputGeotifMulti(r_band,32653,r,nir,g,'img_natural.tif')

出力したGeotifをQGISに読み込ませると、以下の画像が表示されます。

2)NDVI、NDWIの出力

次いで、NDVI（正規化植生指標）とNDWI（正規化水指数）を計算し、Geotif画像を出力します。
NDVI、NDWIについては、「宙畑」の記事「課題に応じて変幻自在？ 衛星データをブレンドして見えるモノ・コト #マンガでわかる衛星データ」がとても参考になりました。

def calNDVI(nir_band,r_band):
    nir = nir_band.astype('float32')
    red = r_band.astype('float32')
    ndvi=np.where(
        (nir+red)==0., 
        0,(nir-red)/(red+nir))
    return ndvi

def calNDWI(lir_band,r_band):
    lir = lir_band.astype('float32')
    red = r_band.astype('float32')
    ndwi=np.where(
        (lir+red)==0., 
        0,(red-lir)/(red+lir))
    return ndwi

def outputGeotifSingle(src,epgs,band1,path):
    tfw=src.GetGeoTransform()
    dtype = gdal.GDT_Float32
    band = 1
    width=src.RasterXSize
    height=src.RasterYSize
    output = gdal.GetDriverByName('GTiff').Create(path, width, height, band, dtype)
    output.SetGeoTransform(tfw)
    crs = osr.SpatialReference()
    crs.ImportFromEPSG(epgs)
    output.SetProjection(crs.ExportToWkt())
    output.GetRasterBand(1).WriteArray(band1)
    output.FlushCache()
    output=None

下記でシングルバンドのGeotif画像を出力します。

ndvi=calNDVI(nir,r)
ndwi=calNDWI(lir,r2)
outputGeotifSingle(r_band,32653,ndvi,'img_ndvi.tif')
outputGeotifSingle(r2_band,32653,ndwi,'img_ndwi.tif')

出力したGeotif画像をQGISで単バンド疑似カラーで表示すると以下のとおりです。

5.衛星画像をJpg/Jgwで出力する

1)カラーコンポジット画像の出力

Geotifで出力すると、カラーコンポジット画像で約700MBと容量が大きく、個人が趣味で扱うにはちょっとキツイサイズとなります。
このため、下記の関数を定義して、画像をuint8のRGB（Jpg）で出力し、座標情報はワールドファイルに出力してみます。
ワールドファイルについては、農林水産研究情報総合センターの「ラスタデータ用の座標ファイル（ワールドファイル）について」を参考ください。

def createRGBImage(band1,band2,band3):
    b1=band1.astype('float32')
    b2=band2.astype('float32')
    b3=band3.astype('float32')
    b1=(b1/b1.max()*255).astype('uint8')
    b2=(b2/b2.max()*255).astype('uint8')
    b3=(b3/b3.max()*255).astype('uint8')
    color=np.dstack((b1,b2,b3))
    img = Image.fromarray(color)
    return ImageOps.equalize(img)

def outputWF(src,path):
    tfw=src.GetGeoTransform()
    with open(path, mode='w') as f:
        f.write(str(tfw[1])+'\n')
        f.write('0\n0\n')
        f.write(str(tfw[5])+'\n')
        f.write(str(tfw[0])+'\n')
        f.write(str(tfw[3])+'\n')

def outputJpgMulti(src,band1,band2,band3,path_name):
    img=createRGBImage(band1,band2,band3)
    img.save(path_name+'.jpg',quality=95)
    outputWF(src,path_name+'.jgw')

下記で各画像のRGBにバンドを割り当て、Jpg画像とjgwファイル（ワールドファイル）を出力します。

outputJpgMulti(r_band,r,g,b,'img_true')
outputJpgMulti(r_band,nir,r,g,'img_false')
outputJpgMulti(r_band,r,nir,g,'img_natural')

出力したJpg画像は以下のとおりです。
容量はGeotifの10分の1程度になりました。
なお、画像生成時にImageOps.equalizeメソッドを使い、RGBヒストグラムを補正しています。

2)NDVI、NDWIの出力

次に、NDVIとNDWIのJpg画像を生成してみます。
colourモジュールでグラデーション配列を生成し、各ピクセルのNDVI値、NDWI値からRGB値を割り当ててJpg画像を生成します。

def createRGBImageSingle(band,color):
    min=band.min()
    dif=band.max()-min
    height,width=band.shape
    col=[]
    for c in color:
        col.append([int(255*c.rgb[0]),int(255*c.rgb[1]),int(255*c.rgb[2])])
    length=len(col)-1
    col=np.array(col)
    rgb=np.zeros((height,width,3)).astype('uint8')
    for y in range(height):
        for x in range(width):
            val=int(((band[y,x]-min)/dif)*length);
            rgb[y,x]=col[val]
    img = Image.fromarray(rgb)
    return ImageOps.equalize(img)

def outputJpgSingle(src,band,color,path_name):
    img=createRGBImageSingle(band,color)
    img.save(path_name+'.jpg',quality=95)
    outputWF(src,path_name+'.jgw')


colors1 = list(Color('#ffffff').range_to(Color("#00ff00"),10))
outputJpgSingle(r2_band,ndvi,colors1,'img_ndvi')

colors2 = list(Color('#ffffff').range_to(Color("#0000ff"),10))
outputJpgSingle(r2_band,ndwi,colors2,'img_ndwi')

出力したJpg画像は以下のとおりです。

6.最後に

個人で衛星画像を取得し遊べる訳ですから、すごい時代になったものです。
11月に参加した「Tellus Bootcamp」のテキストがとても参考になりました。
Sentinelの他、Landsatの衛星画像も「landsatlook」から検索・ダウンロードできるようです。
日本の衛星画像も「G-Portal」から検索・ダウンロードできるようですが、まだ、試していません。

8
6
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
8
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?