14332609845

Committed 08 Apr 2025 11:41AM UTC coverage: 80.084% (+0.2%) from 79.933%

Build # 14332609845

Build Type

Pull #55

github

Committed by

johntruckenbrodt

Commit Message

[docs] add geopandas intersphinx mapping

Pull Request Pull Request #55: [Vector.to_geopandas] new method

Run Details

16 of 19 new or added lines in 2 files covered. (84.21%)

63 existing lines in 2 files now uncovered.

1914 of 2390 relevant lines covered (80.08%)

0.8 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

99.37

/spatialist/tests/test_spatial.py

import os
import shutil
import pytest
import platform
import numpy as np
from osgeo import ogr, gdal
from spatialist.raster import Dtype, png, Raster, stack, rasterize
from spatialist.vector import feature2vector, dissolve, Vector, intersect, bbox, wkt2vector
from spatialist.envi import hdr, HDRobject
from spatialist.sqlite_util import sqlite_setup, __Handler
from spatialist.ancillary import parallel_apply_along_axis
from spatialist.auxil import (crsConvert, haversine, ogr2ogr, gdal_translate, gdal_rasterize, gdalwarp,
                              utm_autodetect, coordinate_reproject, cmap_mpl2gdal)

import logging

logging.basicConfig(level=logging.DEBUG)


def test_crsConvert():
    assert crsConvert(crsConvert(4326, 'wkt'), 'proj4').strip() == '+proj=longlat +datum=WGS84 +no_defs'
    assert crsConvert(crsConvert(4326, 'prettyWkt'), 'opengis') == 'https://www.opengis.net/def/crs/EPSG/0/4326'
    assert crsConvert('https://www.opengis.net/def/crs/EPSG/0/4326', 'epsg') == 4326
    assert crsConvert(crsConvert('https://www.opengis.net/def/crs/EPSG/0/4326', 'osr'), 'epsg') == 4326
    assert crsConvert('EPSG:4326+5773', 'proj4').strip() \
           == '+proj=longlat +datum=WGS84 +geoidgrids=egm96_15.gtx +vunits=m +no_defs' \
           or '+proj=longlat +datum=WGS84 +vunits=m +no_defs'
    with pytest.raises(TypeError):
        crsConvert('xyz', 'epsg')
    with pytest.raises(ValueError):
        crsConvert(4326, 'xyz')


def test_haversine():
    assert haversine(50, 10, 51, 10) == 111194.92664455889


def test_Vector(tmpdir, testdata):
    scene = Raster(testdata['tif'])
    bbox1 = scene.bbox()
    assert bbox1.getArea() == 23262400.0
    assert bbox1.extent == {'ymax': 4830114.70107, 'ymin': 4825774.70107, 'xmin': 620048.241204, 'xmax': 625408.241204}
    assert bbox1.nlayers == 1
    assert bbox1.getProjection('epsg') == 32631
    assert bbox1.proj4.strip() == '+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs'
    assert isinstance(bbox1.getFeatureByIndex(0), ogr.Feature)
    with pytest.raises(IndexError):
        bbox1.getFeatureByIndex(1)
    bbox1.reproject(4326)
    assert bbox1.proj4.strip() == '+proj=longlat +datum=WGS84 +no_defs'
    ext = {key: round(val, 3) for key, val in bbox1.extent.items()}
    assert ext == {'xmax': 4.554, 'xmin': 4.487, 'ymax': 43.614, 'ymin': 43.574}
    assert utm_autodetect(bbox1, 'epsg') == 32631
    assert isinstance(bbox1['fid=0'], Vector)
    with pytest.raises(RuntimeError):
        test = bbox1[0.1]
    assert bbox1.fieldnames == ['area']
    assert bbox1.getUniqueAttributes('area') == [23262400.0]
    feat = bbox1.getFeatureByAttribute('area', 23262400.0)
    assert isinstance(feat, ogr.Feature)
    bbox2 = feature2vector(feat, ref=bbox1)
    bbox2.close()
    feat.Destroy()
    with pytest.raises(KeyError):
        select = bbox1.getFeatureByAttribute('foo', 'bar')
    with pytest.raises(OSError):
        vec = Vector(filename='foobar')
    gdf_out = tmpdir / "test.gpkg"
    gdf = bbox1.to_geopandas()
    gdf.to_file(str(gdf_out))
    assert gdf_out.exists()
    bbox1.close()


def test_dissolve(tmpdir, travis, testdata):
    scene = Raster(testdata['tif'])
    bbox1 = scene.bbox()
    # retrieve extent and shift its coordinates by one unit
    ext = bbox1.extent
    for key in ext.keys():
        ext[key] += 1
    # create new bbox shapefile with modified extent
    bbox2_name = os.path.join(str(tmpdir), 'bbox2.shp')
    bbox(ext, bbox1.srs, bbox2_name)
    # assert intersection between the two bboxes and combine them into one
    with Vector(bbox2_name) as bbox2:
        assert intersect(bbox1, bbox2) is not None
        bbox1.addvector(bbox2)
        # write combined bbox into new shapefile
        bbox3_name = os.path.join(str(tmpdir), 'bbox3.shp')
        bbox1.write(bbox3_name)
    bbox1.close()
    
    if not travis and platform.system() != 'Windows':
        # dissolve the geometries in bbox3 and write the result to new bbox4
        # this test is currently disabled for Travis as the current sqlite3 version on Travis seems to not support
        # loading gdal as extension; Travis CI setup: Ubuntu 14.04 (Trusty), sqlite3 version 3.8.2 (2018-06-04)
        bbox4_name = os.path.join(str(tmpdir), 'bbox4.shp')
        dissolve(bbox3_name, bbox4_name, field='area')
        assert os.path.isfile(bbox4_name)


