Geospatial application in python using OOP

Question

I'm working on a tool that performs some computations using elevations from a DTM file retrieved from a WCS service. These computations will be done on predefined regions of Europe. I'm trying to encapsulate some things in classes, but it is my first time really using them, so I would love to get some feedback. Basically I created a Region class that contains several attributes and internal methods that connects to the WCS service and handles the response. Since the regions are predefined, I also created 2 additional classes that store constants: BboxFlood, that contains the BBOX I want to retrieve from the service (and will be used to do some flooding operations) and Borders, which contains the outline of the region (which is the real region of interest). This outline is used to create a mask (masked_raster) which I will modify using other functions lately.

import matplotlib.pyplot as plt
import numpy as np
    
from requests.exceptions import MissingSchema
from owslib.wcs import WebCoverageService
from owslib.util import ServiceException
from rasterio.io import MemoryFile
from rasterio.errors import RasterioIOError
from rasterio.plot import show
from geopandas import read_file
from rasterio.features import geometry_mask
from shapely.geometry import mapping
from fiona.errors import DriverError


class Region():
    
    
    def __init__(self, name, bbox, epsg, outline):
        """
        Create a region of interest

        Parameters
        ----------
        name : str
            Region name.
        bbox : tuple or list
            Bounding box of the region
        epsg : str
            Spatial reference system
        outline : shapely.geometry.polygon.Polygon or str, optional
            A shapely Polygon defining the border of the region of interest.
            Alternatively, a path to the shapefile. 
            NOTE: The outline epsg should coincide with the raster one.

        Returns
        -------
        None.

        """
            
        self.name = name
        self.bbox = bbox
        self.epsg = epsg
        self.outline = outline
        
        self.raster = None
        self.masked_raster = None
        self.mask = None
        
        self.profile = None
        self.transform = None
        self.resx = None
        self.resy = None
        
        if isinstance(outline, str):
            self._read_outline(outline)
            
         
    def __str__(self):
        return f"Region: {self.name}, Bbox: {self.bbox}, EPSG: {self.epsg}"
    
    
    def __repr__(self):
        return f"Region({self.name!r}, {self.bbox!r}, {self.epsg!r})"
        
    
    def _get_bytes_from_response(self, wcs_response):
        """
        Get tif file from multipart response from DV

        Parameters
        ----------
        wcs_response : owslib.util.ResponseWrapper
            wcs response after getCoverage

        Returns
        -------
        image_bytes : bytes
            image file in bytes

        """

        boundary = b"--wcs"
        
        bytes_parts = wcs_response.read().split(boundary)
        
        for part in bytes_parts:
            
            if part.startswith(b"\nContent-Type: image/tiff"):
                index = part.index(b"II*")
                image_bytes = part[index:]
                break
        
        
        return image_bytes
    
    
    def _bytes_to_raster(self, response_as_bytes):
        # This internal method is used to read the response retrieved from 
        # the WCS service
                
        try:
            # Create a file in memory not to write any file to the hard disk
            with MemoryFile(response_as_bytes) as memfile:
                with memfile.open() as src:
                    raster = src.read(1)
                    profile = src.profile
                    res = src.res
                      
        except RasterioIOError() as e:
            print("Rasterio error:", e)
            
        except Exception as e:
            print('Unexpected error:', e)

        else:
            self.raster = raster
            self.profile = profile
            self.transform = profile['transform']
            self.resx = res[0]
            self.resy = res[1]
            
            
    def _read_outline(self, outline_path):
        # Read outline from path if not given as a Polygon
        try:
            gdf = read_file(outline_path)
     
        except DriverError as e:
            print("Incorrect file name:", e)
            
        else:
            self.outline = gdf.geometry.iloc[0]
    
    
    def _set_masked_raster(self):
        # Create a masked version of the raster based on outline
        mask = geometry_mask([mapping(self.outline)],
                             out_shape=self.raster.shape,
                             transform=self.transform,
                             invert=True
                             )

        self.mask = ~mask
        self.masked_raster = np.where(mask, 
                                      self.raster, 
                                      self.profile['nodata'])
        

    def get_dtm_from_wcs_201(self, 
                         raster_source='https://geo.api.vlaanderen.be/DHMV/wcs',
                         layer='DHMVII_DTM_1m', 
                         format='image/tiff',
                         version="2.0.1",
                         scalefactor=3):
        """
        Get raster and metadata from a WCS service

        Parameters
        ----------
        raster_source : str, optional
            URL of WCS service.
        layer : str, optional
            CoverageID.
        format : str, optional
            File type of the request. The default is 'image/tiff'.
        version : str, optional
            Version of the WCS protocol. The default is "2.0.1".
        scalefactor : int, optional
            Scale. The default is 3.

        Returns
        -------
        None.

        """      
        
        sset = [('x', self.bbox[0], self.bbox[2]), ('y', self.bbox[1], self.bbox[3])]
        
        try: 
            # Connect to WCS and try getting a coverage
            wcs = WebCoverageService(raster_source, version=version)
            response = wcs.getCoverage(identifier=[layer],
                                       bbox=self.bbox,
                                       crs=self.epsg,
                                       subsets=sset,
                                       format=format,
                                       scalefactor=scalefactor, 
                                       )
            
        except MissingSchema as e:
            print('Invalid URL:', e)
            
        except ServiceException as e:
            print("Service error:", e)
            
        except Exception as e:
            print("An unexpected error has ocurred", e)
        
        else:
            image_bytes = self._get_bytes_from_response(response)
            self._bytes_to_raster(image_bytes)
            self._set_masked_raster()
                
    
    def reset_masked_raster(self):
        self._set_masked_raster()



