1. import os

  2. import shutil

  3. import struct

  4. import tempfile

  5. import zipfile

  6. from unittest import mock

  7. 

  8. from django.contrib.gis.gdal import GDALRaster, SpatialReference

  9. from django.contrib.gis.gdal.error import GDALException

  10. from django.contrib.gis.gdal.raster.band import GDALBand

  11. from django.contrib.gis.shortcuts import numpy

  12. from django.test import SimpleTestCase

  13. 

  14. from ..data.rasters.textrasters import JSON_RASTER

  15. 

  16. 

  17. class GDALRasterTests(SimpleTestCase):

  18.     """

  19.     Test a GDALRaster instance created from a file (GeoTiff).

  20.     """

  21. 

  22.     def setUp(self):

  23.         self.rs_path = os.path.join(

  24.             os.path.dirname(__file__), "../data/rasters/raster.tif"

  25.         )

  26.         self.rs = GDALRaster(self.rs_path)

  27. 

  28.     def test_rs_name_repr(self):

  29.         self.assertEqual(self.rs_path, self.rs.name)

  30.         self.assertRegex(repr(self.rs), r"<Raster object at 0x\w+>")

  31. 

  32.     def test_rs_driver(self):

  33.         self.assertEqual(self.rs.driver.name, "GTiff")

  34. 

  35.     def test_rs_size(self):

  36.         self.assertEqual(self.rs.width, 163)

  37.         self.assertEqual(self.rs.height, 174)

  38. 

  39.     def test_rs_srs(self):

  40.         self.assertEqual(self.rs.srs.srid, 3086)

  41.         self.assertEqual(self.rs.srs.units, (1.0, "metre"))

  42. 

  43.     def test_rs_srid(self):

  44.         rast = GDALRaster(

  45.             {

  46.                 "width": 16,

  47.                 "height": 16,

  48.                 "srid": 4326,

  49.             }

  50.         )

  51.         self.assertEqual(rast.srid, 4326)

  52.         rast.srid = 3086

  53.         self.assertEqual(rast.srid, 3086)

  54. 

  55.     def test_geotransform_and_friends(self):

  56.         # Assert correct values for file based raster

  57.         self.assertEqual(

  58.             self.rs.geotransform,

  59.             [511700.4680706557, 100.0, 0.0, 435103.3771231986, 0.0, -100.0],

  60.         )

  61.         self.assertEqual(self.rs.origin, [511700.4680706557, 435103.3771231986])

  62.         self.assertEqual(self.rs.origin.x, 511700.4680706557)

  63.         self.assertEqual(self.rs.origin.y, 435103.3771231986)

  64.         self.assertEqual(self.rs.scale, [100.0, -100.0])

  65.         self.assertEqual(self.rs.scale.x, 100.0)

  66.         self.assertEqual(self.rs.scale.y, -100.0)

  67.         self.assertEqual(self.rs.skew, [0, 0])

  68.         self.assertEqual(self.rs.skew.x, 0)

  69.         self.assertEqual(self.rs.skew.y, 0)

  70.         # Create in-memory rasters and change gtvalues

  71.         rsmem = GDALRaster(JSON_RASTER)

  72.         # geotransform accepts both floats and ints

  73.         rsmem.geotransform = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]

  74.         self.assertEqual(rsmem.geotransform, [0.0, 1.0, 2.0, 3.0, 4.0, 5.0])

  75.         rsmem.geotransform = range(6)

  76.         self.assertEqual(rsmem.geotransform, [float(x) for x in range(6)])

  77.         self.assertEqual(rsmem.origin, [0, 3])

  78.         self.assertEqual(rsmem.origin.x, 0)

  79.         self.assertEqual(rsmem.origin.y, 3)

  80.         self.assertEqual(rsmem.scale, [1, 5])

  81.         self.assertEqual(rsmem.scale.x, 1)

  82.         self.assertEqual(rsmem.scale.y, 5)

  83.         self.assertEqual(rsmem.skew, [2, 4])

  84.         self.assertEqual(rsmem.skew.x, 2)

  85.         self.assertEqual(rsmem.skew.y, 4)

  86.         self.assertEqual(rsmem.width, 5)

  87.         self.assertEqual(rsmem.height, 5)

  88. 

  89.     def test_geotransform_bad_inputs(self):

  90.         rsmem = GDALRaster(JSON_RASTER)

  91.         error_geotransforms = [

  92.             [1, 2],

  93.             [1, 2, 3, 4, 5, "foo"],

  94.             [1, 2, 3, 4, 5, 6, "foo"],

  95.         ]

  96.         msg = "Geotransform must consist of 6 numeric values."

  97.         for geotransform in error_geotransforms:

  98.             with self.subTest(i=geotransform), self.assertRaisesMessage(

  99.                 ValueError, msg

  100.             ):

  101.                 rsmem.geotransform = geotransform

  102. 

  103.     def test_rs_extent(self):

  104.         self.assertEqual(

  105.             self.rs.extent,

  106.             (

  107.                 511700.4680706557,

  108.                 417703.3771231986,

  109.                 528000.4680706557,

  110.                 435103.3771231986,

  111.             ),

  112.         )

  113. 

  114.     def test_rs_bands(self):

  115.         self.assertEqual(len(self.rs.bands), 1)

  116.         self.assertIsInstance(self.rs.bands[0], GDALBand)

  117. 

  118.     def test_memory_based_raster_creation(self):

  119.         # Create uint8 raster with full pixel data range (0-255)

  120.         rast = GDALRaster(

  121.             {

  122.                 "datatype": 1,

  123.                 "width": 16,

  124.                 "height": 16,

  125.                 "srid": 4326,

  126.                 "bands": [

  127.                     {

  128.                         "data": range(256),

  129.                         "nodata_value": 255,

  130.                     }

  131.                 ],

  132.             }

  133.         )

  134. 

  135.         # Get array from raster

  136.         result = rast.bands[0].data()

  137.         if numpy:

  138.             result = result.flatten().tolist()

  139. 

  140.         # Assert data is same as original input

  141.         self.assertEqual(result, list(range(256)))

  142. 

  143.     def test_file_based_raster_creation(self):

  144.         # Prepare tempfile

  145.         rstfile = tempfile.NamedTemporaryFile(suffix=".tif")

  146. 

  147.         # Create file-based raster from scratch

  148.         GDALRaster(

  149.             {

  150.                 "datatype": self.rs.bands[0].datatype(),

  151.                 "driver": "tif",

  152.                 "name": rstfile.name,

  153.                 "width": 163,

  154.                 "height": 174,

  155.                 "nr_of_bands": 1,

  156.                 "srid": self.rs.srs.wkt,

  157.                 "origin": (self.rs.origin.x, self.rs.origin.y),

  158.                 "scale": (self.rs.scale.x, self.rs.scale.y),

  159.                 "skew": (self.rs.skew.x, self.rs.skew.y),

  160.                 "bands": [

  161.                     {

  162.                         "data": self.rs.bands[0].data(),

  163.                         "nodata_value": self.rs.bands[0].nodata_value,

  164.                     }

  165.                 ],

  166.             }

  167.         )

  168. 

  169.         # Reload newly created raster from file

  170.         restored_raster = GDALRaster(rstfile.name)

  171.         # Presence of TOWGS84 depend on GDAL/Proj versions.

  172.         self.assertEqual(

  173.             restored_raster.srs.wkt.replace("TOWGS84[0,0,0,0,0,0,0],", ""),

  174.             self.rs.srs.wkt.replace("TOWGS84[0,0,0,0,0,0,0],", ""),

  175.         )

  176.         self.assertEqual(restored_raster.geotransform, self.rs.geotransform)

  177.         if numpy:

  178.             numpy.testing.assert_equal(

  179.                 restored_raster.bands[0].data(), self.rs.bands[0].data()

  180.             )

  181.         else:

  182.             self.assertEqual(restored_raster.bands[0].data(), self.rs.bands[0].data())

  183. 

  184.     def test_nonexistent_file(self):

  185.         msg = 'Unable to read raster source input "nonexistent.tif".'

  186.         with self.assertRaisesMessage(GDALException, msg):

  187.             GDALRaster("nonexistent.tif")

  188. 

  189.     def test_vsi_raster_creation(self):

  190.         # Open a raster as a file object.

  191.         with open(self.rs_path, "rb") as dat:

  192.             # Instantiate a raster from the file binary buffer.

  193.             vsimem = GDALRaster(dat.read())

  194.         # The data of the in-memory file is equal to the source file.

  195.         result = vsimem.bands[0].data()

  196.         target = self.rs.bands[0].data()

  197.         if numpy:

  198.             result = result.flatten().tolist()

  199.             target = target.flatten().tolist()

  200.         self.assertEqual(result, target)

  201. 

  202.     def test_vsi_raster_deletion(self):

  203.         path = "/vsimem/raster.tif"

  204.         # Create a vsi-based raster from scratch.

  205.         vsimem = GDALRaster(

  206.             {

  207.                 "name": path,

  208.                 "driver": "tif",

  209.                 "width": 4,

  210.                 "height": 4,

  211.                 "srid": 4326,

  212.                 "bands": [

  213.                     {

  214.                         "data": range(16),

  215.                     }

  216.                 ],

  217.             }

  218.         )

  219.         # The virtual file exists.

  220.         rst = GDALRaster(path)

  221.         self.assertEqual(rst.width, 4)

  222.         # Delete GDALRaster.

  223.         del vsimem

  224.         del rst

  225.         # The virtual file has been removed.

  226.         msg = 'Could not open the datasource at "/vsimem/raster.tif"'

  227.         with self.assertRaisesMessage(GDALException, msg):

  228.             GDALRaster(path)

  229. 

  230.     def test_vsi_invalid_buffer_error(self):

  231.         msg = "Failed creating VSI raster from the input buffer."

  232.         with self.assertRaisesMessage(GDALException, msg):

  233.             GDALRaster(b"not-a-raster-buffer")

  234. 

  235.     def test_vsi_buffer_property(self):

  236.         # Create a vsi-based raster from scratch.

  237.         rast = GDALRaster(

  238.             {

  239.                 "name": "/vsimem/raster.tif",

  240.                 "driver": "tif",

  241.                 "width": 4,

  242.                 "height": 4,

  243.                 "srid": 4326,

  244.                 "bands": [

  245.                     {

  246.                         "data": range(16),

  247.                     }

  248.                 ],

  249.             }

  250.         )

  251.         # Do a round trip from raster to buffer to raster.

  252.         result = GDALRaster(rast.vsi_buffer).bands[0].data()

  253.         if numpy:

  254.             result = result.flatten().tolist()

  255.         # Band data is equal to nodata value except on input block of ones.

  256.         self.assertEqual(result, list(range(16)))

  257.         # The vsi buffer is None for rasters that are not vsi based.

  258.         self.assertIsNone(self.rs.vsi_buffer)

  259. 

  260.     def test_vsi_vsizip_filesystem(self):

  261.         rst_zipfile = tempfile.NamedTemporaryFile(suffix=".zip")

  262.         with zipfile.ZipFile(rst_zipfile, mode="w") as zf:

  263.             zf.write(self.rs_path, "raster.tif")

  264.         rst_path = "/vsizip/" + os.path.join(rst_zipfile.name, "raster.tif")

  265.         rst = GDALRaster(rst_path)

  266.         self.assertEqual(rst.driver.name, self.rs.driver.name)

  267.         self.assertEqual(rst.name, rst_path)

  268.         self.assertIs(rst.is_vsi_based, True)

  269.         self.assertIsNone(rst.vsi_buffer)

  270. 

  271.     def test_offset_size_and_shape_on_raster_creation(self):

  272.         rast = GDALRaster(

  273.             {

  274.                 "datatype": 1,

  275.                 "width": 4,

  276.                 "height": 4,

  277.                 "srid": 4326,

  278.                 "bands": [

  279.                     {

  280.                         "data": (1,),

  281.                         "offset": (1, 1),

  282.                         "size": (2, 2),

  283.                         "shape": (1, 1),

  284.                         "nodata_value": 2,

  285.                     }

  286.                 ],

  287.             }

  288.         )

  289.         # Get array from raster.

  290.         result = rast.bands[0].data()

  291.         if numpy:

  292.             result = result.flatten().tolist()

  293.         # Band data is equal to nodata value except on input block of ones.

  294.         self.assertEqual(result, [2, 2, 2, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 2, 2, 2])

  295. 

  296.     def test_set_nodata_value_on_raster_creation(self):

  297.         # Create raster filled with nodata values.

  298.         rast = GDALRaster(

  299.             {

  300.                 "datatype": 1,

  301.                 "width": 2,

  302.                 "height": 2,

  303.                 "srid": 4326,

  304.                 "bands": [{"nodata_value": 23}],

  305.             }

  306.         )

  307.         # Get array from raster.

  308.         result = rast.bands[0].data()

  309.         if numpy:

  310.             result = result.flatten().tolist()

  311.         # All band data is equal to nodata value.

  312.         self.assertEqual(result, [23] * 4)

  313. 

  314.     def test_set_nodata_none_on_raster_creation(self):

  315.         # Create raster without data and without nodata value.

  316.         rast = GDALRaster(

  317.             {

  318.                 "datatype": 1,

  319.                 "width": 2,

  320.                 "height": 2,

  321.                 "srid": 4326,

  322.                 "bands": [{"nodata_value": None}],

  323.             }

  324.         )

  325.         # Get array from raster.

  326.         result = rast.bands[0].data()

  327.         if numpy:

  328.             result = result.flatten().tolist()

  329.         # Band data is equal to zero because no nodata value has been specified.

  330.         self.assertEqual(result, [0] * 4)

  331. 

  332.     def test_raster_metadata_property(self):

  333.         data = self.rs.metadata

  334.         self.assertEqual(data["DEFAULT"], {"AREA_OR_POINT": "Area"})

  335.         self.assertEqual(data["IMAGE_STRUCTURE"], {"INTERLEAVE": "BAND"})

  336. 

  337.         # Create file-based raster from scratch

  338.         source = GDALRaster(

  339.             {

  340.                 "datatype": 1,

  341.                 "width": 2,

  342.                 "height": 2,

  343.                 "srid": 4326,

  344.                 "bands": [{"data": range(4), "nodata_value": 99}],

  345.             }

  346.         )

  347.         # Set metadata on raster and on a band.

  348.         metadata = {

  349.             "DEFAULT": {"OWNER": "Django", "VERSION": "1.0", "AREA_OR_POINT": "Point"},

  350.         }

  351.         source.metadata = metadata

  352.         source.bands[0].metadata = metadata

  353.         self.assertEqual(source.metadata["DEFAULT"], metadata["DEFAULT"])

  354.         self.assertEqual(source.bands[0].metadata["DEFAULT"], metadata["DEFAULT"])

  355.         # Update metadata on raster.

  356.         metadata = {

  357.             "DEFAULT": {"VERSION": "2.0"},

  358.         }

  359.         source.metadata = metadata

  360.         self.assertEqual(source.metadata["DEFAULT"]["VERSION"], "2.0")

  361.         # Remove metadata on raster.

  362.         metadata = {

  363.             "DEFAULT": {"OWNER": None},

  364.         }

  365.         source.metadata = metadata

  366.         self.assertNotIn("OWNER", source.metadata["DEFAULT"])

  367. 

  368.     def test_raster_info_accessor(self):

  369.         infos = self.rs.info

  370.         # Data

  371.         info_lines = [line.strip() for line in infos.split("\n") if line.strip() != ""]

  372.         for line in [

  373.             "Driver: GTiff/GeoTIFF",

  374.             "Files: {}".format(self.rs_path),

  375.             "Size is 163, 174",

  376.             "Origin = (511700.468070655711927,435103.377123198588379)",

  377.             "Pixel Size = (100.000000000000000,-100.000000000000000)",

  378.             "Metadata:",

  379.             "AREA_OR_POINT=Area",

  380.             "Image Structure Metadata:",

  381.             "INTERLEAVE=BAND",

  382.             "Band 1 Block=163x50 Type=Byte, ColorInterp=Gray",

  383.             "NoData Value=15",

  384.         ]:

  385.             self.assertIn(line, info_lines)

  386.         for line in [

  387.             r"Upper Left  \(  511700.468,  435103.377\) "

  388.             r'\( 82d51\'46.1\d"W, 27d55\' 1.5\d"N\)',

  389.             r"Lower Left  \(  511700.468,  417703.377\) "

  390.             r'\( 82d51\'52.0\d"W, 27d45\'37.5\d"N\)',

  391.             r"Upper Right \(  528000.468,  435103.377\) "

  392.             r'\( 82d41\'48.8\d"W, 27d54\'56.3\d"N\)',

  393.             r"Lower Right \(  528000.468,  417703.377\) "

  394.             r'\( 82d41\'55.5\d"W, 27d45\'32.2\d"N\)',

  395.             r"Center      \(  519850.468,  426403.377\) "

  396.             r'\( 82d46\'50.6\d"W, 27d50\'16.9\d"N\)',

  397.         ]:

  398.             self.assertRegex(infos, line)

  399.         # CRS (skip the name because string depends on the GDAL/Proj versions).

  400.         self.assertIn("NAD83 / Florida GDL Albers", infos)

  401. 

  402.     def test_compressed_file_based_raster_creation(self):

  403.         rstfile = tempfile.NamedTemporaryFile(suffix=".tif")

  404.         # Make a compressed copy of an existing raster.

  405.         compressed = self.rs.warp(

  406.             {"papsz_options": {"compress": "packbits"}, "name": rstfile.name}

  407.         )

  408.         # Check physically if compression worked.

  409.         self.assertLess(os.path.getsize(compressed.name), os.path.getsize(self.rs.name))

  410.         # Create file-based raster with options from scratch.

  411.         compressed = GDALRaster(

  412.             {

  413.                 "datatype": 1,

  414.                 "driver": "tif",

  415.                 "name": rstfile.name,

  416.                 "width": 40,

  417.                 "height": 40,

  418.                 "srid": 3086,

  419.                 "origin": (500000, 400000),

  420.                 "scale": (100, -100),

  421.                 "skew": (0, 0),

  422.                 "bands": [

  423.                     {

  424.                         "data": range(40 ^ 2),

  425.                         "nodata_value": 255,

  426.                     }

  427.                 ],

  428.                 "papsz_options": {

  429.                     "compress": "packbits",

  430.                     "pixeltype": "signedbyte",

  431.                     "blockxsize": 23,

  432.                     "blockysize": 23,

  433.                 },

  434.             }

  435.         )

  436.         # Check if options used on creation are stored in metadata.

  437.         # Reopening the raster ensures that all metadata has been written

  438.         # to the file.

  439.         compressed = GDALRaster(compressed.name)

  440.         self.assertEqual(

  441.             compressed.metadata["IMAGE_STRUCTURE"]["COMPRESSION"],

  442.             "PACKBITS",

  443.         )

  444.         self.assertEqual(

  445.             compressed.bands[0].metadata["IMAGE_STRUCTURE"]["PIXELTYPE"], "SIGNEDBYTE"

  446.         )

  447.         self.assertIn("Block=40x23", compressed.info)

  448. 

  449.     def test_raster_warp(self):

  450.         # Create in memory raster

  451.         source = GDALRaster(

  452.             {

  453.                 "datatype": 1,

  454.                 "driver": "MEM",

  455.                 "name": "sourceraster",

  456.                 "width": 4,

  457.                 "height": 4,

  458.                 "nr_of_bands": 1,

  459.                 "srid": 3086,

  460.                 "origin": (500000, 400000),

  461.                 "scale": (100, -100),

  462.                 "skew": (0, 0),

  463.                 "bands": [

  464.                     {

  465.                         "data": range(16),

  466.                         "nodata_value": 255,

  467.                     }

  468.                 ],

  469.             }

  470.         )

  471. 

  472.         # Test altering the scale, width, and height of a raster

  473.         data = {

  474.             "scale": [200, -200],

  475.             "width": 2,

  476.             "height": 2,

  477.         }

  478.         target = source.warp(data)

  479.         self.assertEqual(target.width, data["width"])

  480.         self.assertEqual(target.height, data["height"])

  481.         self.assertEqual(target.scale, data["scale"])

  482.         self.assertEqual(target.bands[0].datatype(), source.bands[0].datatype())

  483.         self.assertEqual(target.name, "sourceraster_copy.MEM")

  484.         result = target.bands[0].data()

  485.         if numpy:

  486.             result = result.flatten().tolist()

  487.         self.assertEqual(result, [5, 7, 13, 15])

  488. 

  489.         # Test altering the name and datatype (to float)

  490.         data = {

  491.             "name": "/path/to/targetraster.tif",

  492.             "datatype": 6,

  493.         }

  494.         target = source.warp(data)

  495.         self.assertEqual(target.bands[0].datatype(), 6)

  496.         self.assertEqual(target.name, "/path/to/targetraster.tif")

  497.         self.assertEqual(target.driver.name, "MEM")

  498.         result = target.bands[0].data()

  499.         if numpy:

  500.             result = result.flatten().tolist()

  501.         self.assertEqual(

  502.             result,

  503.             [

  504.                 0.0,

  505.                 1.0,

  506.                 2.0,

  507.                 3.0,

  508.                 4.0,

  509.                 5.0,

  510.                 6.0,

  511.                 7.0,

  512.                 8.0,

  513.                 9.0,

  514.                 10.0,

  515.                 11.0,

  516.                 12.0,

  517.                 13.0,

  518.                 14.0,

  519.                 15.0,

  520.             ],

  521.         )

  522. 

  523.     def test_raster_warp_nodata_zone(self):

  524.         # Create in memory raster.

  525.         source = GDALRaster(

  526.             {

  527.                 "datatype": 1,

  528.                 "driver": "MEM",

  529.                 "width": 4,

  530.                 "height": 4,

  531.                 "srid": 3086,

  532.                 "origin": (500000, 400000),

  533.                 "scale": (100, -100),

  534.                 "skew": (0, 0),

  535.                 "bands": [

  536.                     {

  537.                         "data": range(16),

  538.                         "nodata_value": 23,

  539.                     }

  540.                 ],

  541.             }

  542.         )

  543.         # Warp raster onto a location that does not cover any pixels of the original.

  544.         result = source.warp({"origin": (200000, 200000)}).bands[0].data()

  545.         if numpy:

  546.             result = result.flatten().tolist()

  547.         # The result is an empty raster filled with the correct nodata value.

  548.         self.assertEqual(result, [23] * 16)

  549. 

  550.     def test_raster_clone(self):

  551.         rstfile = tempfile.NamedTemporaryFile(suffix=".tif")

  552.         tests = [

  553.             ("MEM", "", 23),  # In memory raster.

  554.             ("tif", rstfile.name, 99),  # In file based raster.

  555.         ]

  556.         for driver, name, nodata_value in tests:

  557.             with self.subTest(driver=driver):

  558.                 source = GDALRaster(

  559.                     {

  560.                         "datatype": 1,

  561.                         "driver": driver,

  562.                         "name": name,

  563.                         "width": 4,

  564.                         "height": 4,

  565.                         "srid": 3086,

  566.                         "origin": (500000, 400000),

  567.                         "scale": (100, -100),

  568.                         "skew": (0, 0),

  569.                         "bands": [

  570.                             {

  571.                                 "data": range(16),

  572.                                 "nodata_value": nodata_value,

  573.                             }

  574.                         ],

  575.                     }

  576.                 )

  577.                 clone = source.clone()

  578.                 self.assertNotEqual(clone.name, source.name)

  579.                 self.assertEqual(clone._write, source._write)

  580.                 self.assertEqual(clone.srs.srid, source.srs.srid)

  581.                 self.assertEqual(clone.width, source.width)

  582.                 self.assertEqual(clone.height, source.height)

  583.                 self.assertEqual(clone.origin, source.origin)

  584.                 self.assertEqual(clone.scale, source.scale)

  585.                 self.assertEqual(clone.skew, source.skew)

  586.                 self.assertIsNot(clone, source)

  587. 

  588.     def test_raster_transform(self):

  589.         tests = [

  590.             3086,

  591.             "3086",

  592.             SpatialReference(3086),

  593.         ]

  594.         for srs in tests:

  595.             with self.subTest(srs=srs):

  596.                 # Prepare tempfile and nodata value.

  597.                 rstfile = tempfile.NamedTemporaryFile(suffix=".tif")

  598.                 ndv = 99

  599.                 # Create in file based raster.

  600.                 source = GDALRaster(

  601.                     {

  602.                         "datatype": 1,

  603.                         "driver": "tif",

  604.                         "name": rstfile.name,

  605.                         "width": 5,

  606.                         "height": 5,

  607.                         "nr_of_bands": 1,

  608.                         "srid": 4326,

  609.                         "origin": (-5, 5),

  610.                         "scale": (2, -2),

  611.                         "skew": (0, 0),

  612.                         "bands": [

  613.                             {

  614.                                 "data": range(25),

  615.                                 "nodata_value": ndv,

  616.                             }

  617.                         ],

  618.                     }

  619.                 )

  620. 

  621.                 target = source.transform(srs)

  622. 

  623.                 # Reload data from disk.

  624.                 target = GDALRaster(target.name)

  625.                 self.assertEqual(target.srs.srid, 3086)

  626.                 self.assertEqual(target.width, 7)

  627.                 self.assertEqual(target.height, 7)

  628.                 self.assertEqual(target.bands[0].datatype(), source.bands[0].datatype())

  629.                 self.assertAlmostEqual(target.origin[0], 9124842.791079799, 3)

  630.                 self.assertAlmostEqual(target.origin[1], 1589911.6476407414, 3)

  631.                 self.assertAlmostEqual(target.scale[0], 223824.82664250192, 3)

  632.                 self.assertAlmostEqual(target.scale[1], -223824.82664250192, 3)

  633.                 self.assertEqual(target.skew, [0, 0])

  634. 

  635.                 result = target.bands[0].data()

  636.                 if numpy:

  637.                     result = result.flatten().tolist()

  638.                 # The reprojection of a raster that spans over a large area

  639.                 # skews the data matrix and might introduce nodata values.

  640.                 self.assertEqual(

  641.                     result,

  642.                     [

  643.                         ndv,

  644.                         ndv,

  645.                         ndv,

  646.                         ndv,

  647.                         4,

  648.                         ndv,

  649.                         ndv,

  650.                         ndv,

  651.                         ndv,

  652.                         2,

  653.                         3,

  654.                         9,

  655.                         ndv,

  656.                         ndv,

  657.                         ndv,

  658.                         1,

  659.                         2,

  660.                         8,

  661.                         13,

  662.                         19,

  663.                         ndv,

  664.                         0,

  665.                         6,

  666.                         6,

  667.                         12,

  668.                         18,

  669.                         18,

  670.                         24,

  671.                         ndv,

  672.                         10,

  673.                         11,

  674.                         16,

  675.                         22,

  676.                         23,

  677.                         ndv,

  678.                         ndv,

  679.                         ndv,

  680.                         15,

  681.                         21,

  682.                         22,

  683.                         ndv,

  684.                         ndv,

  685.                         ndv,

  686.                         ndv,

  687.                         20,

  688.                         ndv,

  689.                         ndv,

  690.                         ndv,

  691.                         ndv,

  692.                     ],

  693.                 )

  694. 

  695.     def test_raster_transform_clone(self):

  696.         with mock.patch.object(GDALRaster, "clone") as mocked_clone:

  697.             # Create in file based raster.

  698.             rstfile = tempfile.NamedTemporaryFile(suffix=".tif")

  699.             source = GDALRaster(

  700.                 {

  701.                     "datatype": 1,

  702.                     "driver": "tif",

  703.                     "name": rstfile.name,

  704.                     "width": 5,

  705.                     "height": 5,

  706.                     "nr_of_bands": 1,

  707.                     "srid": 4326,

  708.                     "origin": (-5, 5),

  709.                     "scale": (2, -2),

  710.                     "skew": (0, 0),

  711.                     "bands": [

  712.                         {

  713.                             "data": range(25),

  714.                             "nodata_value": 99,

  715.                         }

  716.                     ],

  717.                 }

  718.             )

  719.             # transform() returns a clone because it is the same SRID and

  720.             # driver.

  721.             source.transform(4326)

  722.             self.assertEqual(mocked_clone.call_count, 1)

  723. 

  724.     def test_raster_transform_clone_name(self):

  725.         # Create in file based raster.

  726.         rstfile = tempfile.NamedTemporaryFile(suffix=".tif")

  727.         source = GDALRaster(

  728.             {

  729.                 "datatype": 1,

  730.                 "driver": "tif",

  731.                 "name": rstfile.name,

  732.                 "width": 5,

  733.                 "height": 5,

  734.                 "nr_of_bands": 1,

  735.                 "srid": 4326,

  736.                 "origin": (-5, 5),

  737.                 "scale": (2, -2),

  738.                 "skew": (0, 0),

  739.                 "bands": [

  740.                     {

  741.                         "data": range(25),

  742.                         "nodata_value": 99,

  743.                     }

  744.                 ],

  745.             }

  746.         )

  747.         clone_name = rstfile.name + "_respect_name.GTiff"

  748.         target = source.transform(4326, name=clone_name)

  749.         self.assertEqual(target.name, clone_name)

  750. 

  751. 

  752. class GDALBandTests(SimpleTestCase):

  753.     rs_path = os.path.join(os.path.dirname(__file__), "../data/rasters/raster.tif")

  754. 

  755.     def test_band_data(self):

  756.         rs = GDALRaster(self.rs_path)

  757.         band = rs.bands[0]

  758.         self.assertEqual(band.width, 163)

  759.         self.assertEqual(band.height, 174)

  760.         self.assertEqual(band.description, "")

  761.         self.assertEqual(band.datatype(), 1)

  762.         self.assertEqual(band.datatype(as_string=True), "GDT_Byte")

  763.         self.assertEqual(band.color_interp(), 1)

  764.         self.assertEqual(band.color_interp(as_string=True), "GCI_GrayIndex")

  765.         self.assertEqual(band.nodata_value, 15)

  766.         if numpy:

  767.             data = band.data()

  768.             assert_array = numpy.loadtxt(

  769.                 os.path.join(

  770.                     os.path.dirname(__file__), "../data/rasters/raster.numpy.txt"

  771.                 )

  772.             )

  773.             numpy.testing.assert_equal(data, assert_array)

  774.             self.assertEqual(data.shape, (band.height, band.width))

  775. 

  776.     def test_band_statistics(self):

  777.         with tempfile.TemporaryDirectory() as tmp_dir:

  778.             rs_path = os.path.join(tmp_dir, "raster.tif")

  779.             shutil.copyfile(self.rs_path, rs_path)

  780.             rs = GDALRaster(rs_path)

  781.             band = rs.bands[0]

  782.             pam_file = rs_path + ".aux.xml"

  783.             smin, smax, smean, sstd = band.statistics(approximate=True)

  784.             self.assertEqual(smin, 0)

  785.             self.assertEqual(smax, 9)

  786.             self.assertAlmostEqual(smean, 2.842331288343558)

  787.             self.assertAlmostEqual(sstd, 2.3965567248965356)

  788. 

  789.             smin, smax, smean, sstd = band.statistics(approximate=False, refresh=True)

  790.             self.assertEqual(smin, 0)

  791.             self.assertEqual(smax, 9)

  792.             self.assertAlmostEqual(smean, 2.828326634228898)

  793.             self.assertAlmostEqual(sstd, 2.4260526986669095)

  794. 

  795.             self.assertEqual(band.min, 0)

  796.             self.assertEqual(band.max, 9)

  797.             self.assertAlmostEqual(band.mean, 2.828326634228898)

  798.             self.assertAlmostEqual(band.std, 2.4260526986669095)

  799. 

  800.             # Statistics are persisted into PAM file on band close

  801.             rs = band = None

  802.             self.assertTrue(os.path.isfile(pam_file))

  803. 

  804.     def _remove_aux_file(self):

  805.         pam_file = self.rs_path + ".aux.xml"

  806.         if os.path.isfile(pam_file):

  807.             os.remove(pam_file)

  808. 

  809.     def test_read_mode_error(self):

  810.         # Open raster in read mode

  811.         rs = GDALRaster(self.rs_path, write=False)

  812.         band = rs.bands[0]

  813.         self.addCleanup(self._remove_aux_file)

  814. 

  815.         # Setting attributes in write mode raises exception in the _flush method

  816.         with self.assertRaises(GDALException):

  817.             setattr(band, "nodata_value", 10)

  818. 

  819.     def test_band_data_setters(self):

  820.         # Create in-memory raster and get band

  821.         rsmem = GDALRaster(

  822.             {

  823.                 "datatype": 1,

  824.                 "driver": "MEM",

  825.                 "name": "mem_rst",

  826.                 "width": 10,

  827.                 "height": 10,

  828.                 "nr_of_bands": 1,

  829.                 "srid": 4326,

  830.             }

  831.         )

  832.         bandmem = rsmem.bands[0]

  833. 

  834.         # Set nodata value

  835.         bandmem.nodata_value = 99

  836.         self.assertEqual(bandmem.nodata_value, 99)

  837. 

  838.         # Set data for entire dataset

  839.         bandmem.data(range(100))

  840.         if numpy:

  841.             numpy.testing.assert_equal(

  842.                 bandmem.data(), numpy.arange(100).reshape(10, 10)

  843.             )

  844.         else:

  845.             self.assertEqual(bandmem.data(), list(range(100)))

  846. 

  847.         # Prepare data for setting values in subsequent tests

  848.         block = list(range(100, 104))

  849.         packed_block = struct.pack("<" + "B B B B", *block)

  850. 

  851.         # Set data from list

  852.         bandmem.data(block, (1, 1), (2, 2))

  853.         result = bandmem.data(offset=(1, 1), size=(2, 2))

  854.         if numpy:

  855.             numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))

  856.         else:

  857.             self.assertEqual(result, block)

  858. 

  859.         # Set data from packed block

  860.         bandmem.data(packed_block, (1, 1), (2, 2))

  861.         result = bandmem.data(offset=(1, 1), size=(2, 2))

  862.         if numpy:

  863.             numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))

  864.         else:

  865.             self.assertEqual(result, block)

  866. 

  867.         # Set data from bytes

  868.         bandmem.data(bytes(packed_block), (1, 1), (2, 2))

  869.         result = bandmem.data(offset=(1, 1), size=(2, 2))

  870.         if numpy:

  871.             numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))

  872.         else:

  873.             self.assertEqual(result, block)

  874. 

  875.         # Set data from bytearray

  876.         bandmem.data(bytearray(packed_block), (1, 1), (2, 2))

  877.         result = bandmem.data(offset=(1, 1), size=(2, 2))

  878.         if numpy:

  879.             numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))

  880.         else:

  881.             self.assertEqual(result, block)

  882. 

  883.         # Set data from memoryview

  884.         bandmem.data(memoryview(packed_block), (1, 1), (2, 2))

  885.         result = bandmem.data(offset=(1, 1), size=(2, 2))

  886.         if numpy:

  887.             numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))

  888.         else:

  889.             self.assertEqual(result, block)

  890. 

  891.         # Set data from numpy array

  892.         if numpy:

  893.             bandmem.data(numpy.array(block, dtype="int8").reshape(2, 2), (1, 1), (2, 2))

  894.             numpy.testing.assert_equal(

  895.                 bandmem.data(offset=(1, 1), size=(2, 2)),

  896.                 numpy.array(block).reshape(2, 2),

  897.             )

  898. 

  899.         # Test json input data

  900.         rsmemjson = GDALRaster(JSON_RASTER)

  901.         bandmemjson = rsmemjson.bands[0]

  902.         if numpy:

  903.             numpy.testing.assert_equal(

  904.                 bandmemjson.data(), numpy.array(range(25)).reshape(5, 5)

  905.             )

  906.         else:

  907.             self.assertEqual(bandmemjson.data(), list(range(25)))

  908. 

  909.     def test_band_statistics_automatic_refresh(self):

  910.         rsmem = GDALRaster(

  911.             {

  912.                 "srid": 4326,

  913.                 "width": 2,

  914.                 "height": 2,

  915.                 "bands": [{"data": [0] * 4, "nodata_value": 99}],

  916.             }

  917.         )

  918.         band = rsmem.bands[0]

  919.         # Populate statistics cache

  920.         self.assertEqual(band.statistics(), (0, 0, 0, 0))

  921.         # Change data

  922.         band.data([1, 1, 0, 0])

  923.         # Statistics are properly updated

  924.         self.assertEqual(band.statistics(), (0.0, 1.0, 0.5, 0.5))

  925.         # Change nodata_value

  926.         band.nodata_value = 0

  927.         # Statistics are properly updated

  928.         self.assertEqual(band.statistics(), (1.0, 1.0, 1.0, 0.0))

  929. 

  930.     def test_band_statistics_empty_band(self):

  931.         rsmem = GDALRaster(

  932.             {

  933.                 "srid": 4326,

  934.                 "width": 1,

  935.                 "height": 1,

  936.                 "bands": [{"data": [0], "nodata_value": 0}],

  937.             }

  938.         )

  939.         self.assertEqual(rsmem.bands[0].statistics(), (None, None, None, None))

  940. 

  941.     def test_band_delete_nodata(self):

  942.         rsmem = GDALRaster(

  943.             {

  944.                 "srid": 4326,

  945.                 "width": 1,

  946.                 "height": 1,

  947.                 "bands": [{"data": [0], "nodata_value": 1}],

  948.             }

  949.         )

  950.         rsmem.bands[0].nodata_value = None

  951.         self.assertIsNone(rsmem.bands[0].nodata_value)

  952. 

  953.     def test_band_data_replication(self):

  954.         band = GDALRaster(

  955.             {

  956.                 "srid": 4326,

  957.                 "width": 3,

  958.                 "height": 3,

  959.                 "bands": [{"data": range(10, 19), "nodata_value": 0}],

  960.             }

  961.         ).bands[0]

  962. 

  963.         # Variations for input (data, shape, expected result).

  964.         combos = (

  965.             ([1], (1, 1), [1] * 9),

  966.             (range(3), (1, 3), [0, 0, 0, 1, 1, 1, 2, 2, 2]),

  967.             (range(3), (3, 1), [0, 1, 2, 0, 1, 2, 0, 1, 2]),

  968.         )

  969.         for combo in combos:

  970.             band.data(combo[0], shape=combo[1])

  971.             if numpy:

  972.                 numpy.testing.assert_equal(

  973.                     band.data(), numpy.array(combo[2]).reshape(3, 3)

  974.                 )

  975.             else:

  976.                 self.assertEqual(band.data(), list(combo[2]))
About Sign Up Leave Feedback Contact