HUB Datacube Cookbook

General

Is HUB Datacube installed

Imports HUB Datacube and exits the program if the modules are not found.

import sys
try:
    import hubdc.core, hubdc.applier
except:
    sys.exit('ERROR: cannot find HUB Datacube modules')

Is the testdata installed

In this guide we use the EnMAP-Box testdata (https://bitbucket.org/hu-geomatics/enmap-box-testdata).

import sys
try:
    import enmapboxtestdata
except:
    sys.exit('ERROR: cannot find EnMAP-Box Testdata modules')

Is QGIS installed (optional)

We use QGIS 3 map canvas for quick looks (https://qgis.org).

try:
    import qgis, qgis.core, qgis.gui
except:
    exit('ERROR: cannot find QGIS modules')

Check versions installed

import hubdc
import enmapboxtestdata
import qgis.utils
print(hubdc.__version__)
print(enmapboxtestdata.__version__)
print(qgis.utils.Qgis.QGIS_VERSION)
0.16.0
0.9
3.4.3-Madeira

Raster dataset

Find API reference here: hubdc.core.RasterDataset

Open a raster dataset from file

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
print(rasterDataset)

Prints something like:

RasterDataset(gdalDataset=<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x000001CD7B004AE0> >)

Open a raster dataset from GDAL dataset

import enmapboxtestdata
from hubdc.core import *
from osgeo import gdal

# init with GDAL Dataset
gdalDataset = gdal.Open(enmapboxtestdata.enmap)
rasterDataset = RasterDataset(gdalDataset=gdalDataset)
print(rasterDataset)

# get the GDAL Dataset handle from a RasterDataset
print(rasterDataset.gdalDataset())

Prints something like:

RasterDataset(gdalDataset=<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x00000110C38B4B40> >)
<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x00000110C38B4B40> >

Close a raster dataset

Closing a raster dataset is useful in the middle of a script, to recover the resources held by accessing the dataset, remove file locks, etc. It is not necessary at the end of the script, as the Python garbage collector will do the same thing automatically when the script exits.

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
rasterDataset.close()
print(rasterDataset)

Prints:

RasterDataset(gdalDataset=None)

Get raster metadata

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)

# get all domains
print('All domains:')
print(rasterDataset.metadataDict())

# get ENVI domain
print('ENVI domain:')
print(rasterDataset.metadataDomain(domain='ENVI'))

# get single item (list with wavelength casted to float)
print('ENVI wavelength:' )
print(rasterDataset.metadataItem(key='wavelength', domain='ENVI', dtype=float))

Prints:

All domains:
{'IMAGE_STRUCTURE':
    {
    'INTERLEAVE': 'BAND'
    },
 'ENVI':
    {
    'bands': '177',
    'band names': ['band 8', 'band 9', 'band 10', 'band 11', 'band 12', 'band 13', 'band 14', 'band 15', 'band 16', 'band 17', 'band 18', 'band 19', 'band 20', 'band 21', 'band 22', 'band 23', 'band 24', 'band 25', 'band 26', 'band 27', 'band 28', 'band 29', 'band 30', 'band 31', 'band 32', 'band 33', 'band 34', 'band 35', 'band 36', 'band 37', 'band 38', 'band 39', 'band 40', 'band 41', 'band 42', 'band 43', 'band 44', 'band 45', 'band 46', 'band 47', 'band 48', 'band 49', 'band 50', 'band 51', 'band 52', 'band 53', 'band 54', 'band 55', 'band 56', 'band 57', 'band 58', 'band 59', 'band 60', 'band 61', 'band 62', 'band 63', 'band 64', 'band 65', 'band 66', 'band 67', 'band 68', 'band 69', 'band 70', 'band 71', 'band 72', 'band 73', 'band 74', 'band 75', 'band 76', 'band 77', 'band 91', 'band 92', 'band 93', 'band 94', 'band 95', 'band 96', 'band 97', 'band 98', 'band 99', 'band 100', 'band 101', 'band 102', 'band 103', 'band 104', 'band 105', 'band 106', 'band 107', 'band 108', 'band 109', 'band 110', 'band 111', 'band 112', 'band 113', 'band 114', 'band 115', 'band 116', 'band 117', 'band 118', 'band 119', 'band 120', 'band 121', 'band 122', 'band 123', 'band 124', 'band 125', 'band 126', 'band 127', 'band 144', 'band 145', 'band 146', 'band 147', 'band 148', 'band 149', 'band 150', 'band 151', 'band 152', 'band 153', 'band 154', 'band 155', 'band 156', 'band 157', 'band 158', 'band 159', 'band 160', 'band 161', 'band 162', 'band 163', 'band 164', 'band 165', 'band 166', 'band 167', 'band 168', 'band 195', 'band 196', 'band 197', 'band 198', 'band 199', 'band 200', 'band 201', 'band 202', 'band 203', 'band 204', 'band 205', 'band 206', 'band 207', 'band 208', 'band 209', 'band 210', 'band 211', 'band 212', 'band 213', 'band 214', 'band 215', 'band 216', 'band 217', 'band 218', 'band 219', 'band 220', 'band 221', 'band 222', 'band 223', 'band 224', 'band 225', 'band 226', 'band 227', 'band 228', 'band 229', 'band 230', 'band 231', 'band 232', 'band 233', 'band 234', 'band 235', 'band 236', 'band 237', 'band 238', 'band 239'],
    'byte order': '0',
    'coordinate system string': ['PROJCS["UTM_Zone_33N"', 'GEOGCS["GCS_WGS_1984"', 'DATUM["D_WGS_1984"', 'SPHEROID["WGS_1984"', '6378137.0', '298.257223563]]', 'PRIMEM["Greenwich"', '0.0]', 'UNIT["Degree"', '0.0174532925199433]]', 'PROJECTION["Transverse_Mercator"]', 'PARAMETER["False_Easting"', '500000.0]', 'PARAMETER["False_Northing"', '0.0]', 'PARAMETER["Central_Meridian"', '15.0]', 'PARAMETER["Scale_Factor"', '0.9996]', 'PARAMETER["Latitude_Of_Origin"', '0.0]', 'UNIT["Meter"', '1.0]]'],
    'data ignore value': '-99',
    'data type': '2',
    'description': ['EnMAP02_Berlin_Urban_Gradient_2009.bsq', 'http://doi.org/10.5880/enmap.2016.008', 'spectral and spatial subset'],
    'file compression': '1',
    'file type': 'ENVI Standard',
    'fwhm': ['0.005800', '0.005800', '0.005800', '0.005800', '0.005800', '0.005800', '0.005800', '0.005800', '0.005800', '0.005800', '0.005900', '0.005900', '0.006000', '0.006000', '0.006100', '0.006100', '0.006200', '0.006200', '0.006300', '0.006400', '0.006400', '0.006500', '0.006600', '0.006600', '0.006700', '0.006800', '0.006900', '0.006900', '0.007000', '0.007100', '0.007200', '0.007300', '0.007300', '0.007400', '0.007500', '0.007600', '0.007700', '0.007800', '0.007900', '0.007900', '0.008000', '0.008100', '0.008200', '0.008300', '0.008400', '0.008400', '0.008500', '0.008600', '0.008700', '0.008700', '0.008800', '0.008900', '0.008900', '0.009000', '0.009100', '0.009100', '0.009200', '0.009300', '0.009300', '0.009400', '0.009400', '0.009500', '0.009500', '0.009600', '0.009600', '0.009600', '0.009600', '0.009700', '0.009700', '0.009700', '0.011800', '0.011900', '0.012100', '0.012200', '0.012400', '0.012500', '0.012700', '0.012800', '0.012900', '0.013100', '0.013200', '0.013300', '0.013400', '0.013500', '0.013600', '0.013700', '0.013800', '0.013900', '0.014000', '0.014000', '0.014100', '0.014100', '0.014200', '0.014200', '0.014300', '0.014300', '0.014300', '0.014400', '0.014400', '0.014400', '0.014400', '0.014400', '0.014400', '0.014400', '0.014400', '0.014400', '0.014400', '0.013700', '0.013600', '0.013600', '0.013500', '0.013500', '0.013400', '0.013400', '0.013300', '0.013200', '0.013200', '0.013100', '0.013100', '0.013000', '0.012900', '0.012900', '0.012800', '0.012800', '0.012700', '0.012700', '0.012600', '0.012500', '0.012500', '0.012400', '0.012400', '0.012300', '0.010900', '0.010800', '0.010800', '0.010700', '0.010700', '0.010600', '0.010600', '0.010500', '0.010500', '0.010400', '0.010400', '0.010400', '0.010300', '0.010300', '0.010200', '0.010200', '0.010100', '0.010100', '0.010100', '0.010000', '0.010000', '0.009900', '0.009900', '0.009900', '0.009800', '0.009800', '0.009700', '0.009700', '0.009700', '0.009600', '0.009600', '0.009600', '0.009500', '0.009500', '0.009400', '0.009400', '0.009400', '0.009300', '0.009300', '0.009300', '0.009200', '0.009200', '0.009100', '0.009100', '0.009100'],
    'header offset': '0',
    'interleave': 'bsq',
    'lines': '400',
    'samples': '220',
    'sensor type': 'Unknown',
    'wavelength': ['0.460000', '0.465000', '0.470000', '0.475000', '0.479000', '0.484000', '0.489000', '0.494000', '0.499000', '0.503000', '0.508000', '0.513000', '0.518000', '0.523000', '0.528000', '0.533000', '0.538000', '0.543000', '0.549000', '0.554000', '0.559000', '0.565000', '0.570000', '0.575000', '0.581000', '0.587000', '0.592000', '0.598000', '0.604000', '0.610000', '0.616000', '0.622000', '0.628000', '0.634000', '0.640000', '0.646000', '0.653000', '0.659000', '0.665000', '0.672000', '0.679000', '0.685000', '0.692000', '0.699000', '0.706000', '0.713000', '0.720000', '0.727000', '0.734000', '0.741000', '0.749000', '0.756000', '0.763000', '0.771000', '0.778000', '0.786000', '0.793000', '0.801000', '0.809000', '0.817000', '0.824000', '0.832000', '0.840000', '0.848000', '0.856000', '0.864000', '0.872000', '0.880000', '0.888000', '0.896000', '0.915000', '0.924000', '0.934000', '0.944000', '0.955000', '0.965000', '0.975000', '0.986000', '0.997000', '1.007000', '1.018000', '1.029000', '1.040000', '1.051000', '1.063000', '1.074000', '1.086000', '1.097000', '1.109000', '1.120000', '1.132000', '1.144000', '1.155000', '1.167000', '1.179000', '1.191000', '1.203000', '1.215000', '1.227000', '1.239000', '1.251000', '1.263000', '1.275000', '1.287000', '1.299000', '1.311000', '1.323000', '1.522000', '1.534000', '1.545000', '1.557000', '1.568000', '1.579000', '1.590000', '1.601000', '1.612000', '1.624000', '1.634000', '1.645000', '1.656000', '1.667000', '1.678000', '1.689000', '1.699000', '1.710000', '1.721000', '1.731000', '1.742000', '1.752000', '1.763000', '1.773000', '1.783000', '2.044000', '2.053000', '2.062000', '2.071000', '2.080000', '2.089000', '2.098000', '2.107000', '2.115000', '2.124000', '2.133000', '2.141000', '2.150000', '2.159000', '2.167000', '2.176000', '2.184000', '2.193000', '2.201000', '2.210000', '2.218000', '2.226000', '2.234000', '2.243000', '2.251000', '2.259000', '2.267000', '2.275000', '2.283000', '2.292000', '2.300000', '2.308000', '2.315000', '2.323000', '2.331000', '2.339000', '2.347000', '2.355000', '2.363000', '2.370000', '2.378000', '2.386000', '2.393000', '2.401000', '2.409000'],
    'wavelength units': 'Micrometers',
    'y start': '24',
    'z plot titles': ['wavelength [!7l!3m]!N', 'reflectance [*10000]']
    }
}

ENVI domain:
{
'bands': '177',
...
'z plot titles': ['wavelength [!7l!3m]!N', 'reflectance [*10000]']
}

ENVI wavelength:
[0.46, 0.465, 0.47, 0.475, 0.479, 0.484, 0.489, 0.494, 0.499, 0.503, 0.508, 0.513, 0.518, 0.523, 0.528, 0.533, 0.538, 0.543, 0.549, 0.554, 0.559, 0.565, 0.57, 0.575, 0.581, 0.587, 0.592, 0.598, 0.604, 0.61, 0.616, 0.622, 0.628, 0.634, 0.64, 0.646, 0.653, 0.659, 0.665, 0.672, 0.679, 0.685, 0.692, 0.699, 0.706, 0.713, 0.72, 0.727, 0.734, 0.741, 0.749, 0.756, 0.763, 0.771, 0.778, 0.786, 0.793, 0.801, 0.809, 0.817, 0.824, 0.832, 0.84, 0.848, 0.856, 0.864, 0.872, 0.88, 0.888, 0.896, 0.915, 0.924, 0.934, 0.944, 0.955, 0.965, 0.975, 0.986, 0.997, 1.007, 1.018, 1.029, 1.04, 1.051, 1.063, 1.074, 1.086, 1.097, 1.109, 1.12, 1.132, 1.144, 1.155, 1.167, 1.179, 1.191, 1.203, 1.215, 1.227, 1.239, 1.251, 1.263, 1.275, 1.287, 1.299, 1.311, 1.323, 1.522, 1.534, 1.545, 1.557, 1.568, 1.579, 1.59, 1.601, 1.612, 1.624, 1.634, 1.645, 1.656, 1.667, 1.678, 1.689, 1.699, 1.71, 1.721, 1.731, 1.742, 1.752, 1.763, 1.773, 1.783, 2.044, 2.053, 2.062, 2.071, 2.08, 2.089, 2.098, 2.107, 2.115, 2.124, 2.133, 2.141, 2.15, 2.159, 2.167, 2.176, 2.184, 2.193, 2.201, 2.21, 2.218, 2.226, 2.234, 2.243, 2.251, 2.259, 2.267, 2.275, 2.283, 2.292, 2.3, 2.308, 2.315, 2.323, 2.331, 2.339, 2.347, 2.355, 2.363, 2.37, 2.378, 2.386, 2.393, 2.401, 2.409]

Set raster metadata

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)