def test_Raster(tmpdir, testdata):
    with pytest.raises(RuntimeError):
        with Raster(1) as ras:
            print(ras)
    with Raster(testdata['tif']) as ras:
        print(ras)
        assert ras.bands == 1
        assert ras.proj4.strip() == '+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs'
        assert ras.cols == 268
        assert ras.rows == 217
        assert ras.dim == (217, 268, 1)
        assert ras.dtype == 'Float32'
        assert ras.epsg == 32631
        assert ras.format == 'GTiff'
        assert ras.geo == {'ymax': 4830114.70107, 'rotation_y': 0.0, 'rotation_x': 0.0, 'xmax': 625408.241204,
                           'xres': 20.0, 'xmin': 620048.241204, 'ymin': 4825774.70107, 'yres': -20.0}
        assert ras.geogcs == 'WGS 84'
        assert ras.is_valid() is True
        assert ras.proj4args == {'units': 'm', 'no_defs': None, 'datum': 'WGS84', 'proj': 'utm', 'zone': '31'}
        assert ras.allstats() == [{'min': -26.65471076965332, 'max': 1.4325850009918213,
                                   'mean': -12.124929534450377, 'sdev': 4.738273594738293}]
        assert ras.bbox().getArea() == 23262400.0
        assert len(ras.layers()) == 1
        assert ras.projcs == 'WGS 84 / UTM zone 31N'
        assert ras.res == (20.0, 20.0)
        
        # test writing a subset with no original data in memory
        outname = os.path.join(str(tmpdir), 'test_sub.tif')
        with ras[0:200, 0:100] as sub:
            sub.write(outname, format='GTiff')
        with Raster(outname) as ras2:
            assert ras2.cols == 100
            assert ras2.rows == 200
        
        ras.load()
        mat = ras.matrix()
        assert isinstance(mat, np.ndarray)
        ras.assign(mat, band=0)
        # ras.reduce()
        ras.rescale(lambda x: 10 * x)
        
        # test writing data with original data in memory
        ras.write(os.path.join(str(tmpdir), 'test'), format='GTiff')
        with pytest.raises(RuntimeError):
            ras.write(os.path.join(str(tmpdir), 'test.tif'), format='GTiff')


def test_Raster_subset(testdata):
    with Raster(testdata['tif']) as ras:
        ext = ras.bbox().extent
        xres, yres = ras.res
        ext['xmin'] += xres
        ext['xmax'] -= xres
        ext['ymin'] += yres
        ext['ymax'] -= yres
        with bbox(ext, ras.projection) as vec:
            with ras[vec] as sub:
                xres, yres = ras.res
                assert sub.geo['xmin'] - ras.geo['xmin'] == xres
                assert ras.geo['xmax'] - sub.geo['xmax'] == xres
                assert sub.geo['ymin'] - ras.geo['ymin'] == xres
                assert ras.geo['ymax'] - sub.geo['ymax'] == xres


def test_Raster_extract(testdata):
    with Raster(testdata['tif']) as ras:
        assert ras.extract(px=624000, py=4830000, radius=5) == -10.48837461270875
        with pytest.raises(RuntimeError):
            ras.extract(1, 4830000)
        with pytest.raises(RuntimeError):
            ras.extract(624000, 1)
        
        # ensure corner extraction capability
        assert ras.extract(px=ras.geo['xmin'], py=ras.geo['ymax']) == -10.147890090942383
        assert ras.extract(px=ras.geo['xmin'], py=ras.geo['ymin']) == -14.640368461608887
        assert ras.extract(px=ras.geo['xmax'], py=ras.geo['ymax']) == -9.599242210388182
        assert ras.extract(px=ras.geo['xmax'], py=ras.geo['ymin']) == -9.406558990478516
        
        # test nodata handling capability and correct indexing
        mat = ras.matrix()
        mat[0:10, 0:10] = ras.nodata
        mat[207:217, 258:268] = ras.nodata
        ras.assign(mat, band=0)
        assert ras.extract(px=ras.geo['xmin'], py=ras.geo['ymax'], radius=5) == ras.nodata
        assert ras.extract(px=ras.geo['xmax'], py=ras.geo['ymin'], radius=5) == ras.nodata


def test_Raster_filestack(testdata):
    with pytest.raises(RuntimeError):
        with Raster([testdata['tif']]) as ras:
            print(ras)
    with Raster([testdata['tif'], testdata['tif2']]) as ras:
        assert ras.bands == 2
        arr = ras.array()
    mean = parallel_apply_along_axis(np.nanmean, axis=2, arr=arr, cores=4)
    assert mean.shape == (217, 268)


def test_dtypes():
    assert Dtype('Float32').gdalint == 6
    assert Dtype(6).gdalstr == 'Float32'
    assert Dtype('uint32').gdalstr == 'UInt32'
    with pytest.raises(ValueError):
        Dtype('foobar')
    with pytest.raises(ValueError):
        Dtype(999)
    with pytest.raises(TypeError):
        Dtype(None)