from shapely.geometry import Polygon

class BboxFlood():
   
    ROMBOUT = {'bbox': (71274, 213729, 75273, 219336), 
                     'epsg': 'EPSG:31370'}

class Borders():
    
    ROMBOUT = {'outline': Polygon([(71275.24244585728, 214206.7075793226),
                                    (71332.89192388157, 214301.1237481891),
                                    (71451.59618105175, 214522.12153573174),
                                    (71542.7001391952, 214700.03943770562),
                                    (71605.07543163445, 214832.82320418613),
                                    (71750.3815106122, 215099.33581733558),
                                    (71788.18471815114, 215162.6561899633),
                                    (71811.33918276873, 215193.37129608868),
                                    (71853.86779125003, 215234.00974419303),
                                    (71886.47305775236, 215265.19739041265),
                                    (71909.15498227572, 215310.56123945938),
                                    (71914.82546340655, 215363.9582701081),
                                    (71920.02340444316, 215416.88276066261),
                                    (71928.05658604519, 215456.57612857848),
                                    (71937.50738792992, 215485.40107432692),
                                    (71965.38725348988, 215542.1058856353),
                                    (72009.33348225389, 215614.87706014776),
                                    (72072.18131478736, 215710.8026992778),
                                    (72097.69847987614, 215755.69400823026),
                                    (72152.04059071335, 215863.43314971623),
                                    (72178.03029589636, 215913.99493979954),
                                    (72187.95363787533, 215954.6333879039),
                                    (72195.04173928888, 215999.05215676213),
                                    (72187.4810977811, 216064.2626897668),
                                    (72169.99711429434, 216140.81418503314),
                                    (72168.57949401163, 216214.0578996398),
                                    (72180.86553646179, 216347.31420621453),
                                    (72189.84379825227, 216402.128857146),
                                    (72197.8769798543, 216429.53618261172),
                                    (72221.97652466036, 216482.93321326046),
                                    (72264.26886309455, 216559.01216843256),
                                    (72300.18191025654, 216626.58540190838),
                                    (72334.20479704157, 216687.54307406492),
                                    (72368.22768382661, 216733.8520033001),
                                    (72415.95423334451, 216787.72157404307),
                                    (72501.95653049558, 216866.16322968635),
                                    (72556.77118142703, 216924.28566127748),
                                    (72574.25516491377, 216948.8577461778),
                                    (73971.55622357152, 219165.54332824182),
                                    (74158.68210088921, 219335.65776216701),
                                    (74694.54256775357, 218468.07414914857),
                                    (74856.1512799825, 218306.46543691962),
                                    (75272.9316430992, 217991.75373415803),
                                    (72559.60642199244, 214402.33917833655),
                                    (71726.04569575902, 213730.38716433206),
                                    (71275.24244585728, 214206.7075793226)]),
                 'epsg':'EPSG:31370'}