# set multiple domains
copy = rasterDataset.translate()
copy.setMetadataDict(metadataDict={'domain1': {'a': 1, 'b': 2},
                                   'domain2': {'c': 3, 'd': 4}})
print({domain: copy.metadataDict()[domain] for domain in ['domain1', 'domain2']})

# set domain
copy = rasterDataset.translate()
copy.setMetadataDomain(metadataDomain={'a': 1, 'b': 2}, domain='domain')
print(copy.metadataDict()['domain'])

# set item
copy = rasterDataset.translate()
copy.setMetadataItem(key='a', value=1, domain='domain')
print(copy.metadataDict()['domain'])

Prints:

{'domain1': {'a': '1', 'b': '2'}, 'domain2': {'c': '3', 'd': '4'}}
{'a': '1', 'b': '2'}
{'a': '1'}

Get and set no data value

Warning

todo

Get raster band

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
rasterBandDataset = rasterDataset.band(index=0)
print(rasterBandDataset)

Prints:

RasterBandDataset(raster=RasterDataset(gdalDataset=<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x000001A6FFF34B40> >), index=0)

Read raster data

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)

# read dataset

# - all bands as 3d array
print(rasterDataset.readAsArray().shape)
# - single band as 2d array
print(rasterDataset.band(index=0).readAsArray().shape)
# - pixel (z) profile as 1d array
print(rasterDataset.zprofile(pixel=Pixel(x=100, y=100)).shape)
# - column (x) profile as 1d array
print(rasterDataset.xprofile(row=Row(y=100, z=0)).shape)
# - row (y) profile as 1d array
print(rasterDataset.yprofile(column=Column(x=100, z=0)).shape)

