storing your sentinel-1 data as COG

Python, Sentinel 1 Leave a comment

Cloud Optimized GeoTIFF ftw

 #!/usr/bin/python

import os
import gdal
import numpy as np

# Specify import
inVV = "../sentinel1Data_S1A_IW_GRDH_1SDV_20180201T115551_20180201T115616_020413_022E26_7CDE_Gamma0_VV_mst_01Feb2018.tif"
inVH = "../sentinel1Data_S1A_IW_GRDH_1SDV_20180201T115551_20180201T115616_020413_022E26_7CDE_Gamma0_VH_mst_01Feb2018.tif"

# Specify output file
outFileName = "../outputFile.tif"

# Dictionary with all meta data
myDict = {"geometry": {"coordinates": [[[91.210281, 23.560665], [93.703537, 23.987156], [93.994301, 22.480576], [91.529236, 22.050838], [91.210281, 23.560665]]], "type": "Polygon"}, "id": 0, "properties": {"acquisitiontype": "NOMINAL", "beginposition": "2018-02-01T11:55:51.548000Z", "endposition": "2018-02-01T11:56:16.547000Z", "filename": "S1A_IW_GRDH_1SDV_20180201T115551_20180201T115616_020413_022E26_7CDE.SAFE", "format": "SAFE", "id": "4a250148-9ed1-4709-b893-ede81d4da033", "identifier": "S1A_IW_GRDH_1SDV_20180201T115551_20180201T115616_020413_022E26_7CDE", "ingestiondate": "2018-02-01T15:48:48.007000Z", "instrumentname": "Synthetic Aperture Radar (C-band)", "instrumentshortname": "SAR-C SAR", "lastorbitnumber": 20413, "lastrelativeorbitnumber": 41, "link": "https://scihub.copernicus.eu/dhus/odata/v1/Products('4a250148-9ed1-4709-b893-ede81d4da033')/$value", "link_alternative": "https://scihub.copernicus.eu/dhus/odata/v1/Products('4a250148-9ed1-4709-b893-ede81d4da033')/", "link_icon": "https://scihub.copernicus.eu/dhus/odata/v1/Products('4a250148-9ed1-4709-b893-ede81d4da033')/Products('Quicklook')/$value", "missiondatatakeid": 142886, "orbitdirection": "ASCENDING", "orbitnumber": 20413, "platformidentifier": "2014-016A", "platformname": "Sentinel-1", "polarisationmode": "VV VH", "productclass": "S", "producttype": "GRD", "relativeorbitnumber": 41, "sensoroperationalmode": "IW", "size": "1.61 GB", "slicenumber": 6, "status": "ARCHIVED", "summary": "Date: 2018-02-01T11:55:51.548Z, Instrument: SAR-C SAR, Mode: VV VH, Satellite: Sentinel-1, Size: 1.61 GB", "swathidentifier": "IW", "title": "S1A_IW_GRDH_1SDV_20180201T115551_20180201T115616_020413_022E26_7CDE", "uuid": "4a250148-9ed1-4709-b893-ede81d4da033"}, "type": "Feature"}

# store properties in a single dictionary
d = {}
for item in myDict['properties']:
	print str(item), str(myDict['properties'][item])
	d.update({str(item) : str(myDict['properties'][item])})

# get VV
dsVV = gdal.Open(inVV)
bandVV = dsVV.GetRasterBand(1)
arrVV = bandVV.ReadAsArray()

# get VH
dsVH = gdal.Open(inVH)
bandVH = dsVH.GetRasterBand(1)
arrVH = bandVH.ReadAsArray()

# get # columns and rows
[cols, rows] = arrVV.shape

# multiply by 10000 to store as integer
arrVH = arrVH*10000
arrVV = arrVV*10000

# filter for edges
arrVH = np.where((arrVH < 0.0001), 0, arrVH)
arrVV = np.where((arrVV < 0.0001), 0, arrVV)

# specify output driver and projection
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(outFileName, rows, cols, 2, gdal.GDT_UInt16,
                         options=["TILED=YES",
                                  "COMPRESS=LZW",
                                  "INTERLEAVE=BAND"])
outdata.SetGeoTransform(dsVV.GetGeoTransform())
outdata.SetProjection(dsVV.GetProjection())

# write the first band
outdata.GetRasterBand(1).WriteArray(arrVH)
outdata.GetRasterBand(1).SetNoDataValue(0)
outdata.GetRasterBand(1).SetDescription("VH")

# write the second band
outdata.GetRasterBand(2).WriteArray(arrVV)
outdata.GetRasterBand(2).SetNoDataValue(0)
outdata.GetRasterBand(2).SetDescription("VV")

# build overview and set metadata
outdata.BuildOverviews("NEAREST", [2, 4, 8, 16, 32, 64])
outdata.SetMetadata(d)

# write the file
outdata.FlushCache()
outdata = None
band=None
ds=None

f

Leave a Reply