# Example

rombout = Region("Romb", **BboxFlood.ROMBOUT, 
                  outline=Borders.ROMBOUT['outline'])
    
rombout.get_dtm_from_wcs_201(scalefactor=5)

I have another module that contains functions that takes a region as parameter (they modify the masked_raster attribute) What do you think? Is it well-structured? One of my biggest concern is how to test it (I have been doing it manually). Any feedback is welcomed.

Answer 1

Pep-8 asks that you organize your imports, perhaps with $ isort *.py. That delayed Polygon import is certainly a surprise.

class documentation

class Region():

You passed up the opportunity to write a """docstring""" telling us the single responsibility of this class. It comes in handy a month down the road when you're tempted to add a kitchen_sink attribute or method, and then upon reading the docstring you realize the new thing belongs in a new class that has a different responsibility.

Dunno what those empty () brackets are doing there. Makes it look like we have a python2 "new style" class that inherits from object. Just elide them.

signature should match docs

You assert that outline is "optional". No, it's not. Caller will get a TypeError if it's missing.

Recommend you give it a keyword default of None.

The epsg description is helpful, thank you. But a naïve human would find it more helpful if "WGS84" was mentioned as a valid value, so we know that uppercase is OK and that "4326" is not what it's looking for.

Even better, leave the docstring alone and give the signature a keyword default of "WGS84".

    def __init__(self, name, bbox, epsg, outline):
        """
        Create a region of interest

        Parameters
        ----------
        name : str
            Region name.
        bbox : tuple or list
            Bounding box of the region
        epsg : str
            Spatial reference system
        outline : shapely.geometry.polygon.Polygon or str, optional

You can tell humans about types in that way. But it's better to tell humans + mypy:

    def __init__(self, name: str, bbox: tuple | list,
                 epsg: str, outline: Polygon | str | None):

Alternatively, a path to the shapefile.

It makes me sad that you chose to use generic str rather than Path to model the concept of a pathname.

Returns ... None.

We can end the signature annotation with -> None:. But I usually elide that on __init__ ctors, since it's obvious to both man and machine. It's certainly helpful on other methods and on def main() -> None:. No need for "None" in the docstring.

I do thank you for the "table of contents", where you assign raster through resy attributes.

tuple vs list

        bbox : tuple or list

That should be just "tuple", please.

In python we use a list for "zero-or-more of same thing", and a tuple for "fixed number of different things" which can have different types. A tuple is like a C struct. Here, bbox[0] means something entirely different from bbox[3] -- the meaning comes from the position. And it's not like a list of five numbers would be valid.

Even nicer would be to define a Bbox named tuple. Then instead of cryptic [3] subscripts we could use (abbreviated) names for the elements.

assign from helper

The ctor ends with

        if isinstance(outline, str):
            self._read_outline(outline)

Prefer

        if isinstance(outline, str):
            self.outline = self._read_outline(outline)

and then the helper ends with return gdf.geometry.iloc[0].

It makes the side effect obvious, rather than hoping the reader will infer a side effect from the spelling of the helper's name.

errors should be fatal

    def _read_outline(self, outline_path):
        ...     
        except DriverError as e:
            print("Incorrect file name:", e)

Either don't catch DriverError, or wrap e and re-raise it with that helpful "incorrect..." diagnostic. As it stands we print() and move on as if self.outline contained something useful. That's just kicking the can down the road, so some other method will get blamed for a fatal error. Accept the blame here, to expedite diagnosis and bug fixes.

Also, consider making f"{outline_path}" part of the diagnostic, for quicker diagnosis.

meaningful identifiers

    def _get_bytes_from_response(self, wcs_response):
        """
        Get tif file from multipart response from DV