# read dataset for given target grid (may include on-the-fly resampling and/or reprojection)

grid = Grid(extent=Extent(xmin=13, xmax=13.5, ymin=52, ymax=52.5, projection=Projection.wgs84()), resolution=0.001)
# - all bands as 3d array
print(rasterDataset.array(grid=grid, resampleAlg=gdal.GRA_Cubic).shape)
# - single band as 2d array
print(rasterDataset.band(index=0).array(grid=grid, resampleAlg=gdal.GRA_Cubic).shape)

Prints:

(177, 400, 220)
(400, 220)
(177,)
(220,)
(400,)
(177, 500, 500)
(500, 500)

Write raster data

Warning

todo

Loop through all raster bands

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
for band in rasterDataset.bands():
    array = band.readAsArray()
    noDataValue = band.noDataValue()
    values = array[array != noDataValue]
    min = values.min()
    max = values.max()
    mean = values.std()
    print('Band {} Stats: Minimum={}, Maximum={}, Mean={}'.format(band.index()+1, min, max, mean))

Prints:

Band 1 Stats: Minimum=110, Maximum=3759, Mean=212.88073662811726
Band 2 Stats: Minimum=105, Maximum=3757, Mean=214.96880571233484
Band 3 Stats: Minimum=96, Maximum=3792, Mean=217.37383030398556
...