def test_stack(tmpdir, testdata):
    name = testdata['tif']
    outname = os.path.join(str(tmpdir), 'test')
    tr = (30, 30)
    # no input files provided
    with pytest.raises(RuntimeError):
        stack(srcfiles=[], resampling='near', targetres=tr,
              srcnodata=-99, dstnodata=-99, dstfile=outname)
    
    # two files, but only one layer name
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name, name], resampling='near', targetres=tr,
              srcnodata=-99, dstnodata=-99, dstfile=outname, layernames=['a'])
    
    # targetres must be a two-entry tuple/list
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name, name], resampling='near', targetres=30,
              srcnodata=-99, dstnodata=-99, dstfile=outname)
    
    # only one file specified
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name], resampling='near', targetres=tr, overwrite=True,
              srcnodata=-99, dstnodata=-99, dstfile=outname)
    
    # targetres must contain two values
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name, name], resampling='near', targetres=(30, 30, 30),
              srcnodata=-99, dstnodata=-99, dstfile=outname)
    
    # unknown resampling method
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name, name], resampling='foobar', targetres=tr,
              srcnodata=-99, dstnodata=-99, dstfile=outname)
    
    # non-existing files
    with pytest.raises(RuntimeError):
        stack(srcfiles=['foo', 'bar'], resampling='near', targetres=tr,
              srcnodata=-99, dstnodata=-99, dstfile=outname)
    
    # create a multi-band stack
    stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=True,
          srcnodata=-99, dstnodata=-99, dstfile=outname, layernames=['test1', 'test2'])
    with Raster(outname) as ras:
        assert ras.bands == 2
        # Raster.rescale currently only supports one band
        with pytest.raises(ValueError):
            ras.rescale(lambda x: x * 10)
    
    # outname exists and overwrite is False
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=False,
              srcnodata=-99, dstnodata=-99, dstfile=outname, layernames=['test1', 'test2'])
    
    # pass shapefile
    outname = os.path.join(str(tmpdir), 'test2')
    with Raster(name).bbox() as box:
        stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=True,
              srcnodata=-99, dstnodata=-99, dstfile=outname, shapefile=box, layernames=['test1', 'test2'])
    with Raster(outname) as ras:
        assert ras.bands == 2
    
    # pass shapefile and do mosaicing
    outname = os.path.join(str(tmpdir), 'test3')
    with Raster(name).bbox() as box:
        stack(srcfiles=[[name, name]], resampling='near', targetres=tr, overwrite=True,
              srcnodata=-99, dstnodata=-99, dstfile=outname, shapefile=box)
    with Raster(outname + '.tif') as ras:
        assert ras.bands == 1
        assert ras.format == 'GTiff'
    
    # projection mismatch
    name2 = os.path.join(str(tmpdir), os.path.basename(name))
    outname = os.path.join(str(tmpdir), 'test4')
    gdalwarp(src=name, dst=name2, dstSRS=crsConvert(4326, 'wkt'))
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name, name2], resampling='near', targetres=tr, overwrite=True,
              srcnodata=-99, dstnodata=-99, dstfile=outname)
    
    # no projection found
    outname = os.path.join(str(tmpdir), 'test5')
    gdal_translate(src=name, dst=name2, options=['-co', 'PROFILE=BASELINE'])
    with Raster(name2) as ras:
        print(ras.projection)
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name2, name2], resampling='near', targetres=tr, overwrite=True,
              srcnodata=-99, dstnodata=-99, dstfile=outname)
    
    # create separate GeoTiffs
    outdir = os.path.join(str(tmpdir), 'subdir')
    stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=True, layernames=['test1', 'test2'],
          srcnodata=-99, dstnodata=-99, dstfile=outdir, separate=True, compress=True)
    
    # repeat with overwrite disabled (no error raised, just a print message)
    stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=False, layernames=['test1', 'test2'],
          srcnodata=-99, dstnodata=-99, dstfile=outdir, separate=True, compress=True)
    
    # repeat without layernames but sortfun
    # bandnames not unique
    outdir = os.path.join(str(tmpdir), 'subdir2')
    with pytest.raises(RuntimeError):
        stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=True, sortfun=os.path.basename,
              srcnodata=-99, dstnodata=-99, dstfile=outdir, separate=True, compress=True)
    
    # repeat without layernames but sortfun
    name2 = os.path.join(str(tmpdir), os.path.basename(name).replace('VV', 'XX'))
    shutil.copyfile(name, name2)
    outdir = os.path.join(str(tmpdir), 'subdir2')
    stack(srcfiles=[name, name2], resampling='near', targetres=tr, overwrite=True, sortfun=os.path.basename,
          srcnodata=-99, dstnodata=-99, dstfile=outdir, separate=True, compress=True)
    
    # shapefile filtering
    outdir = os.path.join(str(tmpdir), 'subdir3')
    files = [testdata['tif'], testdata['tif2'], testdata['tif3']]
    with Raster(files[0]).bbox() as box:
        stack(srcfiles=files, resampling='near', targetres=(30, 30),
              overwrite=False, layernames=['test1', 'test2', 'test3'],
              srcnodata=-99, dstnodata=-99, dstfile=outdir,
              separate=True, compress=True, shapefile=box)
        # repeated run with different scene selection and only one scene after spatial filtering
        stack(srcfiles=files[1:], resampling='near', targetres=(30, 30),
              overwrite=True, layernames=['test2', 'test3'],
              srcnodata=-99, dstnodata=-99, dstfile=outdir,
              separate=True, compress=True, shapefile=box)


def test_auxil(tmpdir, testdata):
    dir = str(tmpdir)
    with Raster(testdata['tif']) as ras:
        bbox = os.path.join(dir, 'bbox.shp')
        ras.bbox(bbox)
        ogr2ogr(src=bbox, dst=os.path.join(dir, 'bbox.gml'), format='GML')
        gdal_translate(src=ras.raster, dst=os.path.join(dir, 'test'), format='ENVI')
    gdal_rasterize(src=bbox, dst=os.path.join(dir, 'test2'), format='GTiff', xRes=20, yRes=20)


def test_auxil_coordinate_reproject():
    point = coordinate_reproject(x=11, y=51, s_crs=4326, t_crs=32632)
    assert round(point[0], 3) == 640333.296
    assert round(point[1], 3) == 5651728.683


def test_auxil_cmap_mpl2gdal():
    cmap = cmap_mpl2gdal(mplcolor='YlGnBu', values=range(0, 100))
    assert type(cmap) == gdal.ColorTable


def test_rasterize(tmpdir, testdata):
    outname = os.path.join(str(tmpdir), 'test.shp')
    with Raster(testdata['tif']) as ras:
        vec = ras.bbox()
        
        # test length mismatch between burn_values and expressions
        with pytest.raises(RuntimeError):
            rasterize(vec, reference=ras, outname=outname, burn_values=[1], expressions=['foo', 'bar'])
        
        # test a faulty expression
        with pytest.raises(RuntimeError):
            rasterize(vec, reference=ras, outname=outname, burn_values=[1], expressions=['foo'])
        
        # test default parametrization
        rasterize(vec, reference=ras, outname=outname)
        assert os.path.isfile(outname)
        
        # test appending to existing file with valid expression
        rasterize(vec, reference=ras, outname=outname, append=True, burn_values=[1], expressions=['area=23262400.0'])
        
        # test wrong input type for reference
        with pytest.raises(RuntimeError):
            rasterize(vec, reference='foobar', outname=outname)


