Source code for flusstools.geotools.geotools

"""``geotools`` is a package for creating, modifying, and transforming geospatial datasets."""


from .kml import *
from .srs_mgmt import *
gdal.UseExceptions()


[docs]def float2int(raster_file_name, band_number=1): """Converts a float number raster to an integer raster (required for converting a raster to a polygon shapefile). Args: raster_file_name (str): Target file name, including directory; must end on ``".tif"``. band_number (int): The raster band number to open (default: ``1``). Returns: str: ``"path/to/ew_raster_file.tif"`` """ raster, array, geo_transform = raster2array( raster_file_name, band_number=band_number) try: array = array.astype(int) except ValueError: logging.error("! Invalid raster pixel values.") return raster_file_name new_name = raster_file_name.split(".tif")[0] + "_int.tif" # get source coordinate system and exit function if not possible src_srs = get_srs(raster) if not src_srs: # ensure consistency return raster_file_name # create integer raster print(" * info: creating integer raster to Polygonize:\n >> %s" % new_name) create_raster(new_name, array, epsg=int(src_srs.GetAuthorityCode(None)), rdtype=gdal.GDT_Int32, geo_info=geo_transform) return new_name
[docs]def raster2line(raster_file_name, out_shp_fn, pixel_value, max_distance_method="simplified"): """Converts a raster to a line shapefile, where ``pixel_value`` determines line start and end points. Args: raster_file_name (str): of input raster file name, including directory; must end on ``".tif"``. out_shp_fn (str): of target shapefile name, including directory; must end on ``".shp"``. pixel_value (``int`` or ``float``): Pixel values to connect. max_distance_method (str): change to (pixel) ``"width"`` or ``"height"`` to force lines to exactly follow pixels (no triangulation). Returns: None: Writes a new shapefile to disk. """ # calculate max. distance between points # ensures correct neighbourhoods for start and end pts of lines raster, array, geo_transform = raster2array(raster_file_name) pixel_width = geo_transform[1] # verify if user provided max_distance_method argument if not max_distance_method == "simplified": if "height" in max_distance_method: max_distance = geo_transform[1] else: # assume that user want pixel width max_distance = pixel_width else: max_distance = np.ceil(np.sqrt(2 * pixel_width**2)) # extract pixels with the user-defined pixel value from the raster array trajectory = np.where(array == pixel_value) if np.count_nonzero(trajectory) == 0: logging.error( "! The defined pixel_value (%s) does not occur in the raster band." % str(pixel_value)) return None # convert pixel offset to coordinates and append to nested list of points points = [] count = 0 for offset_y in trajectory[0]: offset_x = trajectory[1][count] points.append(offset2coords(geo_transform, offset_x, offset_y)) count += 1 # create multiline (write points dictionary to line geometry (wkbMultiLineString) multi_line = ogr.Geometry(ogr.wkbMultiLineString) for i in itertools.combinations(points, 2): point1 = ogr.Geometry(ogr.wkbPoint) point1.AddPoint(i[0][0], i[0][1]) point2 = ogr.Geometry(ogr.wkbPoint) point2.AddPoint(i[1][0], i[1][1]) distance = point1.Distance(point2) if distance < max_distance: line = ogr.Geometry(ogr.wkbLineString) line.AddPoint(i[0][0], i[0][1]) line.AddPoint(i[1][0], i[1][1]) multi_line.AddGeometry(line) # write multiline (wkbMultiLineString2shp) to shapefile new_shp = create_shp( out_shp_fn, layer_name="raster_pts", layer_type="line") lyr = new_shp.GetLayer() feature_def = lyr.GetLayerDefn() new_line_feat = ogr.Feature(feature_def) new_line_feat.SetGeometry(multi_line) lyr.CreateFeature(new_line_feat) # create projection file srs = get_srs(raster) make_prj(out_shp_fn, int(srs.GetAuthorityCode(None))) print(" * success (raster2line): wrote %s" % str(out_shp_fn))
[docs]def raster2polygon(file_name, out_shp_fn, band_number=1, field_name="values"): """Converts a raster to a polygon shapefile. Args: file_name (str): Target file name, including directory; must end on ``".tif"`` out_shp_fn (str): Shapefile name (with directory e.g., ``"C:/temp/poly.shp"``) band_number (int): Raster band number to open (default: ``1``) field_name (str): Field name where raster pixel values will be stored (default: ``"values"``) add_area (bool): If ``True``, an "area" field will be added, where the area in the shapefiles unit system is calculated (default: ``False``) Returns: osgeo.ogr.DataSource: Python object of the provided ``out_shp_fn``. """ logging.info(" * Polygonizing %s ..." % str(file_name)) # ensure that the input raster contains integer values only and open the input raster file_name = float2int(file_name) raster, raster_band = open_raster(file_name, band_number=band_number) # create new shapefile with the create_shp function new_shp = create_shp( out_shp_fn, layer_name="raster_data", layer_type="polygon") dst_layer = new_shp.GetLayer() # create new field to define values new_field = ogr.FieldDefn(field_name, ogr.OFTInteger) dst_layer.CreateField(new_field) # Polygonize(band, hMaskBand[optional]=None, destination lyr, field ID, papszOptions=[], callback=None) gdal.Polygonize(raster_band, None, dst_layer, 0, [], callback=None) # create projection file srs = get_srs(raster) make_prj(out_shp_fn, int(srs.GetAuthorityCode(None))) logging.info(" * success (Polygonize): wrote %s" % str(out_shp_fn)) return new_shp
[docs]def rasterize(in_shp_file_name, out_raster_file_name, pixel_size=10, no_data_value=-9999, rdtype=gdal.GDT_Float32, overwrite=True, interpolate_gap_pixels=False, **kwargs): """Converts any ESRI shapefile to a raster. Args: in_shp_file_name (str): A shapefile name (with directory e.g., ``"C:/temp/poly.shp"``) out_raster_file_name (str): Target file name, including directory; must end on ``".tif"`` pixel_size (float): Pixel size as multiple of length units defined in the spatial reference (default: ``10``) no_data_value (int OR float): Numeric value for no-data pixels (default: ``-9999``) rdtype (gdal.GDALDataType): The raster data type (default: ``gdal.GDT_Float32`` (32 bit floating point) overwrite (bool): Overwrite existing files (default: ``True``) interpolate_gap_pixels (bool): Fill empty pixels that are not touched by a shapefile element with interpolated values (default: ``False``) Keyword Args: field_name (str): Name of the shapefile's field with values to burn to raster pixel values. radius1 (float): Define the x-radius for interpolating pixels (default: ``-1``, corresponding to infinity). Only applicable ``with interpolate_gap_pixels``. radius2 (float): Define the y-radius for interpolating pixels (default: ``-1``, corresponding to infinity). Only applicable ``with interpolate_gap_pixels``. power (float): Power of the function for interpolating pixel values (default: ``1.0``, corresponding to linear). smoothing (float): Smoothing parameter for interpolating pixel values (default: ``0.0``). min_points (int): Minimum number of points to use for interpolation. If the interpolator cannot find at least ``min_points`` for a pixel, it assigns a ``no_data`` value to that pixel (default: ``0``). max_points (int): Maximum number of points to use for interpolation. The interpolator will not use more than ``max_points`` closest points to interpolate a pixel value (default: ``0``). Hints: More information on pixel value interpolation: ``interpolate_gap_pixels=True`` interpolates values at pixels that are not touched by any las point. The pixel value interpolation uses ``gdal_grid`` (i.e., its Python bindings through ``gdal.Grid()``). Control the interpolation parameters with the keyword arguments ``radius1``, ``radius2``, ``power``, ``max_points``, ``min_points``, and ``smoothing``.. Returns: int: Creates the GeoTIFF raster defined with ``out_raster_file_name`` (success: ``0``, otherwise ``None``). """ default_keys = {"radius1": -1, "radius2": -1, "power": 1.0, "smoothing": 0.0, "min_points": 0, "max_points": 0, } for k in default_keys.keys(): if kwargs.get(k): default_keys[k] = str(kwargs.get(k)) # check if any action is required if os.path.isfile(out_raster_file_name) and not overwrite: logging.info(" * %s already exists. Nothing to do." % out_raster_file_name) return None # open data source try: source_ds = ogr.Open(in_shp_file_name) except RuntimeError as err: logging.error("! Could not open %s." % str(in_shp_file_name)) return None source_lyr = source_ds.GetLayer() # read extent x_min, x_max, y_min, y_max = source_lyr.GetExtent() # get x and y resolution in number of pixel x_res = int((x_max - x_min) / pixel_size) y_res = int((y_max - y_min) / pixel_size) # get spatial reference system and assign to raster srs = get_srs(source_ds) try: srs.ImportFromEPSG(int(srs.GetAuthorityCode(None))) except RuntimeError as err: logging.error(e) return None if float(pixel_size) < 1.0: logging.info( " -- Yeek! This will be a high resolution raster. Be prepared that your system resources will be occupied for a while.") # use gdal.Grid if gap interpolation (fill void pixels) is True if interpolate_gap_pixels: logging.info( " * Creating gridded raster with interpolated values for empty pixels from neighbouring pixels ...") logging.info( " -- Note: to deactivate pixel value interpolation option use interpolate_gap_pixels=False") try: algorithm = "invdist:power={0}:radius1={1}:radius2={2}:smoothing={3}:min_points={4}:max_points={5}".format( str(default_keys["power"]), str( default_keys["radius1"]), str(default_keys["radius2"]), str(default_keys["smoothing"]), str( default_keys["min_points"]), str(default_keys["max_points"]) ) gdal.Grid(out_raster_file_name, in_shp_file_name, algorithm=algorithm, zfield=kwargs.get("field_name"), outputType=rdtype, outputSRS=srs, width=x_res, height=y_res, outputBounds=[x_min, y_min, x_max, y_max]) return 0 except KeyError: logging.error("! Invalid gdal.Grid options provided.") return None except RuntimeError as err: logging.error("! %s." % str(err)) # create destination data source (GeoTIff raster) try: target_ds = gdal.GetDriverByName('GTiff').Create( out_raster_file_name, x_res, y_res, 1, eType=rdtype) except RuntimeError as err: logging.error("! Could not create %s." % str(out_raster_file_name)) return None target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size)) band = target_ds.GetRasterBand(1) band.Fill(no_data_value) band.SetNoDataValue(no_data_value) # assign spatial reference target_ds.SetProjection(srs.ExportToWkt()) # RasterizeLayer(Dataset dataset, int bands, Layer layer, pfnTransformer=None, pTransformArg=None, # int burn_values=0, options=None, GDALProgressFunc callback=0, callback_data=None) try: if kwargs.get("field_name"): gdal.RasterizeLayer(target_ds, [1], source_lyr, None, None, burn_values=[0], options=["ALL_TOUCHED=TRUE", "ATTRIBUTE=" + str(kwargs.get("field_name"))]) else: gdal.RasterizeLayer(target_ds, [1], source_lyr, None, None, burn_values=[0], options=["ALL_TOUCHED=TRUE"]) except RuntimeError as err: logging.error("! Could not rasterize (burn values from %s)." % str(in_shp_file_name)) return None # release raster band band.FlushCache()