Convert a vector to a raster

import enmapboxtestdata
from hubdc.core import *

vectorDataset = openVectorDataset(filename=enmapboxtestdata.landcover_polygons)
grid = vectorDataset.grid(resolution=5)
rasterDataset = vectorDataset.rasterize(grid, noDataValue=-9999, initValue=-9999,
                                        burnAttribute='level_3_id',
                                        filename='raster.tif', driver=GTiffDriver())
../../../../_images/vector_to_raster.png

Clip a raster with a vector

Clip a raster with the extent from a vector.

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
vectorDataset = openVectorDataset(filename=enmapboxtestdata.landcover_polygons)
grid = rasterDataset.grid().clip(extent=vectorDataset.extent())
clipped = rasterDataset.translate(grid=grid, filename='raster.tif', driver=GTiffDriver())

Note that the result raster grid is snapped to the original raster grid to prevent subpixel shifts. Because of this, some vector geometries may slightly lap over the grid borders.

../../../../_images/clip_raster_with_vector.png

Calculate zonal statistics

Calculates statistics on values (i.e. mean value) of a raster band within the zones given by a vector attribute.

In this example we use the level_3_id attribute as zones.

import enmapboxtestdata
from hubdc.core import *
import numpy as np

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
vectorDataset = openVectorDataset(filename=enmapboxtestdata.landcover_polygons)

# rasterize zones attribute
zonesRasterDataset = vectorDataset.rasterize(grid=rasterDataset.grid(),
                                             burnAttribute='level_3_id')