def test_envi(tmpdir):
    with pytest.raises(RuntimeError):
        obj = HDRobject(1)
    with pytest.raises(RuntimeError):
        obj = HDRobject('foobar')
    outname = os.path.join(str(tmpdir), 'test')
    with HDRobject() as header:
        header.band_names = ['band1']
        header.write(outname)
    outname += '.hdr'
    with HDRobject(outname) as header:
        assert header.band_names == ['band1']
        vals = vars(header)
    with HDRobject(vals) as header:
        assert header.byte_order == 0
    hdr(vals, outname + '2')


def test_sqlite():
    with pytest.raises(RuntimeError):
        con = sqlite_setup(extensions='spatialite')
    con = sqlite_setup(extensions=['spatialite'])
    con.close()
    con = __Handler()
    assert sorted(con.version.keys()) == ['sqlite']
    
    con = __Handler(extensions=['spatialite'])
    assert sorted(con.version.keys()) == ['spatialite', 'sqlite']
    assert 'spatial_ref_sys' in con.get_tablenames()


def test_png(tmpdir, testdata):
    outname = os.path.join(str(tmpdir), 'test')
    with Raster(testdata['tif']) as ras:
        png(src=ras, dst=outname, percent=100, scale=(2, 98), worldfile=True)
    assert os.path.isfile(outname + '.png')
    
    with pytest.raises(TypeError):
        png(src=testdata['tif'], dst=outname, percent=100, scale=(2, 98), worldfile=True)
    
    src = [testdata['tif'], testdata['tif2']]
    with pytest.raises(ValueError):
        with Raster(src) as ras:
            png(src=ras, dst=outname, percent=100, scale=(2, 98), worldfile=True)
    
    src.append(testdata['tif3'])
    outname = os.path.join(str(tmpdir), 'test_rgb.png')
    with Raster(src) as ras:
        png(src=ras, dst=outname, percent=100, scale=(2, 98), worldfile=True)


def test_addfield():
    extent = {'xmin': 10, 'xmax': 11, 'ymin': 50, 'ymax': 51}
    with bbox(coordinates=extent, crs=4326) as box:
        box.addfield(name='test1', type=ogr.OFTString, values=['a'])
        box.addfield(name='test2', type=ogr.OFTStringList, values=[['a', 'b']])
        box.addfield(name='test3', type=ogr.OFTInteger, values=[1])
        box.addfield(name='test4', type=ogr.OFTIntegerList, values=[[1, 2]])
        box.addfield(name='test5', type=ogr.OFTInteger64, values=[1])
        box.addfield(name='test6', type=ogr.OFTInteger64List, values=[[1, 2]])
        box.addfield(name='test7', type=ogr.OFTReal, values=[1])
        box.addfield(name='test8', type=ogr.OFTRealList, values=[[1., 2.]])
        box.addfield(name='test9', type=ogr.OFTBinary, values=[b'1'])


def test_wkt2vector():
    wkt1 = 'POLYGON ((0. 0., 0. 1., 1. 1., 1. 0., 0. 0.))'
    wkt2 = 'POLYGON ((1. 1., 1. 2., 2. 2., 2. 1., 1. 1.))'
    with wkt2vector(wkt1, srs=4326) as vec:
        assert vec.getArea() == 1.
    with wkt2vector([wkt1, wkt2], srs=4326) as vec:
        assert vec.getArea() == 2.