This is pretty good.

Much of the docstring would be superfluous if it was named

    def _get_image_from_coverage( ... ) -> bytes: 

I don't know what "DV" denotes, but that's ok.

Consider using the variant spelling of "tiff", or "TIFF", given that we're not talking about a *.tif file extension here.

        bytes_parts = wcs_response.read().split(boundary)

It would be enough to just call this parts. Add an annotation if you feel that bytes isn't clear to the reader.

            if part.startswith(b"\nContent-Type: image/tiff"):

This is on the fragile side. Yes, I understand that WCS uses this spelling ATM. The header case could change yet it would mean the same thing. Consider case smashing to .lower() before the comparison.

docstring rather than comment

This is good.

    def _bytes_to_raster( ... ):
        # This internal method is used to read the response retrieved from 
        # the WCS service

It would be better if it appeared in """docstring""" form, just a single English sentence. (And then add a type annotation, please.)

        except RasterioIOError() as e:
            print("Rasterio error:", e)
            
        except Exception as e:
            print('Unexpected error:', e)

No, no, please just let those bubble up the call stack as fatal. As it stands, you're swallowing the exception, carrying on, and not even showing a maintenance engineer the source line number in a stack trace.

a new dataclass

This is kind of interesting:

        else:
            self.raster = raster
            self.profile = profile
            self.transform = profile['transform']
            self.resx = res[0]
            self.resy = res[1]

There's an atomic all-or-nothing aspect to this, as if you were concerned de-referencing the src.profile or src.res properties could raise.

Up in the ctor I was thinking "gosh, seems like there's a bunch of attributes on this object." This stanza of five assignments suggests that maybe we should hang onto a single self.raster object, which in turn has five attributes. The @dataclass decorator is your friend.

The (mask, masked_raster) tuple might be viewed in a similar light.

good defaults

I found get_dtm_from_wcs_201() extremely clear, thank you for those helpful default values. The "dtm" part is obscure, so the only part of the docstring I really need is where it explains that a DTM is apparently raster + metadata.

Also layer is perhaps really a layer_id.

Usual grumble: things like MissingSchema and ServiceException should re-raise as fatal.

new module

The ROMBOUT data dump should probably go in a new rombout.py module.

"main" guard

rombout.get_dtm_from_wcs_201(scalefactor=5)

Please use the if __name__ == "__main__": idiom, so e.g. a unit test could safely import this.

---

This code achieves its design goals.

I would be willing to delegate or accept maintenance tasks on it.

---

EDIT

I read somewhere that for private methods commenting instead of docstring was the standard.

I suspect someone had tooling that would accidentally add _my_helper to published ReadTheDocs web pages if it had a docstring.

Most folks use sphinx, which by default won't do that. I've never seen that it's a concern.

Every method, whether it's part of your Public API or not, has a specification, a contract. It can be poorly specified and implicit. Or you can use a clear English sentence to spell it out in a docstring. We want to know its single responsibility. We want to know what it promises will be true upon return. A docstring is an excellent way to communicate that to your colleagues. For public() and for _private() methods.

The DTM is a digital terrain model. I will change.

No biggie. I mean, I knew that it was some Term of Art, and I knew that I wasn't familiar enough with the literature but probably a maintenance engineer a few months from now would be. So I wasn't worried. I was just offering that a # comment which expands the TLA would have been welcome, is all. And sure, getting elevation would be a helpful phrase in the identifier.

Use "tiff", or "TIFF",

It was just a spelling nit I was quibbling over. We say *.tif in a comment when talking about file extensions. But we say "tiff", or more properly "TIFF", when we mean Tagged Image File Format. (There's no such thing as TIF.)

I didn't understand the dataclass part neither. You meant making other classes?

Yes.

In OO design we look for nouns, "objects", that exist in the Business Domain and that we can model within the machine. Now, a noun could be "big", like a Kitchen and all it contains. Or our modeling could be more granular, and we might model a Stove, a Microwave, and so on.

Your Region has eleven attributes, which isn't necessarily Too Many, but it's a code smell, it's starting to look a little on the big side. It makes you wonder whether there could be some grouping of things, which collectively are their own thing.

So when I got down to that atomic assignment of five attributes I said "Aha!", I had been looking for such an opportunity, and maybe I just found it. Now, it's up to you. I don't know what makes the most sense in your environment. But I offered it as something to consider -- maybe you'd like to group them together. And similarly for that masking 2-tuple.

assignment vs side-effect

If a ctor is creating a House object, I would much rather see

        self.kitchen = self._get_kitchen( ... )

than

        self._build_kitchen( ... )

Now, maybe the "side effect" aspect of that helper is apparent from the name. But I worry that it might update self.plumbing and some other things. Plus, it's hard to unit test. If I can reveal the side effecting via an assignment like that, it will be my preference.

Answer 2

manifesto

Code should be bug-free.

Design by Contract is one path toward that. It promotes composability of components.

What should I do with the fatal errors? Make the program crash?

Yes.

Or make it do the Right Thing, if you prefer. (The current implementation just ignores the error and stumbles forward.)

total function

A function is specified to accept inputs and produce some result: return value and/or side effects.

The inputs can be explicit in the signature, they can also be global variables, filesystem fullness, wallclock time, PRNG state, and so on.

Some functions are total -- they will always produce a correct result, for any input. Here is one, in python:

def increment(n: int) -> int:
    return n + 1

Assuming we can malloc() and the power is on, this function pretty much always returns what you expect, when given the proper object type.

Here is a similar function in C:

uchar increment(uchar n) {
    return n + 1;
}

Calling increment(254) does the Right Thing. But probably the caller is being too adventurous upon supplying an unsigned arg of 255. For modulo operations this is well defined, but if the contract is to "return a bigger number", then this is a partial function. If this were C++ we might choose to throw an exception for that.

The familiar sqrt() function offers a similar example. If we rule out NaN and complex / imaginary numbers, then caller should avoid offering a negative argument. This is called a pre-condition on acceptable inputs. If we do see negative input, it is appropriate for sqrt to immediately raise a fatal error so caller will fix their broken code.

post-condition

The contract for the root finder looks like this: "Give me a non-negative real. I promise to hand back a real which, when squared, equals the original argument." (We are in FP-land, here, so we finesse "equals" to mean "low relative error", typically specified by some epsilon parameter.)

Now, in a case like sqrt(-3), there is simply no way to fulfill that contract. What should the root finder do? Return a sentinel like 99? (Nevermind that it's indistinguishable from sqrt(9801).) Maybe "better" sentinel values happen to be available, like -99? Sure, that's one possibility that folks have tried. It makes the contract more complex, and some callers will forget to check for that condition. Or, we could throw up our hands!

The much better alternative, when faced with an "impossible" contract, is to give up, with raise. Did we return some float value? No. We returned nothing, the computation didn't happen.

Throwing an exception interrupts the flow so caller won't accidentally trust some arbitrary sentinel value as being "correct". It might be a fatal exception, which halts with a stack trace. Or the caller might have anticipated this situation, and have some strategy for coping with it. Caller might try an alternative call (flip the sign), or pause while it asks a human for help, or re-raise so its higher level caller is aware of the situation.

Make the program crash?

Yes. Or no. It's not your responsibility. A pre-condition wasn't satisfied or you can't make your post-condition true for whatever reason, so you bail out and make the caller aware of it. What the caller does is up to him, it's not your job.