# calculate zonal statistics for the first raster band
values = rasterDataset.band(index=0).readAsArray()
zones = zonesRasterDataset.band(index=0).readAsArray()

for i in np.unique(zones):
    mean = np.mean(values[zones == i]) # subset pixels for current zone and calculate mean
    print('Zone {} mean: {}'.format(i+1, mean))

Prints:

Zone 1.0 mean: 331.07156981342763
Zone 2.0 mean: 546.2222222222222
Zone 3.0 mean: 546.7610921501706
Zone 4.0 mean: 387.565625
Zone 5.0 mean: 365.40816326530614
Zone 6.0 mean: 838.8627450980392
Zone 7.0 mean: 191.05405405405406

Create raster from array

from hubdc.core import *
import numpy as np

array = np.array([[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1],
                  [ 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1],
                  [ 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1],
                  [ 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1],
                  [ 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])

rasterDataset = RasterDataset.fromArray(array=array, filename='raster.tif', driver=GTiffDriver())
../../../../_images/raster_from_array.png

Create memory raster

Copy a raster into memory.

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(enmapboxtestdata.enmap)

# option a) with MemDriver
driver = MemDriver()
copy = rasterDataset.translate(driver=driver)

# option b) with in-memory files and any driver
driver = GTiffDriver()
options = [driver.Option.INTERLEAVE.BAND, driver.Option.COMPRESS.LZW]
copy = rasterDataset.translate(filename='/vsimem/raster.tif', driver=driver, options=options)

Replace no data value of raster with new value

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)

# read data
array = rasterDataset.readAsArray()

# replace no data value
oldNoDataValue = rasterDataset.noDataValue()
newNoDataValue = -1
array[array==oldNoDataValue] = newNoDataValue

# write data
newRasterDataset = RasterDataset.fromArray(array=array, grid=rasterDataset.grid(),
                                           filename='raster.tif', driver=GTiffDriver())
newRasterDataset.setNoDataValue(value=newNoDataValue)

Raster band dataset

Get raster band information

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.createClassification(30))
rasterBand = rasterDataset.band(index=0)
print('no data value: {}'.format(rasterBand.noDataValue()))
print('description: {}'.format(rasterBand.description()))
print('category names: {}'.format(rasterBand.categoryNames()))
print('category colors: {}'.format(rasterBand.categoryColors()))

Prints:

no data value: 0.0
description: Classification
category names: ['unclassified', 'impervious', 'low vegetation', 'tree', 'soil', 'water']
category colors: [(0, 0, 0, 255), (230, 0, 0, 255), (152, 230, 0, 255), (38, 115, 0, 255), (168, 112, 0, 255), (0, 100, 255, 255)]

Set raster band category names and colors

from hubdc.core import *

rasterDataset = RasterDataset.fromArray(array=[[[0, 1, 2, 3]]], filename='raster.bsq', driver=EnviDriver())
band = rasterDataset.band(index=0)
band.setCategoryNames(names=['unclassified', 'class 1', 'class 2', 'class 3'])
band.setCategoryColors(colors=[(0,0,0), (255, 0, 0), (0, 255, 0), (0, 0, 255)]) # list of rgb or rgba tuples
print(band.categoryNames())
print(band.categoryColors())

Prints:

['unclassified', 'class 1', 'class 2', 'class 3']
[(0, 0, 0, 255), (255, 0, 0, 255), (0, 255, 0, 255), (0, 0, 255, 255)]
..... Classes
class 1
class 2
class 3
../../../../_images/set_raster_band_categories.png

VectorData

Find API reference here: hubdc.core.VectorDataset

Warning

todo

Grid

Create grid

from hubflow.core import *

# create a grid in WGS84 for the whole world with 1 degree resolution
extent=Extent(xmin=-180.0, ymin=-90.0, xmax=180.0, ymax=90.0, projection=Projection.wgs84())
grid = Grid(extent=extent, resolution=1.0)
print(grid)

Prints:

Grid(extent=Extent(xmin=-180.0, xmax=180.0, ymin=-90.0, ymax=90.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])

Get grid information

print(grid.shape())
print(grid.size())
print(grid.geoTransform())

Prints:

(180, 360)
RasterSize(x=360, y=180)
(-180.0, 1.0, 0.0, 90.0, 0.0, -1.0)

Subset grid by pixel offset and size

subgrid = grid.subset(offset=(0, 0), size=(10, 10))
print(subgrid)
print(subgrid.size())

Prints:

Grid(extent=Extent(xmin=-180.0, xmax=-170.0, ymin=80.0, ymax=90.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])
RasterSize(x=10, y=10)

Create systematic subgrids (tiling scheme)

subgrids = grid.subgrids(size=(90, 90))
for subgrid, i, iy, ix in subgrids:
    print('x{}, y{}: {}'.format(ix, iy, subgrid))