1	import os	1✔
2	import shutil	1✔
3	import pytest	1✔
4	import platform	1✔
5	import numpy as np	1✔
6	from osgeo import ogr, gdal	1✔
7	from spatialist.raster import Dtype, png, Raster, stack, rasterize	1✔
8	from spatialist.vector import feature2vector, dissolve, Vector, intersect, bbox, wkt2vector	1✔
9	from spatialist.envi import hdr, HDRobject	1✔
10	from spatialist.sqlite_util import sqlite_setup, __Handler	1✔
11	from spatialist.ancillary import parallel_apply_along_axis	1✔
12	from spatialist.auxil import (crsConvert, haversine, ogr2ogr, gdal_translate, gdal_rasterize, gdalwarp,	1✔
13	utm_autodetect, coordinate_reproject, cmap_mpl2gdal)
14
15	import logging	1✔
16
17	logging.basicConfig(level=logging.DEBUG)	1✔
18
19
20	def test_crsConvert():	1✔
21	assert crsConvert(crsConvert(4326, 'wkt'), 'proj4').strip() == '+proj=longlat +datum=WGS84 +no_defs'	1✔
22	assert crsConvert(crsConvert(4326, 'prettyWkt'), 'opengis') == 'https://www.opengis.net/def/crs/EPSG/0/4326'	1✔
23	assert crsConvert('https://www.opengis.net/def/crs/EPSG/0/4326', 'epsg') == 4326	1✔
24	assert crsConvert(crsConvert('https://www.opengis.net/def/crs/EPSG/0/4326', 'osr'), 'epsg') == 4326	1✔
25	assert crsConvert('EPSG:4326+5773', 'proj4').strip() \	1✔
26	== '+proj=longlat +datum=WGS84 +geoidgrids=egm96_15.gtx +vunits=m +no_defs' \
27	or '+proj=longlat +datum=WGS84 +vunits=m +no_defs'
28	with pytest.raises(TypeError):	1✔
29	crsConvert('xyz', 'epsg')	1✔
30	with pytest.raises(ValueError):	1✔
31	crsConvert(4326, 'xyz')	1✔
32
33
34	def test_haversine():	1✔
35	assert haversine(50, 10, 51, 10) == 111194.92664455889	1✔
36
37
38	def test_Vector(tmpdir, testdata):	1✔
39	scene = Raster(testdata['tif'])	1✔
40	bbox1 = scene.bbox()	1✔
41	assert bbox1.getArea() == 23262400.0	1✔
42	assert bbox1.extent == {'ymax': 4830114.70107, 'ymin': 4825774.70107, 'xmin': 620048.241204, 'xmax': 625408.241204}	1✔
43	assert bbox1.nlayers == 1	1✔
44	assert bbox1.getProjection('epsg') == 32631	1✔
45	assert bbox1.proj4.strip() == '+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs'	1✔
46	assert isinstance(bbox1.getFeatureByIndex(0), ogr.Feature)	1✔
47	with pytest.raises(IndexError):	1✔
48	bbox1.getFeatureByIndex(1)	1✔
49	bbox1.reproject(4326)	1✔
50	assert bbox1.proj4.strip() == '+proj=longlat +datum=WGS84 +no_defs'	1✔
51	ext = {key: round(val, 3) for key, val in bbox1.extent.items()}	1✔
52	assert ext == {'xmax': 4.554, 'xmin': 4.487, 'ymax': 43.614, 'ymin': 43.574}	1✔
53	assert utm_autodetect(bbox1, 'epsg') == 32631	1✔
54	assert isinstance(bbox1['fid=0'], Vector)	1✔
55	with pytest.raises(RuntimeError):	1✔
56	test = bbox1[0.1]	1✔
57	assert bbox1.fieldnames == ['area']	1✔
58	assert bbox1.getUniqueAttributes('area') == [23262400.0]	1✔
59	feat = bbox1.getFeatureByAttribute('area', 23262400.0)	1✔
60	assert isinstance(feat, ogr.Feature)	1✔
61	bbox2 = feature2vector(feat, ref=bbox1)	1✔
62	bbox2.close()	1✔
63	feat.Destroy()	1✔
64	with pytest.raises(KeyError):	1✔
65	select = bbox1.getFeatureByAttribute('foo', 'bar')	1✔
66	with pytest.raises(OSError):	1✔
67	vec = Vector(filename='foobar')	1✔
68	gdf_out = tmpdir / "test.gpkg"	1✔
69	gdf = bbox1.to_geopandas()	1✔
70	gdf.to_file(str(gdf_out))	1✔
71	assert gdf_out.exists()	1✔
72	bbox1.close()	1✔
73
74
75	def test_dissolve(tmpdir, travis, testdata):	1✔
76	scene = Raster(testdata['tif'])	1✔
77	bbox1 = scene.bbox()	1✔
78	# retrieve extent and shift its coordinates by one unit
79	ext = bbox1.extent	1✔
80	for key in ext.keys():	1✔
81	ext[key] += 1	1✔
82	# create new bbox shapefile with modified extent
83	bbox2_name = os.path.join(str(tmpdir), 'bbox2.shp')	1✔
84	bbox(ext, bbox1.srs, bbox2_name)	1✔
85	# assert intersection between the two bboxes and combine them into one
86	with Vector(bbox2_name) as bbox2:	1✔
87	assert intersect(bbox1, bbox2) is not None	1✔
88	bbox1.addvector(bbox2)	1✔
89	# write combined bbox into new shapefile
90	bbox3_name = os.path.join(str(tmpdir), 'bbox3.shp')	1✔
91	bbox1.write(bbox3_name)	1✔
92	bbox1.close()	1✔
93
94	if not travis and platform.system() != 'Windows':	1✔
95	# dissolve the geometries in bbox3 and write the result to new bbox4
96	# this test is currently disabled for Travis as the current sqlite3 version on Travis seems to not support
97	# loading gdal as extension; Travis CI setup: Ubuntu 14.04 (Trusty), sqlite3 version 3.8.2 (2018-06-04)
98	bbox4_name = os.path.join(str(tmpdir), 'bbox4.shp')	1✔
99	dissolve(bbox3_name, bbox4_name, field='area')	1✔
100	assert os.path.isfile(bbox4_name)	1✔
101
102
103	def test_Raster(tmpdir, testdata):	1✔
104	with pytest.raises(RuntimeError):	1✔
105	with Raster(1) as ras:	1✔
UNCOV 106	print(ras)	×
107	with Raster(testdata['tif']) as ras:	1✔
108	print(ras)	1✔
109	assert ras.bands == 1	1✔
110	assert ras.proj4.strip() == '+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs'	1✔
111	assert ras.cols == 268	1✔
112	assert ras.rows == 217	1✔
113	assert ras.dim == (217, 268, 1)	1✔
114	assert ras.dtype == 'Float32'	1✔
115	assert ras.epsg == 32631	1✔
116	assert ras.format == 'GTiff'	1✔
117	assert ras.geo == {'ymax': 4830114.70107, 'rotation_y': 0.0, 'rotation_x': 0.0, 'xmax': 625408.241204,	1✔
118	'xres': 20.0, 'xmin': 620048.241204, 'ymin': 4825774.70107, 'yres': -20.0}
119	assert ras.geogcs == 'WGS 84'	1✔
120	assert ras.is_valid() is True	1✔
121	assert ras.proj4args == {'units': 'm', 'no_defs': None, 'datum': 'WGS84', 'proj': 'utm', 'zone': '31'}	1✔
122	assert ras.allstats() == [{'min': -26.65471076965332, 'max': 1.4325850009918213,	1✔
123	'mean': -12.124929534450377, 'sdev': 4.738273594738293}]
124	assert ras.bbox().getArea() == 23262400.0	1✔
125	assert len(ras.layers()) == 1	1✔
126	assert ras.projcs == 'WGS 84 / UTM zone 31N'	1✔
127	assert ras.res == (20.0, 20.0)	1✔
128
129	# test writing a subset with no original data in memory
130	outname = os.path.join(str(tmpdir), 'test_sub.tif')	1✔
131	with ras[0:200, 0:100] as sub:	1✔
132	sub.write(outname, format='GTiff')	1✔
133	with Raster(outname) as ras2:	1✔
134	assert ras2.cols == 100	1✔
135	assert ras2.rows == 200	1✔
136
137	ras.load()	1✔
138	mat = ras.matrix()	1✔
139	assert isinstance(mat, np.ndarray)	1✔
140	ras.assign(mat, band=0)	1✔
141	# ras.reduce()
142	ras.rescale(lambda x: 10 * x)	1✔
143
144	# test writing data with original data in memory
145	ras.write(os.path.join(str(tmpdir), 'test'), format='GTiff')	1✔
146	with pytest.raises(RuntimeError):	1✔
147	ras.write(os.path.join(str(tmpdir), 'test.tif'), format='GTiff')	1✔
148
149
150	def test_Raster_subset(testdata):	1✔
151	with Raster(testdata['tif']) as ras:	1✔
152	ext = ras.bbox().extent	1✔
153	xres, yres = ras.res	1✔
154	ext['xmin'] += xres	1✔
155	ext['xmax'] -= xres	1✔
156	ext['ymin'] += yres	1✔
157	ext['ymax'] -= yres	1✔
158	with bbox(ext, ras.projection) as vec:	1✔
159	with ras[vec] as sub:	1✔
160	xres, yres = ras.res	1✔
161	assert sub.geo['xmin'] - ras.geo['xmin'] == xres	1✔
162	assert ras.geo['xmax'] - sub.geo['xmax'] == xres	1✔
163	assert sub.geo['ymin'] - ras.geo['ymin'] == xres	1✔
164	assert ras.geo['ymax'] - sub.geo['ymax'] == xres	1✔
165
166
167	def test_Raster_extract(testdata):	1✔
168	with Raster(testdata['tif']) as ras:	1✔
169	assert ras.extract(px=624000, py=4830000, radius=5) == -10.48837461270875	1✔
170	with pytest.raises(RuntimeError):	1✔
171	ras.extract(1, 4830000)	1✔
172	with pytest.raises(RuntimeError):	1✔
173	ras.extract(624000, 1)	1✔
174
175	# ensure corner extraction capability
176	assert ras.extract(px=ras.geo['xmin'], py=ras.geo['ymax']) == -10.147890090942383	1✔
177	assert ras.extract(px=ras.geo['xmin'], py=ras.geo['ymin']) == -14.640368461608887	1✔
178	assert ras.extract(px=ras.geo['xmax'], py=ras.geo['ymax']) == -9.599242210388182	1✔
179	assert ras.extract(px=ras.geo['xmax'], py=ras.geo['ymin']) == -9.406558990478516	1✔
180
181	# test nodata handling capability and correct indexing
182	mat = ras.matrix()	1✔
183	mat[0:10, 0:10] = ras.nodata	1✔
184	mat[207:217, 258:268] = ras.nodata	1✔
185	ras.assign(mat, band=0)	1✔
186	assert ras.extract(px=ras.geo['xmin'], py=ras.geo['ymax'], radius=5) == ras.nodata	1✔
187	assert ras.extract(px=ras.geo['xmax'], py=ras.geo['ymin'], radius=5) == ras.nodata	1✔
188
189
190	def test_Raster_filestack(testdata):	1✔
191	with pytest.raises(RuntimeError):	1✔
192	with Raster([testdata['tif']]) as ras:	1✔
UNCOV 193	print(ras)	×
194	with Raster([testdata['tif'], testdata['tif2']]) as ras:	1✔
195	assert ras.bands == 2	1✔
196	arr = ras.array()	1✔
197	mean = parallel_apply_along_axis(np.nanmean, axis=2, arr=arr, cores=4)	1✔
198	assert mean.shape == (217, 268)	1✔
199
200
201	def test_dtypes():	1✔
202	assert Dtype('Float32').gdalint == 6	1✔
203	assert Dtype(6).gdalstr == 'Float32'	1✔
204	assert Dtype('uint32').gdalstr == 'UInt32'	1✔
205	with pytest.raises(ValueError):	1✔
206	Dtype('foobar')	1✔
207	with pytest.raises(ValueError):	1✔
208	Dtype(999)	1✔
209	with pytest.raises(TypeError):	1✔
210	Dtype(None)	1✔
211
212
213	def test_stack(tmpdir, testdata):	1✔
214	name = testdata['tif']	1✔
215	outname = os.path.join(str(tmpdir), 'test')	1✔
216	tr = (30, 30)	1✔
217	# no input files provided
218	with pytest.raises(RuntimeError):	1✔
219	stack(srcfiles=[], resampling='near', targetres=tr,	1✔
220	srcnodata=-99, dstnodata=-99, dstfile=outname)
221
222	# two files, but only one layer name
223	with pytest.raises(RuntimeError):	1✔
224	stack(srcfiles=[name, name], resampling='near', targetres=tr,	1✔
225	srcnodata=-99, dstnodata=-99, dstfile=outname, layernames=['a'])
226
227	# targetres must be a two-entry tuple/list
228	with pytest.raises(RuntimeError):	1✔
229	stack(srcfiles=[name, name], resampling='near', targetres=30,	1✔
230	srcnodata=-99, dstnodata=-99, dstfile=outname)
231
232	# only one file specified
233	with pytest.raises(RuntimeError):	1✔
234	stack(srcfiles=[name], resampling='near', targetres=tr, overwrite=True,	1✔
235	srcnodata=-99, dstnodata=-99, dstfile=outname)
236
237	# targetres must contain two values
238	with pytest.raises(RuntimeError):	1✔
239	stack(srcfiles=[name, name], resampling='near', targetres=(30, 30, 30),	1✔
240	srcnodata=-99, dstnodata=-99, dstfile=outname)
241
242	# unknown resampling method
243	with pytest.raises(RuntimeError):	1✔
244	stack(srcfiles=[name, name], resampling='foobar', targetres=tr,	1✔
245	srcnodata=-99, dstnodata=-99, dstfile=outname)
246
247	# non-existing files
248	with pytest.raises(RuntimeError):	1✔
249	stack(srcfiles=['foo', 'bar'], resampling='near', targetres=tr,	1✔
250	srcnodata=-99, dstnodata=-99, dstfile=outname)
251
252	# create a multi-band stack
253	stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=True,	1✔
254	srcnodata=-99, dstnodata=-99, dstfile=outname, layernames=['test1', 'test2'])
255	with Raster(outname) as ras:	1✔
256	assert ras.bands == 2	1✔
257	# Raster.rescale currently only supports one band
258	with pytest.raises(ValueError):	1✔
259	ras.rescale(lambda x: x * 10)	1✔
260
261	# outname exists and overwrite is False
262	with pytest.raises(RuntimeError):	1✔
263	stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=False,	1✔
264	srcnodata=-99, dstnodata=-99, dstfile=outname, layernames=['test1', 'test2'])
265
266	# pass shapefile
267	outname = os.path.join(str(tmpdir), 'test2')	1✔
268	with Raster(name).bbox() as box:	1✔
269	stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=True,	1✔
270	srcnodata=-99, dstnodata=-99, dstfile=outname, shapefile=box, layernames=['test1', 'test2'])
271	with Raster(outname) as ras:	1✔
272	assert ras.bands == 2	1✔
273
274	# pass shapefile and do mosaicing
275	outname = os.path.join(str(tmpdir), 'test3')	1✔
276	with Raster(name).bbox() as box:	1✔
277	stack(srcfiles=[[name, name]], resampling='near', targetres=tr, overwrite=True,	1✔
278	srcnodata=-99, dstnodata=-99, dstfile=outname, shapefile=box)
279	with Raster(outname + '.tif') as ras:	1✔
280	assert ras.bands == 1	1✔
281	assert ras.format == 'GTiff'	1✔
282
283	# projection mismatch
284	name2 = os.path.join(str(tmpdir), os.path.basename(name))	1✔
285	outname = os.path.join(str(tmpdir), 'test4')	1✔
286	gdalwarp(src=name, dst=name2, dstSRS=crsConvert(4326, 'wkt'))	1✔
287	with pytest.raises(RuntimeError):	1✔
288	stack(srcfiles=[name, name2], resampling='near', targetres=tr, overwrite=True,	1✔
289	srcnodata=-99, dstnodata=-99, dstfile=outname)
290
291	# no projection found
292	outname = os.path.join(str(tmpdir), 'test5')	1✔
293	gdal_translate(src=name, dst=name2, options=['-co', 'PROFILE=BASELINE'])	1✔
294	with Raster(name2) as ras:	1✔
295	print(ras.projection)	1✔
296	with pytest.raises(RuntimeError):	1✔
297	stack(srcfiles=[name2, name2], resampling='near', targetres=tr, overwrite=True,	1✔
298	srcnodata=-99, dstnodata=-99, dstfile=outname)
299
300	# create separate GeoTiffs
301	outdir = os.path.join(str(tmpdir), 'subdir')	1✔
302	stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=True, layernames=['test1', 'test2'],	1✔
303	srcnodata=-99, dstnodata=-99, dstfile=outdir, separate=True, compress=True)
304
305	# repeat with overwrite disabled (no error raised, just a print message)
306	stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=False, layernames=['test1', 'test2'],	1✔
307	srcnodata=-99, dstnodata=-99, dstfile=outdir, separate=True, compress=True)
308
309	# repeat without layernames but sortfun
310	# bandnames not unique
311	outdir = os.path.join(str(tmpdir), 'subdir2')	1✔
312	with pytest.raises(RuntimeError):	1✔
313	stack(srcfiles=[name, name], resampling='near', targetres=tr, overwrite=True, sortfun=os.path.basename,	1✔
314	srcnodata=-99, dstnodata=-99, dstfile=outdir, separate=True, compress=True)
315
316	# repeat without layernames but sortfun
317	name2 = os.path.join(str(tmpdir), os.path.basename(name).replace('VV', 'XX'))	1✔
318	shutil.copyfile(name, name2)	1✔
319	outdir = os.path.join(str(tmpdir), 'subdir2')	1✔
320	stack(srcfiles=[name, name2], resampling='near', targetres=tr, overwrite=True, sortfun=os.path.basename,	1✔
321	srcnodata=-99, dstnodata=-99, dstfile=outdir, separate=True, compress=True)
322
323	# shapefile filtering
324	outdir = os.path.join(str(tmpdir), 'subdir3')	1✔
325	files = [testdata['tif'], testdata['tif2'], testdata['tif3']]	1✔
326	with Raster(files[0]).bbox() as box:	1✔
327	stack(srcfiles=files, resampling='near', targetres=(30, 30),	1✔
328	overwrite=False, layernames=['test1', 'test2', 'test3'],
329	srcnodata=-99, dstnodata=-99, dstfile=outdir,
330	separate=True, compress=True, shapefile=box)
331	# repeated run with different scene selection and only one scene after spatial filtering
332	stack(srcfiles=files[1:], resampling='near', targetres=(30, 30),	1✔
333	overwrite=True, layernames=['test2', 'test3'],
334	srcnodata=-99, dstnodata=-99, dstfile=outdir,
335	separate=True, compress=True, shapefile=box)
336
337
338	def test_auxil(tmpdir, testdata):	1✔
339	dir = str(tmpdir)	1✔
340	with Raster(testdata['tif']) as ras:	1✔
341	bbox = os.path.join(dir, 'bbox.shp')	1✔
342	ras.bbox(bbox)	1✔
343	ogr2ogr(src=bbox, dst=os.path.join(dir, 'bbox.gml'), format='GML')	1✔
344	gdal_translate(src=ras.raster, dst=os.path.join(dir, 'test'), format='ENVI')	1✔
345	gdal_rasterize(src=bbox, dst=os.path.join(dir, 'test2'), format='GTiff', xRes=20, yRes=20)	1✔
346
347
348	def test_auxil_coordinate_reproject():	1✔
349	point = coordinate_reproject(x=11, y=51, s_crs=4326, t_crs=32632)	1✔
350	assert round(point[0], 3) == 640333.296	1✔
351	assert round(point[1], 3) == 5651728.683	1✔
352
353
354	def test_auxil_cmap_mpl2gdal():	1✔
355	cmap = cmap_mpl2gdal(mplcolor='YlGnBu', values=range(0, 100))	1✔
356	assert type(cmap) == gdal.ColorTable	1✔
357
358
359	def test_rasterize(tmpdir, testdata):	1✔
360	outname = os.path.join(str(tmpdir), 'test.shp')	1✔
361	with Raster(testdata['tif']) as ras:	1✔
362	vec = ras.bbox()	1✔
363
364	# test length mismatch between burn_values and expressions
365	with pytest.raises(RuntimeError):	1✔
366	rasterize(vec, reference=ras, outname=outname, burn_values=[1], expressions=['foo', 'bar'])	1✔
367
368	# test a faulty expression
369	with pytest.raises(RuntimeError):	1✔
370	rasterize(vec, reference=ras, outname=outname, burn_values=[1], expressions=['foo'])	1✔
371
372	# test default parametrization
373	rasterize(vec, reference=ras, outname=outname)	1✔
374	assert os.path.isfile(outname)	1✔
375
376	# test appending to existing file with valid expression
377	rasterize(vec, reference=ras, outname=outname, append=True, burn_values=[1], expressions=['area=23262400.0'])	1✔
378
379	# test wrong input type for reference
380	with pytest.raises(RuntimeError):	1✔
381	rasterize(vec, reference='foobar', outname=outname)	1✔
382
383
384	def test_envi(tmpdir):	1✔
385	with pytest.raises(RuntimeError):	1✔
386	obj = HDRobject(1)	1✔
387	with pytest.raises(RuntimeError):	1✔
388	obj = HDRobject('foobar')	1✔
389	outname = os.path.join(str(tmpdir), 'test')	1✔
390	with HDRobject() as header:	1✔
391	header.band_names = ['band1']	1✔
392	header.write(outname)	1✔
393	outname += '.hdr'	1✔
394	with HDRobject(outname) as header:	1✔
395	assert header.band_names == ['band1']	1✔
396	vals = vars(header)	1✔
397	with HDRobject(vals) as header:	1✔
398	assert header.byte_order == 0	1✔
399	hdr(vals, outname + '2')	1✔
400
401
402	def test_sqlite():	1✔
403	with pytest.raises(RuntimeError):	1✔
404	con = sqlite_setup(extensions='spatialite')	1✔
405	con = sqlite_setup(extensions=['spatialite'])	1✔
406	con.close()	1✔
407	con = __Handler()	1✔
408	assert sorted(con.version.keys()) == ['sqlite']	1✔
409
410	con = __Handler(extensions=['spatialite'])	1✔
411	assert sorted(con.version.keys()) == ['spatialite', 'sqlite']	1✔
412	assert 'spatial_ref_sys' in con.get_tablenames()	1✔
413
414
415	def test_png(tmpdir, testdata):	1✔
416	outname = os.path.join(str(tmpdir), 'test')	1✔
417	with Raster(testdata['tif']) as ras:	1✔
418	png(src=ras, dst=outname, percent=100, scale=(2, 98), worldfile=True)	1✔
419	assert os.path.isfile(outname + '.png')	1✔
420
421	with pytest.raises(TypeError):	1✔
422	png(src=testdata['tif'], dst=outname, percent=100, scale=(2, 98), worldfile=True)	1✔
423
424	src = [testdata['tif'], testdata['tif2']]	1✔
425	with pytest.raises(ValueError):	1✔
426	with Raster(src) as ras:	1✔
427	png(src=ras, dst=outname, percent=100, scale=(2, 98), worldfile=True)	1✔
428
429	src.append(testdata['tif3'])	1✔
430	outname = os.path.join(str(tmpdir), 'test_rgb.png')	1✔
431	with Raster(src) as ras:	1✔
432	png(src=ras, dst=outname, percent=100, scale=(2, 98), worldfile=True)	1✔
433
434
435	def test_addfield():	1✔
436	extent = {'xmin': 10, 'xmax': 11, 'ymin': 50, 'ymax': 51}	1✔
437	with bbox(coordinates=extent, crs=4326) as box:	1✔
438	box.addfield(name='test1', type=ogr.OFTString, values=['a'])	1✔
439	box.addfield(name='test2', type=ogr.OFTStringList, values=[['a', 'b']])	1✔
440	box.addfield(name='test3', type=ogr.OFTInteger, values=[1])	1✔
441	box.addfield(name='test4', type=ogr.OFTIntegerList, values=[[1, 2]])	1✔
442	box.addfield(name='test5', type=ogr.OFTInteger64, values=[1])	1✔
443	box.addfield(name='test6', type=ogr.OFTInteger64List, values=[[1, 2]])	1✔
444	box.addfield(name='test7', type=ogr.OFTReal, values=[1])	1✔
445	box.addfield(name='test8', type=ogr.OFTRealList, values=[[1., 2.]])	1✔
446	box.addfield(name='test9', type=ogr.OFTBinary, values=[b'1'])	1✔
447
448
449	def test_wkt2vector():	1✔
450	wkt1 = 'POLYGON ((0. 0., 0. 1., 1. 1., 1. 0., 0. 0.))'	1✔
451	wkt2 = 'POLYGON ((1. 1., 1. 2., 2. 2., 2. 1., 1. 1.))'	1✔
452	with wkt2vector(wkt1, srs=4326) as vec:	1✔
453	assert vec.getArea() == 1.	1✔
454	with wkt2vector([wkt1, wkt2], srs=4326) as vec:	1✔
455	assert vec.getArea() == 2.	1✔

johntruckenbrodt / spatialist / 14332609845

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous