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

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

## 2.衛星画像のダウンロードと参考資料
衛星画像は「[無料で最新の衛星画像を入手する方法](https://qiita.com/nigo1973/items/9bb6a11caac8e3e1e850)」を参考に「[Copernicus Open Data Hub](https://scihub.copernicus.eu/dhus/#/home)」からダウンロードしました。
なお、今回は関東周辺の2019/10/30観測レベル1Cプロダクトをダウンロードしています。
Sentinel-2のプロダクトについては「[Sentinel User Guide](https://earth.esa.int/web/sentinel/user-guides)」の「[product-types](https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/product-types)」を参考ください。
その他、Sentinel-2の情報は、(一財)リモート・センシング技術センターの「[衛星総覧](https://www.restec.or.jp/satellite/sentinel-2-a-2-b)」を、GDAL/OGRについては、「 [osgeo -  Python Api](https://gdal.org/python/index.html)」、「[Python GDAL/OGR Cookbook](https://pcjericks.github.io/py-gdalogr-cookbook/)」を参考にしました。

## 3.GDALでの衛星画像の読み込み
今回使用したモジュールは下記のとおりです。
なお、Jpeg2000画像を扱うので、[OpenJPEG](https://www.openjpeg.org/)のインストールが必要です。

```python
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
```
ダウンロードしたファイルを解凍し、パスを設定します。

```python
path='S2B_MSIL1C_20191030T012759_N0208_R074_T54SUE_20191030T030910.SAFE\
      GRANULE\L1C_T54SUE_A013828_20191030T012755\IMG_DATA/T54SUE_20191030T012759_'
```
下記で解像度10ｍの赤・緑・青・近赤外バンド、解像度20ｍの赤・中間赤外バンドの画像を読み込みます。

```python
#解像度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）であることがわかります。
また、平面座標系でのアフィン変換を確認します。

```python
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「[衛星画像](https://ja.wikipedia.org/wiki/%E8%A1%9B%E6%98%9F%E7%94%BB%E5%83%8F)」を参考ください。

```python
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を出力します。

```python
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に読み込ませると、以下の画像が表示されます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/f3040b55-2630-48b5-93f4-4d46f859870b.jpeg"  height="300" />

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

```python
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画像を出力します。

```python
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で単バンド疑似カラーで表示すると以下のとおりです。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/76d7b6f4-a541-5523-e8c7-6f7b6ff042ff.jpeg"  height="300" />

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

```python
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ファイル（ワールドファイル）を出力します。

```python
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ヒストグラムを補正しています。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/b725c669-f2e5-b3db-c808-3b286b2b1905.jpeg"  height="300" />
### 2)NDVI、NDWIの出力
次に、NDVIとNDWIのJpg画像を生成してみます。
[colour](https://www.colour-science.org/)モジュールでグラデーション配列を生成し、各ピクセルのNDVI値、NDWI値からRGB値を割り当ててJpg画像を生成します。

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

```


```python
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画像は以下のとおりです。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/698c728b-877f-af8c-d636-2eca60999490.jpeg"  height="300" />

## 6.最後に
個人で衛星画像を取得し遊べる訳ですから、すごい時代になったものです。
11月に参加した「[Tellus Bootcamp](https://sorabatake.jp/7696/)」のテキストがとても参考になりました。
Sentinelの他、Landsatの衛星画像も「[landsatlook](https://landsatlook.usgs.gov/viewer.html)」から検索・ダウンロードできるようです。
日本の衛星画像も「[G-Portal](https://gportal.jaxa.jp/gpr/)」から検索・ダウンロードできるようですが、まだ、試していません。