Prints:

x0, y0: Grid(extent=Extent(xmin=-180.0, xmax=-90.0, ymin=0.0, ymax=90.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])
x1, y0: Grid(extent=Extent(xmin=-90.0, xmax=0.0, ymin=0.0, ymax=90.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])
x2, y0: Grid(extent=Extent(xmin=0.0, xmax=90.0, ymin=0.0, ymax=90.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])
x3, y0: Grid(extent=Extent(xmin=90.0, xmax=180.0, ymin=0.0, ymax=90.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])
x0, y1: Grid(extent=Extent(xmin=-180.0, xmax=-90.0, ymin=-90.0, ymax=0.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])
x1, y1: Grid(extent=Extent(xmin=-90.0, xmax=0.0, ymin=-90.0, ymax=0.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])
x2, y1: Grid(extent=Extent(xmin=0.0, xmax=90.0, ymin=-90.0, ymax=0.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])
x3, y1: Grid(extent=Extent(xmin=90.0, xmax=180.0, ymin=-90.0, ymax=0.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=1.0, y=1.0), projection=Projection(wkt=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"]])

Anchor a grid to a point

from hubflow.core import *
import enmapboxtestdata

# Get raster grid.
rasterDataset = openRasterDataset(enmapboxtestdata.enmap)
print(rasterDataset.grid())

# Reproject that grid into WGS84 with 0.005 degree resolution.
extentWgs84 = rasterDataset.grid().extent().reproject(projection=Projection.wgs84())
gridWgs84 = Grid(extent=extentWgs84, resolution=0.005)
print(gridWgs84)

# Anchor the grid to (0, 0) degrees.
# This ensures that the upper left and lower right point coordinates are multiples of the 0.005 degree resolution

gridWgs84Anchored = gridWgs84.anchor(Point(x=0, y=0, projection=Projection.wgs84()))
print(gridWgs84Anchored)

Prints:

Get grid coordinates

from hubflow.core import *

extent=Extent(xmin=-180.0, ymin=-90.0, xmax=180.0, ymax=90.0, projection=Projection.wgs84())
grid = Grid(extent=extent, resolution=10)
print(grid)
print(grid.size())

# Get x center coordinats.
print(grid.xMapCoordinates())

# Get y center coordinats.
print(grid.yMapCoordinates())

Prints:

Grid(extent=Extent(xmin=-180.0, xmax=180.0, ymin=-90.0, ymax=90.0, projection=Projection(wkt=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"]])), resolution=Resolution(x=10.0, y=10.0), projection=Projection(wkt=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"]])
RasterSize(x=36, y=18)
[-175.0, -165.0, -155.0, -145.0, -135.0, -125.0, -115.0, -105.0, -95.0, -85.0, -75.0, -65.0, -55.0, -45.0, -35.0, -25.0, -15.0, -5.0, 5.0, 15.0, 25.0, 35.0, 45.0, 55.0, 65.0, 75.0, 85.0, 95.0, 105.0, 115.0, 125.0, 135.0, 145.0, 155.0, 165.0, 175.0]
[85.0, 75.0, 65.0, 55.0, 45.0, 35.0, 25.0, 15.0, 5.0, -5.0, -15.0, -25.0, -35.0, -45.0, -55.0, -65.0, -75.0, -85.0]

Extent

Create extent

from hubflow.core import *

# create a extent in WGS84 for the whole world
extent=Extent(xmin=-180.0, ymin=-90.0, xmax=180.0, ymax=90.0, projection=Projection.wgs84())
print(extent)

Prints:

Extent(xmin=-180.0, xmax=180.0, ymin=-90.0, ymax=90.0, projection=Projection(wkt=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"]]))

Get extent information

Warning

todo

Geometric calculation

Warning

todo

RasterDriver

Create raster driver

from hubdc.core import *

# driver by name
print(RasterDriver(name='GTiff'))

# some predined driver
print(GTiffDriver())
print(EnviDriver())
print(ErdasDriver())

# by file extension
print(RasterDriver.fromFilename('raster.vrt'))

Prints:

RasterDriver(name='GTiff')
GTiffDriver()
EnviDriver()
ErdasDriver()
VrtDriver()

Specify GeoTiff creation options

from hubdc.core import *

driver = GTiffDriver()

# GTiff default options
print(driver.options())

# options for LZW compressed GTiff
print(driver.options() + [driver.Option.COMPRESS.LZW])

# options for JPEG compressed GTiff
print(driver.options() + [driver.Option.COMPRESS.JPEG, driver.Option.JPEG_QUALITY(75)])

# options for tiled GTiff
print(driver.options() + [driver.Option.TILED.YES, driver.Option.BLOCKXSIZE(256), driver.Option.BLOCKYSIZE(256)])

Prints:

['INTERLEAVE=BAND']
['INTERLEAVE=BAND', 'COMPRESS=LZW']
['INTERLEAVE=BAND', 'COMPRESS=JPEG', 'JPEG_QUALITY=75']
['INTERLEAVE=BAND', 'TILED=YES', 'BLOCKXSIZE=256', 'BLOCKYSIZE=256']

Projection

Create projection

from hubdc.core import *

# from well known text
print(Projection(wkt='GEOGCS["WGS84", DATUM["WGS_1984", SPHEROID["WGS84",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"]]'))

# from EPSG code
print(Projection.fromEpsg(epsg=4326))

# some predefined projections
print(Projection.wgs84())
print(Projection.wgs84WebMercator())
print(Projection.utm(zone=33, north=True))
print(Projection.utm(zone=33, north=False))

Prints:

Projection(wkt=GEOGCS["WGS84", DATUM["WGS_1984", SPHEROID["WGS84",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(wkt=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(wkt=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(wkt=PROJCS["WGS 84 / Pseudo-Mercator",     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["Mercator_1SP"],     PARAMETER["central_meridian",0],     PARAMETER["scale_factor",1],     PARAMETER["false_easting",0],     PARAMETER["false_northing",0],     UNIT["metre",1,         AUTHORITY["EPSG","9001"]],     AXIS["X",EAST],     AXIS["Y",NORTH],     EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],     AUTHORITY["EPSG","3857"]])
Projection(wkt=PROJCS["WGS 84 / UTM zone 33N",     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",15],     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","32633"]])
Projection(wkt=PROJCS["WGS 84 / UTM zone 33S",     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",15],     PARAMETER["scale_factor",0.9996],     PARAMETER["false_easting",500000],     PARAMETER["false_northing",10000000],     UNIT["metre",1,         AUTHORITY["EPSG","9001"]],     AXIS["Easting",EAST],     AXIS["Northing",NORTH],     AUTHORITY["EPSG","32733"]])

Reproject point, geometry or extent

from hubdc.core import *

utm = Projection.utm(zone=33)
wgs84 = Projection.wgs84()

point = Point(x=380000, y=5800000, projection=utm)
extent = Extent(xmin=380000, xmax=390000, ymin=5800000, ymax=5830000, projection=utm)
geometry = Geometry(wkt='POLYGON((380000 5830000,390000 5830000,390000 5800000,380000 5800000,380000 5830000))', projection=utm)

print(point.reproject(projection=wgs84))
print(extent.reproject(projection=wgs84))
print(geometry.reproject(projection=wgs84))

Prints:

from hubdc.docutils import createDocPrint

print = createDocPrint(__file__)

# START
from hubdc.core import *

utm = Projection.utm(zone=33)
wgs84 = Projection.wgs84()

point = Point(x=380000, y=5800000, projection=utm)
extent = Extent(xmin=380000, xmax=390000, ymin=5800000, ymax=5830000, projection=utm)
geometry = Geometry(wkt='POLYGON((380000 5830000,390000 5830000,390000 5800000,380000 5800000,380000 5830000))', projection=utm)

print(point.reproject(projection=wgs84))
print(extent.reproject(projection=wgs84))
print(geometry.reproject(projection=wgs84))
# END

Get projection

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
vectorDataset = openVectorDataset(filename=enmapboxtestdata.landcover_polygons)
extent = rasterDataset.extent()
point = extent.upperLeft()
geometry = extent.geometry()

# from raster dataset
projection = rasterDataset.projection()

# from vector dataset
projection = vectorDataset.projection()

# from extent
projection = extent.projection()

# from point
projection = point.projection()

# from geometry
projection = geometry.projection()

Reproject a vector dataset

import enmapboxtestdata
from hubdc.core import *

vectorDataset = openVectorDataset(filename=enmapboxtestdata.landcover_polygons)
vectorDataset.reproject(projection=Projection.wgs84(), filename='vector.gpkg', driver=GeoPackageDriver())

Reproject a raster dataset

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)

# define target grid with extent of original raster, but in WGS 84 projection and 0.0001 degree resolution
grid = Grid(extent=rasterDataset.extent().reproject(Projection.wgs84()), resolution=0.0001)

# reproject raster into target grid
warped = rasterDataset.warp(grid=grid, filename='raster.bsq', driver=EnviDriver())

Export projection

from hubdc.core import *

projection = Projection.wgs84()

# to well known text
print(projection.wkt())

# to OSR spatial reference
print(type(projection.osrSpatialReference()))

Prints:

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"]]
<class 'osgeo.osr.SpatialReference'>

MapViewer

The map viewer is an interactive application for exploring maps.

View single map layer

Shortcuts for initialising a map viewer directly from a raster or vector dataset.

Multi band raster

import enmapboxtestdata
from hubdc.core import *

# Starting a map viewer from a multiband raster results in a MultiBandColorRenderer representation.
# Note that enmapboxtestdata.enmap has a default style defined (see .qml file next to it).
rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
MapViewer().addLayer(rasterDataset.mapLayer()).show()

../../../../_images/view_single_map1.png

Layer with MultiBandColorRenderer

Single band raster

# Starting a map viewer from a raster band results in a SingleBandGrayRenderer representation.
rasterBandDataset = rasterDataset.band(index=0)
MapViewer().addLayer(rasterBandDataset.mapLayer()).show()

../../../../_images/view_single_map2.png

Layer with SingleBandGrayRenderer

Vector

# Starting a map viewer from a vector results in a SingleSymbolRenderer representation by default.
# Note that enmapboxtestdata.landcover_points uses a CategorizedSymbolRenderer as default (see .qml file next to it).
vectorDataset = openVectorDataset(filename=enmapboxtestdata.landcover_points)
#vectorDataset.mapLayer().show()
#vectorDataset.mapViewer().show()
MapViewer().addLayer(vectorDataset.mapLayer()).show()

../../../../_images/view_single_map3.png

Layer with …

Paletted raster

# Starting a map viewer from a raster band with a color lookup table results in a PalettedRasterRenderer representation.
rasterDatasetWithLookupTable = openRasterDataset(enmapboxtestdata.createClassification(gridOrResolution=10, level='level_3_id'))
MapViewer().addLayer(rasterDatasetWithLookupTable.mapLayer()).show()
../../../../_images/view_single_map4.png

Layer with PalettedRasterRenderer

Add multiple map layer

import enmapboxtestdata
from hubdc.core import *

vectorDataset = openVectorDataset(filename=enmapboxtestdata.landcover_points)
rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)

mapViewer = MapViewer()
mapViewer.addLayer(vectorDataset.mapLayer())
mapViewer.addLayer(rasterDataset.mapLayer())
mapViewer.show()
../../../../_images/view_multi_map.png

MapViewer with multiple map layer

Set viewer extent and projection

import enmapboxtestdata
from hubdc.core import *

rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)

mapViewer = rasterDataset.mapViewer()

# change projection
mapViewer.setProjection(projection=Projection.wgs84())
mapViewer.show()
../../../../_images/set_extent1.png
# change extent
mapViewer.setExtent(Extent(xmin=13.29, xmax=13.32, ymin=52.47, ymax=52.49, projection=Projection.wgs84()))
mapViewer.show()
../../../../_images/set_extent2.png

Rendering configuration

MultiBandColorRenderer

import enmapboxtestdata
from hubdc.core import *

# setup map viewer
rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
mapViewer = MapViewer()
mapLayer = rasterDataset.mapLayer()
mapViewer.addLayer(mapLayer)

# setup layer rendering

# - default is to apply a 2% - 98% stretch on the selected bands
mapLayer.initMultiBandColorRenderer(redIndex=63, greenIndex=37, blueIndex=23)

# - but it's also possible to set min-max stretch values explicitely
mapLayer.initMultiBandColorRenderer(redIndex=63, greenIndex=37, blueIndex=23,
                                    redMin=809, redMax=2861,
                                    greenMin=183, greenMax=1370,
                                    blueMin=279, blueMax=1275)

mapViewer.show()
../../../../_images/rendering_MultiBandColorRenderer.png

SingleBandGrayRenderer

import enmapboxtestdata
from hubdc.core import *
from qgis.core import QgsSingleBandGrayRenderer, QgsContrastEnhancement

# setup map viewer
rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)
rasterBandDataset = rasterDataset.band(index=0)
mapViewer = MapViewer()
mapLayer = rasterDataset.band(index=0).mapLayer()
mapViewer.addLayer(mapLayer)

# setup layer rendering

# - default is to apply a 2% - 98% stretch on the selected band
mapLayer.initSingleBandGrayRenderer(grayIndex=0)

# - but it's also possible to set min-max stretch values explicitely
mapLayer.initSingleBandGrayRenderer(grayIndex=0, grayMin=179, grayMax=1026)

mapViewer.show()
../../../../_images/rendering_QgsSingleBandGrayRenderer.png

PalettedRasterRenderer

Warning

todo

SingleSymbolRenderer

Warning

todo

CategorizedSymbolRenderer

Warning

todo

Save map viewer content to PNG file

import enmapboxtestdata
from hubdc.core import *

vectorDataset = openVectorDataset(filename=enmapboxtestdata.landcover_points)
rasterDataset = openRasterDataset(filename=enmapboxtestdata.enmap)

mapViewer = MapViewer()
mapViewer.addLayer(vectorDataset.mapLayer())
mapViewer.addLayer(rasterDataset.mapLayer())
mapViewer.save('image.png')