The number of pixels to expand around a center pixel
* A 'size' of 0 would make a processing matrix with size 1x1. As in,
just the value at each point. This would be silly to call...
* A 'size' of 1 would make a processing matrix of size 3x3. As in, one
pixel around the center pixel in all directions
* Processed matrix size is equal to 2*size+1
Whether or not to pass the x and y index to the processing function
* If True, the decorated function must accept an input called 'xi' and
'yi' in addition to the matrix
* The xi and yi correspond to the index of the center pixel in the
original matrix
Returns:
function
–
Example:
Say we want to make a processor which calculates the average of pixels
which are within a distance of 2 indices. In other words, we want the
average of a 5x5 matrix centered around each pixel.
Assume that we can use the value -9999 as a no data value
@KernelProcessor(2, edgeValue=-9999)
def getMean( mat ):
# Get only good values
goodValues = mat[mat!=-9999]
defKernelProcessor(size,edgeValue=0,outputType=None,passIndex=False):"""A decorator which automates the production of kernel processors for use in mutateRaster (although it could really used for processing any matrix). Parameters ---------- size : int The number of pixels to expand around a center pixel * A 'size' of 0 would make a processing matrix with size 1x1. As in, just the value at each point. This would be silly to call... * A 'size' of 1 would make a processing matrix of size 3x3. As in, one pixel around the center pixel in all directions * Processed matrix size is equal to 2*size+1 edgeValue : numeric; optional The value to apply to the edges of the matrix before applying the kernel * Will be factored into the kernelling when processing near the edges outputType : np.dtype; optional The datatype of the processed values * Only useful if the output type of the kerneling step is different from the matrix input type passIndex : bool Whether or not to pass the x and y index to the processing function * If True, the decorated function must accept an input called 'xi' and 'yi' in addition to the matrix * The xi and yi correspond to the index of the center pixel in the original matrix Returns ------- function Example: -------- * Say we want to make a processor which calculates the average of pixels which are within a distance of 2 indices. In other words, we want the average of a 5x5 matrix centered around each pixel. * Assume that we can use the value -9999 as a no data value >>> @KernelProcessor(2, edgeValue=-9999) >>> def getMean( mat ): >>> # Get only good values >>> goodValues = mat[mat!=-9999] >>> >>> # Return the mean >>> return goodValues.mean() """passdefwrapper1(kernel):passdefwrapper2(matrix):# get the original matrix sizesyN,xN=matrix.shape# make a padded version of the matrix# paddedMatrix = (# np.ones((yN + 2 * size, xN + 2 * size),) * edgeValue# )paddedMatrix=np.full(shape=(yN+2*size,xN+2*size),fill_value=edgeValue)paddedMatrix[size:-size,size:-size]=matrix# apply kernel to each pixeloutput=np.zeros((yN,xN),dtype=matrix.dtypeifoutputTypeisNoneelseoutputType)foryiinrange(yN):forxiinrange(xN):slicedMatrix=paddedMatrix[yi:2*size+yi+1,xi:2*size+xi+1]ifpassIndex:output[yi,xi]=kernel(slicedMatrix,xi=xi,yi=yi)else:output[yi,xi]=kernel(slicedMatrix)# done!returnoutputreturnwrapper2returnwrapper1
This function applies a buffer to any geom, avoiding edge issues with geoms
near the SRS bounds. By shifting the geom to a zero longitude, geometry
distortions are avoided when the buffered geom exceeds the bounds (i.e.
antimeridian or latitudes of +/-90° e.g. in the case of EPSG:4326). Buffered
geom areas extending over the bounds can either be clipped off or shifted to
the respective "other end of the map". If the buffer is applied in a
different (e.g. metric) EPSG, latitudinal overlaps will always be clipped.
geom : osgeo.ogr.Geometry
Geometry to be buffered.
buffer : int, float
The buffer value to be applied to the geom, in unit of the SRS unless
'bufferInEPSG6933' is True, then always in meters.
srs : Anything acceptable to geokit.srs.loadSRS(); optional
Allows to specify an EPSG integer code or an osgeo.osr.SpatialReference
instance to define the SRS in which the buffer will be applied, then in
the unit of the specified EPSG. If e.g. 6933 is given, the buffer will
be applied in meters in a metric system. Passing "laea" (case-insensitive)
will generate a geometry-centered, metric Lambert Azimuthal Equal Area
system in which the buffer will be applied. By default False, i.e. the
original SRS of the geom will be used.
split : str, optional
'shift' : shift areas that exceed the antimeridian line to the other end (default)
'clip' : remove/clip polygon parts that exceed the antimeridian
'none' : do not split geoms at all that cross the antimeridian
defapplyBuffer(geom,buffer,srs=None,split="shift"):""" This function applies a buffer to any geom, avoiding edge issues with geoms near the SRS bounds. By shifting the geom to a zero longitude, geometry distortions are avoided when the buffered geom exceeds the bounds (i.e. antimeridian or latitudes of +/-90° e.g. in the case of EPSG:4326). Buffered geom areas extending over the bounds can either be clipped off or shifted to the respective "other end of the map". If the buffer is applied in a different (e.g. metric) EPSG, latitudinal overlaps will always be clipped. geom : osgeo.ogr.Geometry Geometry to be buffered. buffer : int, float The buffer value to be applied to the geom, in unit of the SRS unless 'bufferInEPSG6933' is True, then always in meters. srs : Anything acceptable to geokit.srs.loadSRS(); optional Allows to specify an EPSG integer code or an osgeo.osr.SpatialReference instance to define the SRS in which the buffer will be applied, then in the unit of the specified EPSG. If e.g. 6933 is given, the buffer will be applied in meters in a metric system. Passing "laea" (case-insensitive) will generate a geometry-centered, metric Lambert Azimuthal Equal Area system in which the buffer will be applied. By default False, i.e. the original SRS of the geom will be used. split : str, optional 'shift' : shift areas that exceed the antimeridian line to the other end (default) 'clip' : remove/clip polygon parts that exceed the antimeridian 'none' : do not split geoms at all that cross the antimeridian """# create copy and extract SRS_srs_orig=geom.GetSpatialReference()_geom=geom.Clone()ifsrsisnotNone:# generate own geom-centered LAEA or load SRStry:# allow LAEA or EPSG or SRS objectsrs=SRS.loadSRS(srs,geom=_geom)exceptExceptionase:raiseValueError(f"Invalid SRS '{srs}': {e}")# make sure the poles are far enough to allow the planned bufferbuffer_mtrs=buffer*srs.GetLinearUnits()north_srs=SRS.loadSRS("+proj=aeqd +lat_0=90 +lon_0=0 +datum=WGS84 +units=m +no_defs")north_dist=transform(_geom,toSRS=north_srs).Distance(point(0,0,srs=north_srs))south_srs=SRS.loadSRS("+proj=aeqd +lat_0=-90 +lon_0=0 +datum=WGS84 +units=m +no_defs")south_dist=transform(_geom,toSRS=south_srs).Distance(point(0,0,srs=south_srs))ifmin([north_dist,south_dist])<=buffer_mtrs:raiseGeoKitGeomError(f"buffered geometry would intersect with North or South Pole.")ifnotsrs.IsSame(_srs_orig):# convert geom to srs_geom=transform(_geom,toSRS=srs)else:srs=_srs_orig# apply buffer_geom_buf=_geom.Buffer(buffer)assert_geom_buf.IsValid(),f"buffered geom invalid after applying buffer."# make sure that the buffered geom was not already partially projected by 360° if crossing antimeridian# envelope diffs on either side must be equal to buffer with 1% toleranceassertnp.isclose(_geom.GetEnvelope()[0]-_geom_buf.GetEnvelope()[0],buffer,atol=0,rtol=0.01)andnp.isclose(_geom_buf.GetEnvelope()[1]-_geom.GetEnvelope()[1],buffer,atol=0,rtol=0.01),f"buffered geom envelope does not match unbuffered geom envelope plus buffers on both sides."# now retransform to original srs if neededifsrsisnotNoneandnotsrs.IsSame(_srs_orig):_geom_buf=transform(_geom_buf,toSRS=_srs_orig,revert360degProj=True)# now split if neededifnot(splitisNoneorisinstance(split,str)andsplit.upper()=="NONE"):_geom_buf=fixOutOfBoundsGeoms(_geom_buf,how=split)assert_geom_buf.IsValid(),f"buffered and re-transformed geom invalid after '{split}' operation"return_geom_buf
Applies an arbitrary function to the data frames provided.
Convenience function to apply geopandas methods to a geokit-style
dataframe with 'geom' column with osgeo.ogr.Geometry objects.
NOTE: All arguments besides **kwargs must be passed as positional
arguments.
geopandasMethod : str, executable
Geopandas method to apply to the dataframe, either str-formatted
method name or method as a callable function.
dfs : pd.DataFrames
One or multiple comma-separated pd.DataFrames with 'geom' column
with osgeo.ogr.Geometry objects. Will be passed to the geopandas
method as positional arguments, starting from the first position.
*kwargs
Will be passed on to the geopandas function.
defapplyGeopandasMethod(geopandasMethod,*dfs,**kwargs):""" Applies an arbitrary function to the data frames provided. Convenience function to apply geopandas methods to a geokit-style dataframe with 'geom' column with osgeo.ogr.Geometry objects. NOTE: All arguments besides **kwargs must be passed as positional arguments. geopandasMethod : str, executable Geopandas method to apply to the dataframe, either str-formatted method name or method as a callable function. *dfs : pd.DataFrames One or multiple comma-separated pd.DataFrames with 'geom' column with osgeo.ogr.Geometry objects. Will be passed to the geopandas method as positional arguments, starting from the first position. **kwargs Will be passed on to the geopandas function. """# load geopandas# try:# import geopandas as gpd# except:# raise ImportError(# "'geopandas' is required for geokit.vector.createGeoDataFrame() but is not installed in the current environment."# )# get the method as callableifcallable(geopandasMethod):# we have a callable function already, just make sure its gpdifnotgeopandasMethod.__module__.split(".")[0]=="geopandas":raiseAttributeError(f"'{geopandasMethod}' is not a callable method from geopandas!")else:# not an executable yet, get the executable function via its nameifnotisinstance(geopandasMethod,str):# must be a str formatted nameraiseTypeError(f"'geopandasMethod' must be str-formatted geopandas method name if not callable.")geopandasMethod=getattr(gpd,geopandasMethod,None)# create geodataframes from input geokit dfsassertall([isinstance(_df,pd.DataFrame)and"geom"in_df.columnsfor_dfindfs]),(f"positional *dfs arguments must be geokit-style pd.DataFrames with 'geom' column.")gdfs=[createGeoDataFrame(dfGeokit=_df)for_dfindfs]# now apply geopandas method onto gdfsgdf=geopandasMethod(*gdfs,**kwargs)# now convert the result back to geokit style and returndf=createDataFrameFromGeoDataFrame(gdf)returndf
defbox(*args,srs:srs_input=4326)->ogr.Geometry:"""Make an ogr polygon object from extents. Parameters ---------- *args : 4 numeric argument, or one tuple argument with 4 numerics The X_Min, Y_Min, X_Max and Y_Max bounds of the box to create srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the point to create * If not given, longitude/latitude is assumed * srs MUST be given as a keyword argument Returns ------- ogr.Geometry Example: -------- box(xMin, yMin, xMax, yMax [, srs]) box( (xMin, yMin, xMax, yMax) [, srs]) """iflen(args)==1:xMin,yMin,xMax,yMax=args[0]eliflen(args)==4:xMin,yMin,xMax,yMax=argselse:raiseGeoKitGeomError('Incorrect number positional inputs (only accepts 1 or 4). Did you mean to specify "srs="?')# make sure inputs are goodxMin=float(xMin)xMax=float(xMax)yMin=float(yMin)yMax=float(yMax)assertxMin<xMax,f"xMin must be less than xMax"assertyMin<yMax,f"yMin must be less than yMax"# make boxoutBox=ogr.Geometry(ogr.wkbPolygon)ring=ogr.Geometry(ogr.wkbLinearRing)forx,yin[(xMin,yMin),(xMax,yMin),(xMax,yMax),(xMin,yMax),(xMin,yMin)]:ring.AddPoint(x,y)outBox.AddGeometry(ring)ifnotsrsisNone:srs=SRS.loadSRS(srs)outBox.AssignSpatialReference(srs)returnoutBox
Load a Lambert-Azimuthal-Equal_Area spatial reference system (SRS) centered
on a given set of latitude and longitude coordinates. Alternatively, a geom
can be passed to center the LAEA on.
The region shape to center the LAEA in. If given, lat and lon must not
be given, instead they will be defined automatically as the coordinates
of the region centroid.
defcenteredLAEA(lon=None,lat=None,name="unnamed_m",geom=None):""" Load a Lambert-Azimuthal-Equal_Area spatial reference system (SRS) centered on a given set of latitude and longitude coordinates. Alternatively, a geom can be passed to center the LAEA on. Parameters ---------- lon : float The longitude of the projection's center. Required if no geom is given. lat : float The latitude of the projection's center. Required if no geom is given. geom: osgeo.ogr.Geometry The region shape to center the LAEA in. If given, lat and lon must not be given, instead they will be defined automatically as the coordinates of the region centroid. Returns ------- osr.SpatialReference """ifgeomisNone:assertisinstance(lat,numbers.Number)andisinstance(lon,numbers.Number),("If geom is not passed, lat and lon must be given as float values.")else:assertisinstance(geom,ogr.Geometry),"geom must be given as osgeo.ogr.Geometry class object if not None."assertlatisNoneandlonisNone,"If geom is given, lat and lon must not be given."# check if lat/lon can be used or if it must be extracted from geom firstifnotgeomisNone:# transform to EPSG:4326 in case it not already is lat/lon projectiongeom=GEOM.transform(geom,toSRS=4326)# extract lat/lon centroid coordinates to center LAEA uponlon=geom.Centroid().GetX()lat=geom.Centroid().GetY()srs=osr.SpatialReference()srs.ImportFromWkt('PROJCS["{}",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",{}],PARAMETER["longitude_of_center",{}],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'.format(name,float(lat),float(lon)))ifgdal.__version__>="3.0.0":srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)# assert that the srs is valid (may be invalid if e.g. wrong integer codes were passed)assertsrs.Validate()==0,f"Created srs is invalid."returnsrs
Filepath to output raster file. If it is an existing file, datasets will
be added to output. Recommended to create a new file every time though.
If None, no output dataset will be loaded or created on disk and output
dataset kept in memory only, by default None
If True, metadata of output dataset will be a combination of all input
rasters, by default False.
NOTE: In the case of multiple values for the metadata keys, the last
dataset metadata will take precedence.
defcombineSimilarRasters(datasets,output=None,combiningFunc=None,verbose=True,updateMeta=False,allowNumericMismatch=False,**kwargs,):""" Combines several similar raster files into one single raster file. Parameters ---------- datasets : string or list glob string path describing datasets to combine, alternatively list of gdal.Datasets or iterable object with paths. output : string, optional Filepath to output raster file. If it is an existing file, datasets will be added to output. Recommended to create a new file every time though. If None, no output dataset will be loaded or created on disk and output dataset kept in memory only, by default None combiningFunc : [type], optional Allows customized functions to combine matrices, by default None verbose : bool, optional If True, additional status print stamenets will be issued, by default True updateMeta : bool, optional If True, metadata of output dataset will be a combination of all input rasters, by default False. NOTE: In the case of multiple values for the metadata keys, the last dataset metadata will take precedence. allowNumericMismatch : bool, optional If True, minor deviations in raster context will be ignored/corrected. By default False, i.e. only exactly similar rasters will be combined. **kwargs Will be passed on to geokit.raster.createRaster(). Returns ------- output dataset: osgeo.gdal.Dataset Raster file containing the combined matrices of all input datasets. """ifnotisinstance(allowNumericMismatch,bool):raiseTypeError(f"allowNumericMismatch must be boolean.")# CHECK AND PREPROCESS INPUT DATASETS# Ensure we have a list of raster datasetsifisinstance(datasets,str):datasets=glob(datasets)datasets.sort()elifisinstance(datasets,gdal.Dataset):datasets=[datasets,]else:# assume datasets is iterabledatasets=list(datasets)# make sure we do actually have rasters and they are all indeed "similar"iflen(datasets)==0:raiseGeoKitError("No datasets given")rtol=0ifallowNumericMismatch:rtol=0.0001# allow a numeric deviation of 0.01%datasets=checkSimilarRasters(datasets=datasets,rtol=rtol)# GET REFERENCE CONTEXT FOR THE OUTPUT RASTER# determine info for all datasetsraster_info_list=[rasterInfo(sourceDS=d,compute_statistics=True)fordindatasets]# get reference srs - are all the same thanks to checkSimilarRasterssrs_ref=raster_info_list[0].srs# get the unique actual dtypes in input rastersdata_type_list=[]minimum_and_maximum_values=[]forcurrent_raster_infoinraster_info_list:ifcurrent_raster_info.dtypeisNone:raiseGeoKitError("Input raster has no dtype.")data_type_list.append(current_raster_info.data_type_name_str)minimum_and_maximum_values.append(current_raster_info.minimum_value)minimum_and_maximum_values.append(current_raster_info.maximum_value)# now get the most lightweight commonly usable dtypegdal_data_type_as_string=MinimumCDataTypeHandler.get_valid_gdal_data_type_as_string(list_of_numbers=minimum_and_maximum_values,minimum_gdal_type_list=data_type_list)# get_common_dtype(dtypes=dtypes, fallback=None)# get the reference resolution in x and y dir as the most commonly used valuedx_ref=statistics.mode([_i.pixelWidthfor_iinraster_info_list])dy_ref=statistics.mode([_i.pixelHeightfor_iinraster_info_list])# try to align all bounds to the first matching raster which has correct x_ref and y_ref resolutioni_match=next((ifori,(x,y)inenumerate(zip([_i.pixelWidthfor_iinraster_info_list],[_i.pixelHeightfor_iinraster_info_list],))ifx==dx_refandy==dy_ref),None,)ifi_matchisnotNone:# we have a "perfect" raster, use the min. bounds of that one as reference for all other rastersboundsXmin_ref=raster_info_list[i_match].bounds[0]boundsYmin_ref=raster_info_list[i_match].bounds[2]else:# we do not have any raster which matches both r_ref any y_ref in its resolution. Simply use the first raster.boundsXmin_ref=raster_info_list[0].bounds[0]boundsYmin_ref=raster_info_list[0].bounds[2]ifverbose:print(datetime.datetime.now(),f"NOTE: None of the rasters matches both reference resolutions in x and y direction. Use first raster bounds as reference.",flush=True,)# calculate the possibly adapted bounds for all datasetsboundsSet=[]for_infoinraster_info_list:# calculate the new bounds by aligning bottom left corner with boundsXmin_ref/boundsYmin_ref + multiple of cell size_bounds_Xmin=boundsXmin_ref+round((_info.bounds[0]-boundsXmin_ref)/dx_ref)*dx_ref_bounds_Ymin=boundsYmin_ref+round((_info.bounds[1]-boundsYmin_ref)/dy_ref)*dy_ref_bounds=(_bounds_Xmin,_bounds_Ymin,_bounds_Xmin+_info.xWinSize*dx_ref,_bounds_Ymin+_info.yWinSize*dy_ref,)# make sure the number of cells would remain the same in the new boundsassertround(_info.xWinSize*dx_ref/_info.dx,0)==_info.xWinSize,(f"The change in bounds width would lead to a different amount of cell columns in the original raster.")assertround(_info.yWinSize*dy_ref/_info.dy,0)==_info.yWinSize,(f"The change in bounds height would lead to a different amount of cell rows in the original raster.")boundsSet.append(_bounds)# get summary info about the whole dataset groupdataXMin=min([i[0]foriinboundsSet])dataXMax=max([i[2]foriinboundsSet])dataYMin=min([i[1]foriinboundsSet])dataYMax=max([i[3]foriinboundsSet])# get noData value from kwargs, else take from dataset infosnoData_ref=kwargs.pop("noData",None)ifnoData_refisNone:noDataSet=set([i.noDataforiinraster_info_list])assertlen(noDataSet)==1# make sure, is enforced by checkSimilarRastersnoData_ref=noDataSet.pop()# Maybe create a new output datasetifisinstance(output,str):ifnotos.path.isfile(output):# we will need to create a output sourcecreateRaster(bounds=(dataXMin,dataYMin,dataXMax,dataYMax),output=output,dtype=gdal_data_type_as_string,pixelWidth=dx_ref,pixelHeight=dy_ref,noData=noData_ref,srs=srs_ref,fill=noData_ref,**kwargs,)else:warn("WARNING: Overwriting existing output file. Sometimes writing to an non empty output fails. Recommended to write to a non existing location instead and include maser into datasets.")elifoutputisNone:# create raster in memoryoutputDS=createRaster(bounds=(dataXMin,dataYMin,dataXMax,dataYMax),dtype=gdal_data_type_as_string,pixelWidth=dx_ref,pixelHeight=dy_ref,noData=noData_ref,srs=srs_ref,fill=noData_ref,**kwargs,)else:raiseTypeError("output must be None or a str formatted file path to an existing output file or a file to be created.")# open output dataset if required and check parametersifnotoutputisNone:outputDS=gdal.Open(output,gdal.GA_Update)mInfo=rasterInfo(outputDS)mExtent=Extent(mInfo.bounds,srs=mInfo.srs)outputBand=outputDS.GetRasterBand(1)# Make a meta containerifupdateMeta:meta=outputDS.GetMetadata_Dict()# Add each dataset to outputforiinrange(len(datasets)):ifverbose:print(datetime.datetime.now(),f"Now adding raster No. {i+1}/{len(datasets)}",)# create dataset extentdExtent=Extent(boundsSet[i],srs=srs_ref)# extract the dataset's matrixdMatrix=extractMatrix(datasets[i])ifnotraster_info_list[i].yAtTop:dMatrix=dMatrix[::-1,:]# Calculate starting indicesidx=mExtent.findWithin(dExtent,(mInfo.dx,mInfo.dy),yAtTop=mInfo.yAtTop)# Get output datamMatrix=outputBand.ReadAsArray(xoff=idx.xStart,yoff=idx.yStart,win_xsize=idx.xWin,win_ysize=idx.yWin)ifmMatrixisNone:raiseGeoKitError("mMatrix is None")# create selectorifnotcombiningFuncisNone:# update rasterInfo since we might have slightly changed cell/bounds inforInfo_dict=raster_info_list[i]._asdict()rInfo_dict["bounds"]=boundsSet[i]rInfo_dict["pixelWidth"]=dx_refrInfo_dict["pixelHeight"]=dy_refrInfo_dict["dx"]=dx_refrInfo_dict["dy"]=dy_refrInfo_dict["xMin"]=boundsSet[i][0]rInfo_dict["yMin"]=boundsSet[i][1]rInfo_dict["xMax"]=boundsSet[i][2]rInfo_dict["yMax"]=boundsSet[i][3]rInfo_upd=RasterInfo(**rInfo_dict)writeMatrix=combiningFunc(mMatrix=mMatrix,mInfo=mInfo,dMatrix=dMatrix,dInfo=rInfo_upd)elifnotraster_info_list[i].noDataisNone:sel=dMatrix!=raster_info_list[i].noDatamMatrix[sel]=dMatrix[sel]writeMatrix=mMatrixelse:writeMatrix=dMatrix# Add to outputoutputBand.WriteArray(writeMatrix,idx.xStart,idx.yStart)outputBand.FlushCache()# update metaData, maybeifupdateMeta:meta.update(raster_info_list[i].meta)ifupdateMeta:outputDS.SetMetadata(meta)# Write final rasteroutputDS.FlushCache()outputBand.ComputeRasterMinMax(0)outputBand.ComputeBandStats(0)returnoutputDS
Compare two lists of geometries and return a list of booleans indicating whether each pair of geometries are equal.
The order of the lists is important, as the first geometry in the first list will be compared to the first geometry in the second list, and so on.
Returns:
list
–
A list of booleans indicating whether each pair of geometries are equal.
defcompare_geoms(geoms_1,geoms_2):"""Compare two lists of geometries and return a list of booleans indicating whether each pair of geometries are equal. The order of the lists is important, as the first geometry in the first list will be compared to the first geometry in the second list, and so on. Returns ------- list A list of booleans indicating whether each pair of geometries are equal. """equal=map(lambdag1,g2:g1.Equals(g2),geoms_1,geoms_2)returnlist(equal)
Create contour geometries at specified edges for the given raster data.
Notes
This function is similar to geokit.geom.polygonizeMatrix, although it only
operates on the user-specified edges AND applies the 'Marching Squares'
algorithm
See the gdal function "GDALContourGenerateEx" for more information on the
specifics of this algorithm
The edges to search for within the raster dataset
* This parameter can be set as "None", in which case an additional
argument should be given to specify how the edges should be determined
- See the documentation of "GDALContourGenerateEx"
- Ex. "LEVEL_INTERVAL=10", contourEdges=None
defcontours(source:load_raster_input,contourEdges:list[float]|None,polygonize:bool=True,unpack:bool=True,raster_band_index:int=1,**kwargs,)->pd.DataFrame:"""Create contour geometries at specified edges for the given raster data. Notes ----- This function is similar to geokit.geom.polygonizeMatrix, although it only operates on the user-specified edges AND applies the 'Marching Squares' algorithm See the gdal function "GDALContourGenerateEx" for more information on the specifics of this algorithm Parameters ---------- source : Anything acceptable by loadRaster() The raster datasource to operate on contourEdges : list[float] | None The edges to search for within the raster dataset * This parameter can be set as "None", in which case an additional argument should be given to specify how the edges should be determined - See the documentation of "GDALContourGenerateEx" - Ex. "LEVEL_INTERVAL=10", contourEdges=None polygons : bool If true, contours are returned as polygons instead of linstrings unpack : bool If True, Multipolygon/MultiLinestring objects are decomposed raster_band_index : int The index of the raster which should be loaded from the raster band. **kwargs: * All keyword arguments are passed on to a call to gdal.ContourGenerateEx * They are used to construct the 'options' parameter * Example keys include: LEVEL_INTERVAL, LEVEL_BASE, LEVEL_EXP_BASE, NODATA * Do not use the key "ID_FIELD", since this is employed already Returns ------- pandas.DataFrame * The column 'geom' corresponds to generated geometry objects * The columns 'ID' corresponds to the associated contour edge for each object """warnings.warn(message="The current behavior of geokits's contours function is deprecated. GDAL has changed"" how contours are drawn close to minimum and maximum values, and will discontinue"" the current behavior in GDAL 3.11 or 3.12. Geokit will also drop the current"" behavior with the next GDAL update. For more information, please see the"" following discussion: https://github.com/OSGeo/gdal/issues/12938.",category=DeprecationWarning,)# Open rasterraster=loadRaster(source)band=raster.GetRasterBand(raster_band_index)rasterSRS=SRS.loadSRS(raster.GetProjectionRef())# Make temporary vectordriver:Driver=gdal.GetDriverByName("Memory")source:gdal.Dataset=driver.Create("",0,0,0,gdal.GDT_Unknown)layer:ogr.Layer=source.CreateLayer("",rasterSRS,ogr.wkbPolygonifpolygonizeelseogr.wkbLineString)field=ogr.FieldDefn("DN",ogr.OFTInteger)layer.CreateField(field)# Setup contour functionargs=["{}={}".format(a,b)fora,binkwargs.items()]args.append("ID_FIELD=DN")ifpolygonize:args.append("POLYGONIZE=YES")ifcontourEdgesisnotNone:opt="FIXED_LEVELS="foredgeincontourEdges:opt+=str(edge)+","args.append(opt[:-1])result=gdal.ContourGenerateEx(band,layer,options=args)ifnotresult==gdal.CE_None:raiseGeoKitRasterError("Failed to compute raster contours")layer.CommitTransaction()IDs=[]geoms=[]iterator_n=0forftridinrange(layer.GetFeatureCount()):ftr:ogr.Feature=layer.GetFeature(ftrid)geom=ftr.GetGeometryRef()value=ftr.GetField(0)print(iterator_n)iterator_n=iterator_n+1ifunpack:forgiinrange(geom.GetGeometryCount()):geoms.append(geom.GetGeometryRef(gi).Clone())IDs.append(value)else:geoms.append(geom.Clone())IDs.append(value)countour_data_frame=pd.DataFrame(dict(geom=geoms,ID=IDs))# return geomsreturncountour_data_frame
defconvertGeoJson(geojson,srs=3857):"""Make a geometry from a well known text (WKT) string. TODO: UPDATE!!! Parameters ---------- wkt : str The WKT string to convert srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the geometry to create """geom=ogr.CreateGeometryFromJson(geojson)# Create new geometry from stringifgeomisNone:# test for successraiseGeoKitGeomError("Failed to create geometry")ifsrs:geom.AssignSpatialReference(SRS.loadSRS(srs))# Assign the given srsreturngeom
defconvertWKT(wkt,srs=None):"""Make a geometry from a well known text (WKT) string. Parameters ---------- wkt : str The WKT string to convert srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the geometry to create """geom=ogr.CreateGeometryFromWkt(wkt)# Create new geometry from stringifgeomisNone:# test for successraiseGeoKitGeomError("Failed to create geometry")ifsrs:geom.AssignSpatialReference(SRS.loadSRS(srs))# Assign the given srsreturngeom
An SQL-like where statement to apply to the source
* Feature attribute name do not need quotes
* String values should be wrapped in 'single quotes'
Example: If the source vector has a string attribute called "ISO" and
a integer attribute called "POP", you could use....
defcountFeatures(source,geom=None,where=None):"""Returns the number of features found in the given source and within a given geometry and/or where-statement. Parameters ---------- source : Anything acceptable by loadVector() The vector datasource to read from geom : ogr.Geometry; optional The geometry to search within * All features are extracted which touch this geometry where : str; optional An SQL-like where statement to apply to the source * Feature attribute name do not need quotes * String values should be wrapped in 'single quotes' Example: If the source vector has a string attribute called "ISO" and a integer attribute called "POP", you could use.... where = "ISO='DEU' AND POP>1000" Returns ------- int """ds=loadVector(source)layer=ds.GetLayer()filterLayer(layer,geom,where)returnlayer.GetFeatureCount()
defcreateDataFrameFromGeoDataFrame(gdf:gpd.GeoDataFrame)->pd.DataFrame:"""Creates a geokit-style dataframe from a geopandas geodataframe. Parameters ---------- gdf : pd.DataFrame geopandas-style pd.DataFrame, need a 'geometry' column. Returns ------- pd.DataFrame Same as the previous, just as an geokit-style dataframe with 'geom' column with osgeo.ogr.Geometry objects. """assertisinstance(gdf,pd.DataFrame)assert"geometry"ingdf.columns# # import the required external packages - these are not part of the requirements.yml and are possibly not installed# try:# import geopandas as gpd# except:# raise ImportError(# "'geopandas' is required for geokit.vector.createGeoDataFrame() but is not installed in the current environment."# )# get the CRS/SRS of the gdfcrs_wkt=gdf.crs.to_wkt()srs=SRS.loadSRS(crs_wkt)# create WKT geometry and then convert to osgeo.ogr.Geometrygdf["geom"]=gdf["geometry"].apply(lambdax:GEOM.convertWKT(x.wkt,srs=srs))# remove geometry column and return geokit dfdf=gdf.drop(columns="geometry")returndf
defcreateGeoDataFrame(dfGeokit:pd.DataFrame)->gpd.GeoDataFrame:"""Creates a gdf from an Reskit shape pd.DataFrame. Parameters ---------- dfGeokit : pd.DataFrame Reskit shape pd.DataFrame, need a 'geom' column. Returns ------- gpd.GeoDataFrame Same as the previous, just as an GeodataFrame """assertisinstance(dfGeokit,pd.DataFrame)assert"geom"indfGeokit.columns# import the required external packages - these are not part of the requirements.yml and are possibly not installedtry:importshapelyexcept:raiseImportError("'shapely' is required for geokit.vector.createGeoDataFrame() but is not installed in the current environment.")try:importgeopandasasgpdexcept:raiseImportError("'geopandas' is required for geokit.vector.createGeoDataFrame() but is not installed in the current environment.")# get valuesvalues={}forcolindfGeokit.columns:# geoms need to be converted to shapelyifcol=="geom":values["geometry"]=[shapely.wkt.loads(g.ExportToWkt())forgindfGeokit.geom]# this takes some time :/# other values are just stored as they areelse:values[col]=list(dfGeokit[col])# get srs as Well known textcrs=dfGeokit.geom.iloc[0].GetSpatialReference().ExportToWkt()# create gdf and set indexgdf=gpd.GeoDataFrame(values,crs=crs)gdf=gdf.set_index(dfGeokit.index,drop=True)# over and out!returngdf
defcreateGeoJson(geoms,output=None,srs=4326,topo=False,fill=""):"""Convert a set of geometries to a geoJSON object."""ifsrs:srs=SRS.loadSRS(srs)# arrange geom, index, and dataifisinstance(geoms,ogr.Geometry):# geoms is a single geometryfinalGeoms=[geoms,]data=Noneindex=[0,]elifisinstance(geoms,pd.Series):index=geoms.indexfinalGeoms=geoms.valuesdata=Noneelifisinstance(geoms,pd.DataFrame):index=geoms.indexfinalGeoms=geoms.geom.valuesdata=geoms.loc[:,geoms.columns!="geom"]data["_index"]=indexelse:finalGeoms=list(geoms)data=Noneindex=list(range(len(finalGeoms)))iflen(finalGeoms)==0:raiseGeoKitVectorError("Empty geometry list given")# Transform?ifnotsrsisNone:finalGeoms=GEOM.transform(finalGeoms,toSRS=srs)# Make JSON objectfromioimportBytesIOifnotoutputisNoneandnotisinstance(output,str):ifnotoutput.writable():raiseGeoKitVectorError("Output object is not writable")iftopo:fo=BytesIO()else:fo=outputelifisinstance(output,str)andnottopo:fo=open(output,"wb")else:fo=BytesIO()fo.write(bytes('{"type":"FeatureCollection","features":[',encoding="utf-8"))forj,i,ginzip(range(len(index)),index,finalGeoms):fo.write(bytes('%s{"type":"Feature",'%(""ifj==0else","),encoding="utf-8"))ifdataisNone:fo.write(bytes('"properties":{"_index":%s},'%str(i),encoding="utf-8"))else:fo.write(bytes('"properties":%s,'%data.loc[i].fillna(fill).to_json(),encoding="utf-8",))fo.write(bytes('"geometry":%s}'%g.ExportToJson(),encoding="utf-8"))# fo.write(bytes('"geometry": {"type": "Point","coordinates": [125.6, 10.1] }}', encoding='utf-8'))fo.write(bytes("]}",encoding="utf-8"))fo.flush()# Put in the right formatiftopo:fromioimportTextIOWrapperfromtopojsonimportconversionfo.seek(0)topo=conversion.convert(TextIOWrapper(fo),object_name="primary")# automatically closes fotopo=str(topo).replace("'",'"')# Done!ifoutputisNone:iftopo:returntopoelse:fo.seek(0)geojson=fo.read()fo.close()returngeojson.decode("utf-8")elifisinstance(output,str):iftopo:withopen(output,"w")asfo:fo.write(topo)else:pass# we already wrote to the file!returnoutputelse:iftopo:output.write(bytes(topo,encoding="utf-8"))else:pass# We already wrote to the file!returnNone
Raster datasets are always written in the 'yAtTop' orientation. Meaning that
the first row of data values (either written to or read from the dataset) will
refer to the TOP of the defined boundary, and will then move downward from
there
If a data matrix is given, and a negative pixelWidth is defined, the data
will be flipped automatically
The pixel width of the raster in units of the input srs
* The keyword 'dx' can be used as well and will override anything given
assigned to 'pixelWidth'
The pixel height of the raster in units of the input srs
* The keyword 'dy' can be used as well and will override anything given
assigned to 'pixelHeight'
A path to an output file
* If output is None, the raster will be created in memory and a dataset
handle will be returned
* If output is given, the raster will be written to disk and nothing will
be returned
The datatype of the represented by the created raster's band
* Options are: Byte, Int16, Int32, Int64, Float32, Float64
* If dtype is None and data is None, the assumed datatype is a 'Byte'
* If dtype is None and data is not None, the datatype will be inferred
from the given data
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs of the point to create
* If not given, no srs will be assigned to the created raster. Please be
aware that some operations may not work correctly if no SRS is given.
* If 'bounds' is an Extent object, the bounds' internal srs will override
this input
A flag instructing the output raster to use a compression algorithm
* only useful if 'output' has been defined
* "DEFLATE" used for Linux/Mac, "LZW" used for Windows
Specifies which value should be considered as 'no data' in the created
raster
* Must be the same datatype as the 'dtype' input (or that which is derived)
defcreateRaster(bounds:tuple[numeric,numeric,numeric,numeric],output:None|str|pathlib.Path=None,pixelWidth:numeric=100,pixelHeight:numeric=100,dtype:None|geokit_c_data_types_literal=None,srs:srs_input|None=None,compress:bool=True,noData:numeric|None=None,overwrite:bool=True,fill:numeric|None=None,data:np.ndarray|None=None,meta:dict|None=None,scale:numeric|None=1,offset:numeric=0,creationOptions:dict=dict(),raster_band_index:int=1,)->gdal.Dataset|str:"""Create a raster file. NOTE: ----- Raster datasets are always written in the 'yAtTop' orientation. Meaning that the first row of data values (either written to or read from the dataset) will refer to the TOP of the defined boundary, and will then move downward from there If a data matrix is given, and a negative pixelWidth is defined, the data will be flipped automatically Parameters ---------- bounds : (xMin, yMin, xMax, yMax) or Extent The geographic extents spanned by the raster pixelWidth : numeric The pixel width of the raster in units of the input srs * The keyword 'dx' can be used as well and will override anything given assigned to 'pixelWidth' pixelHeight : numeric The pixel height of the raster in units of the input srs * The keyword 'dy' can be used as well and will override anything given assigned to 'pixelHeight' output : str, pathlib.Path, None A path to an output file * If output is None, the raster will be created in memory and a dataset handle will be returned * If output is given, the raster will be written to disk and nothing will be returned dtype : geokit_c_data_types_literal, none The datatype of the represented by the created raster's band * Options are: Byte, Int16, Int32, Int64, Float32, Float64 * If dtype is None and data is None, the assumed datatype is a 'Byte' * If dtype is None and data is not None, the datatype will be inferred from the given data srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the point to create * If not given, no srs will be assigned to the created raster. Please be aware that some operations may not work correctly if no SRS is given. * If 'bounds' is an Extent object, the bounds' internal srs will override this input compress : bool A flag instructing the output raster to use a compression algorithm * only useful if 'output' has been defined * "DEFLATE" used for Linux/Mac, "LZW" used for Windows noData : numeric; optional Specifies which value should be considered as 'no data' in the created raster * Must be the same datatype as the 'dtype' input (or that which is derived) fill : numeric; optional The initial value given to all pixels in the created raster band - numeric * Must be the same datatype as the 'dtype' input (or that which is derived) overwrite : bool A flag to overwrite a pre-existing output file * If set to False and an 'output' is specified which already exists, an error will be raised data : matrix_like A 2D matrix to write into the resulting raster * array dimensions must fit raster dimensions as calculated by the bounds and the pixel resolution meta : dictionary, None of key values pairs that is stored in the raster as meta data. Is written to the raster using the scale : numeric; optional The scaling value given to apply to all values - numeric * Must be the same datatype as the 'dtype' input (or that which is derived) offset : numeric; optional The offset value given to apply to all values - numeric * Must be the same datatype as the 'dtype' input (or that which is derived) raster_band_index: int, defaults to 1 Determines which band is written to in the output raster dataset. Returns ------- * If 'output' is None: gdal.Dataset * If 'output' is a string: The path to the output is returned (for easy opening). It has to be saved as geotiff with .tif suffix. """# Check for existing fileifoutputisnotNone:output=str(output)output_path=pathlib.Path(output)ifoutput_path.exists():ifoverwriteisTrue:output_path.unlink()output_path_with_aux_extension=output_path.with_suffix(".aux.xml")ifoutput_path_with_aux_extension.exists():output_path_with_aux_extension.unlink()else:raiseGeoKitRasterError("Output file already exists: %s"%output)# check if the directory existselifnotoutput_path.parent.is_dir():raiseFileNotFoundError(f"Output directory does not exist: {os.path.dirname(output)}")# check if writeable:ifnotos.access(os.path.dirname(output),os.W_OK):print(os.path.dirname(output))raisePermissionError(f"Writing permission error for path: {os.path.dirname(output)}")# Ensure bounds is okay# bounds = UTIL.fitBoundsTo(bounds, pixelWidth, pixelHeight)# Make a raster dataset and pull the band/maskBand objectsx_min=bounds[0]y_min=bounds[1]x_max=bounds[2]y_max=bounds[3]# Always use the "Y-at-Top" orientationcols=int(round((x_max-x_min)/pixelWidth))rows=int(round((y_max-y_min)/abs(pixelHeight)))list_of_numbers=[]minimum_gdal_type_list=[]ifisinstance(dtype,str):minimum_gdal_type_list.append(dtype)ifisinstance(noData,numeric):list_of_numbers.append(noData)ifisinstance(fill,numeric):list_of_numbers.append(fill)ifisinstance(data,np.ndarray):numpy_data_type=str(data.dtype)minimum_gdal_type_list.append(numpy_data_type)list_of_numbers.append(data.min())list_of_numbers.append(data.max())elifdataisNone:passelse:raiseGeoKitRasterError("Data must be given as a numpy ndarray or None")data_type_constant=MinimumCDataTypeHandler.get_valid_gdal_data_type_as_constant(list_of_numbers=list_of_numbers,minimum_gdal_type_list=minimum_gdal_type_list)# # Get DataType# if dtype is not None: # a dtype was given, use it!# dtype = gdalType(dtype)# elif data is not None: # a data matrix was give, use it's dtype! (assume a numpy array or derivative)# dtype = gdalType(data.dtype)# else: # Otherwise, just assume we want a Byte# dtype = "GDT_Byte"# Open the driveropts=OrderedDict()ifcompressandoutputisnotNone:opts["COMPRESS"]=COMPRESSION_OPTION_STRifcreationOptionsisnotNone:opts.update(creationOptions)opts=["{}={}".format(k,v)fork,vinopts.items()]ifoutputisNone:driver:gdal.Driver=gdal.GetDriverByName("Mem")# create a raster in memoryraster:gdal.Dataset=driver.Create("",cols,rows,1,data_type_constant,opts)else:driver:gdal.Driver=gdal.GetDriverByName("GTiff")# Create a raster in storageraster:gdal.Dataset=driver.Create(output,cols,rows,1,data_type_constant,opts)ifrasterisNone:raiseGeoKitRasterError("Failed to create raster")# Do the rest in a "try" statement so that a failure won't bind the sourcetry:raster.SetGeoTransform((x_min,abs(pixelWidth),0,y_max,0,-1*abs(pixelHeight)))# Set the SRSifsrsisnotNone:rasterSRS=SRS.loadSRS(srs)raster.SetProjection(rasterSRS.ExportToWkt())else:warnings.warn(message="No srs given when creating raster. Please be aware that some operations may not work correctly.",category=UserWarning,)# Fill the raster will zeros, null values, or initial values (if given)band:gdal.Band=raster.GetRasterBand(raster_band_index)ifscaleisnotNone:band.SetScale(scale)ifoffsetisnotNone:band.SetOffset(offset)ifnoDataisnotNone:band.SetNoDataValue(noData)iffillisNoneanddataisNone:band.Fill(noData)ifdataisNone:iffillisNone:band.Fill(0)else:band.Fill(fill)else:# make sure dimension size is goodifnot(data.shape[0]==rowsanddata.shape[1]==cols):raiseGeoKitRasterError("Raster dimensions and input data dimensions do not match.The data has rows="+str(data.shape[0])+" and columns="+str(data.shape[1])+". The raster specification expects rows="+str(rows)+" and columns="+str(cols))# See if data needs flippingifpixelHeight<0:data=data[::-1,:]# Write it!band.WriteArray(data)band.FlushCache()band.ComputeRasterMinMax(0)band.ComputeBandStats(0)# Write MetaData, maybeifmetaisnotNone:fork,vinmeta.items():raster.SetMetadataItem(k,v)# writes the raster to the hard driveraster.FlushCache()# Return raster if in memoryifoutputisNone:returnraster# Donereturnoutput# Handle the fail caseexceptExceptionase:raster=Noneraisee
defcreateRasterLike(source:load_raster_input|RasterInfo,copyMetadata:bool=True,metadata:dict|None=None,data:None|np.ndarray=None,dtype:None|geokit_c_data_types_literal=None,**kwargs,):"""Create a raster described by the given raster info (as returned from a call to rasterInfo() ). * This copies all characteristics of the given raster, including: bounds, pixelWidth, pixelHeight, dtype, srs, noData, and meta. * Any keyword argument which is given will override values found in the source """ifUTIL.isRaster(source):raster_info=rasterInfo(source)elifisinstance(source,RasterInfo):raster_info=sourceelse:raiseGeoKitRasterError("Could not understand source")ifcopyMetadataandmetadataisnotNone:raiseGeoKitRasterError("If metadata is given, copyMetadata cannot be True!")bounds=kwargs.pop("bounds",raster_info.bounds)pixelWidth=kwargs.pop("pixelWidth",raster_info.pixelWidth)pixelHeight=kwargs.pop("pixelHeight",raster_info.pixelHeight)srs=kwargs.pop("srs",raster_info.srs)noData=kwargs.pop("noData",raster_info.noData)data_type_as_string=kwargs.pop("data_type_as_string",dtype)ifcopyMetadata:meta=kwargs.pop("meta",raster_info.meta)else:meta=metadatareturncreateRaster(bounds=bounds,pixelWidth=pixelWidth,pixelHeight=pixelHeight,srs=srs,noData=noData,meta=meta,dtype=data_type_as_string,data=data,**kwargs,)
The geometries to write into the vector file
* If a DataFrame is given, it must have a column called 'geom'
* All geometries must share the same type (point, line, polygon, etc...)
* All geometries must share the same SRS
* If geometry SRS differs from the 'srs' input, then all geometries will
be projected to the input srs
A path on disk to create the output vector
* If output is None, the vector dataset will be created in memory
* Assumed to be of "ESRI Shapefile" format
* Will create a number of files with different extensions
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs of the vector to create
* If not given, the geometries' inherent srs is used
* If srs does not match the inherent srs, all geometries will be
transformed
(dict of lists or pandas.DataFrame; optional, default:
None
)
–
Explicit attribute values to assign to each geometry
* If a dict is given, must be a dictionary of the attribute names and
associated lists of the attribute values
* If a DataFrame is given, the field names are taken from column names
and attribute values from the corresponding column data
* The order of each column/list will correspond to the order geometries
are written into the dataset
* The length of each column/list MUST match the number of geometries
* All values in a single column/list must share the same type
- Options are int, float, or str
A dictionary specifying the datatype of each attribute when written into
the final dataset
* Options are defined from ogr.OFT[...]
- ex. Integer, Real, String
* The ogrType() function can be used to map typical python and numpy types
to appropriate ogr types
If True, all geoms will be asserted in object type and exact srs. Else, only
the first geom in the geom column/iterable will be checked for performance reasons.
By default False.
defcreateVector(# geoms: ogr.Geometry | list[ogr.Geometry | gdal.Dataset | str] | pd.DataFrame | np.ndarray | str | gdal.Dataset,geoms,output:str|None=None,srs=None,driverName:str="ESRI Shapefile",layerName:str="default",fieldVals=None,fieldDef=None,checkAllGeoms=False,overwrite:bool=True,):""" Create a vector on disk from geometries or a DataFrame with 'geom' column. Parameters ---------- geoms : ogr.Geometry or [ogr.Geometry, ] or pandas.DataFrane The geometries to write into the vector file * If a DataFrame is given, it must have a column called 'geom' * All geometries must share the same type (point, line, polygon, etc...) * All geometries must share the same SRS * If geometry SRS differs from the 'srs' input, then all geometries will be projected to the input srs output : str; optional A path on disk to create the output vector * If output is None, the vector dataset will be created in memory * Assumed to be of "ESRI Shapefile" format * Will create a number of files with different extensions srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the vector to create * If not given, the geometries' inherent srs is used * If srs does not match the inherent srs, all geometries will be transformed driverName : str; optional The name of the driver to use when creating the vector. Currently supported options are: - ESRI Shapefile - GPKG For a list of all supported vector drivers by OGR, see: https://gdal.org/drivers/vector/index.html layerName : str; optional The name of the layer to create within the vector. (Only applicable for GeoPackages) * If the layer already exists, it will be overwritten fieldVals : dict of lists or pandas.DataFrame; optional Explicit attribute values to assign to each geometry * If a dict is given, must be a dictionary of the attribute names and associated lists of the attribute values * If a DataFrame is given, the field names are taken from column names and attribute values from the corresponding column data * The order of each column/list will correspond to the order geometries are written into the dataset * The length of each column/list MUST match the number of geometries * All values in a single column/list must share the same type - Options are int, float, or str fieldDef : dict; optional A dictionary specifying the datatype of each attribute when written into the final dataset * Options are defined from ogr.OFT[...] - ex. Integer, Real, String * The ogrType() function can be used to map typical python and numpy types to appropriate ogr types checkAllGeoms : bool, optional If True, all geoms will be asserted in object type and exact srs. Else, only the first geom in the geom column/iterable will be checked for performance reasons. By default False. overwrite : bool; optional Determines whether the preexisting files should be overwritten * Only used when output is not None Returns ------- * If 'output' is None: gdal.Dataset * If 'output' is given: None """ifsrs:srs=SRS.loadSRS(srs)# make geom or wkt list into a list of ogr.Geometry objectsfinalGeoms=[]geomIndex=Noneifisinstance(geoms,str)orisinstance(geoms,ogr.Geometry):# geoms is a single geometrygeoms=[geoms,]elifisinstance(geoms,pd.Series):geomIndex=geoms.indexgeoms=geoms.valueselifisinstance(geoms,pd.DataFrame):ifnotfieldValsisNone:raiseGeoKitVectorError("fieldVals must be None when geoms input is a DataFrame")fieldVals=geoms.copy()geoms=geoms.geom.valuesfieldVals.drop("geom",inplace=True,axis=1)iflen(geoms)==0:raiseGeoKitVectorError("Empty geometry list given")# Test if the first geometry is an ogr-Geometry typeifisinstance(geoms[0],ogr.Geometry):# (Assume all geometries in the array have the same type)geomSRS=geoms[0].GetSpatialReference()ifcheckAllGeoms:# check if all other geoms have the same SRSassertall([geomSRS.IsSame(g.GetSpatialReference())forgingeoms]),(f"Not all geoms have the same SRS, srs of first geom: {geomSRS}")# Set up some variables for the upcoming loopdoTransform=FalsesetSRS=FalseifsrsandgeomSRSandnotsrs.IsSame(geomSRS):# Test if a transformation is neededtrx=osr.CoordinateTransformation(geomSRS,srs)doTransform=True# Test if an srs was NOT given, but the incoming geometries have an SRS alreadyelifsrsisNoneandgeomSRS:srs=geomSRS# Test if an srs WAS given, but the incoming geometries do NOT have an srs alreadyelifsrsandgeomSRSisNone:# In which case, assume the geometries are meant to have the given srssetSRS=True# Create geomsforiinrange(len(geoms)):# clone just in case the geometry is tethered outside of functionfinalGeoms.append(geoms[i].Clone())ifdoTransform:finalGeoms[i].Transform(trx)# Do transform if necessaryifsetSRS:finalGeoms[i].AssignSpatialReference(srs)# Set srs if necessary# Otherwise, the geometry array should contain only WKT stringselifisinstance(geoms[0],str):ifsrsisNone:raiseValueError("srs must be given when passing wkt strings")# Create geomsfinalGeoms=[GEOM.convertWKT(wkt,srs)forwktingeoms]else:raiseValueError("Geometry inputs must be ogr.Geometry objects or WKT strings")geoms=finalGeoms# Overwrite geoms array with the conditioned finalGeoms# Determine geometry typestypes=set()forgingeoms:# Get a set of all geometry type-names (POINT, POLYGON, etc...)types.add(g.GetGeometryName())iftypes.issubset({"POINT","MULTIPOINT"}):geomType=ogr.wkbPointeliftypes.issubset({"LINESTRING","MULTILINESTRING"}):geomType=ogr.wkbLineStringeliftypes.issubset({"POLYGON","MULTIPOLYGON"}):geomType=ogr.wkbPolygonelse:# geomType = ogr.wkbGeometryCollectionraiseRuntimeError("Could not determine output shape's geometry type")# Create a driver and datasource# driver = gdal.GetDriverByName("ESRI Shapefile")# dataSource = driver.Create(output, 0, 0)ifoutputisNone:# Create datasource if no output path is provideddriver=gdal.GetDriverByName("Memory")# Using 'Create' from a Memory driver leads to an error. But creating# a temporary shapefile driver (it does not actually produce a file, I think)# and then using 'CreateCopy' seems to worktmp_driver=gdal.GetDriverByName("ESRI Shapefile")t=TemporaryDirectory()tmp_dataSource=tmp_driver.Create(t.name+"tmp.shp",0,0)dataSource=driver.CreateCopy("MEMORY",tmp_dataSource)t.cleanup()deltmp_dataSource,tmp_driver,telifoutputisnotNone:# Create datasource if output path is providedifoverwriteisTrue:# Search for directoryifos.path.dirname(output)=="":# If no directory is given, assume current directoryoutput=os.path.join(os.getcwd(),output)elifnotos.path.isdir(os.path.dirname(output)):# If directory does not exist, raise errorraiseFileNotFoundError(f"Directory {os.path.dirname(output)} does not exist")# Remove file if it existsifos.path.isfile(output):os.remove(output)driver=ogr.GetDriverByName(driverName)dataSource=driver.CreateDataSource(output)elifoverwriteisFalse:dataSource=ogr.Open(output,1)assertdataSourceisnotNone,f"Could not open {output}"# Wrap the whole writing function in a 'try' statement in case it failstry:# Create the layerifoutputisnotNoneandoverwriteisFalse:layerName=layerNameassertlayerNamenotinlistLayers(output),(f"Layer name '{layerName}' already exists in {output}. Please Specify a new layer name or set overwrite=True.")else:layerName=layerNamelayer=dataSource.CreateLayer(layerName,srs,geomType)assertlayerisnotNone,"Could not create layer!"# Setup fieldVals and fieldDef dictsifnotfieldValsisNone:# Ensure fieldVals is a dataframeifnotisinstance(fieldVals,pd.DataFrame):# If not, try converting itfieldVals=pd.DataFrame(fieldVals)ifnotgeomIndexisNone:fieldVals=fieldVals.loc[geomIndex]# check if length is goodiffieldVals.shape[0]!=len(geoms):raiseGeoKitVectorError("Values table length does not equal geometry count")# Ensure fieldDefs exists and is a dictiffieldDefisNone:# Try to determine data typesfieldDef={}fork,vinfieldVals.items():fieldDef[k]=ogrType(v.dtype)elifisinstance(fieldDef,dict):# Assume fieldDef is a dict mapping fields to inappropriate# datatypes. Fix these to ogr indicators!forkinfieldVals.columns:try:fieldDef[k]=ogrType(fieldDef[k])exceptKeyErrorase:print('"%s" not in attribute definition table'%k)raiseeelse:# Assume a single data type was given. Apply to all fields_type=ogrType(fieldDef)fieldDef={}forkinfieldVals.keys():fieldDef[k]=_type# Write field definitions to layerforfieldName,dtypeinfieldDef.items():layer.CreateField(ogr.FieldDefn(str(fieldName),getattr(ogr,dtype)))# Ensure list lengths match geom lengthfork,vinfieldVals.items():iflen(v)!=len(geoms):raiseRuntimeError("'{}' length does not match geom list".format(k))# Create featuresforgiinrange(len(geoms)):# Create a blank featurefeature=ogr.Feature(layer.GetLayerDefn())# Fill the attributes, if requiredifnotfieldValsisNone:forfieldName,valueinfieldVals.items():_type=fieldDef[fieldName]# cast to basic typeif_type=="OFTString":val=str(value.iloc[gi])elif_type=="OFTInteger"or_type=="OFTInteger64":val=int(value.iloc[gi])else:val=float(value.iloc[gi])# Write to featurefeature.SetField(str(fieldName),val)# Set the Geometryfeature.SetGeometry(geoms[gi])# Create the featurelayer.CreateFeature(feature)feature.Destroy()# Free resources (probably not necessary here)# Finishifoutput:returnoutputelse:returndataSource# Delete the datasource in case it failedexceptExceptionase:dataSource=Noneraisee
Multipolygons spanning the antimeridian (this includes polygons that are
split at the antimeridian, with the Western half shifted Eastwards by 360°
longitude) are separated into a part East and West of the antimeridian by
identifying the largest longitudinal gap between any of the sub polygons and
dividing the sub polys into one Eastern and one Western polygon list which
is returned as multipolygons.
NOTE: This function only works for already shifted subpolygons with an
overall envelope between -180° and +180° longitude.
'left' or 'right' to return the left or right side of the antimeridian
'main' to return the side with the largest area
'both' to return both sides as a tuple (left, right)
defdivideMultipolygonIntoEasternAndWesternPart(geom,side="both"):""" Multipolygons spanning the antimeridian (this includes polygons that are split at the antimeridian, with the Western half shifted Eastwards by 360° longitude) are separated into a part East and West of the antimeridian by identifying the largest longitudinal gap between any of the sub polygons and dividing the sub polys into one Eastern and one Western polygon list which is returned as multipolygons. NOTE: This function only works for already shifted subpolygons with an overall envelope between -180° and +180° longitude. Parameters ---------- geom : ogr.Geometry The geometry to split. Must be a MultiPolygon. side : str, optional 'left' or 'right' to return the left or right side of the antimeridian 'main' to return the side with the largest area 'both' to return both sides as a tuple (left, right) """# check inputsassertsidein("both","left","right","main",),"side must be 'left', 'right', 'main' or 'both'"assertisinstance(geom,ogr.Geometry),"geom must be of type osgeo.ogr.Geometry"assertgeom.GetGeometryName()=="MULTIPOLYGON","Only MultiPolygon supported"assertgeom.GetSpatialReference().IsSame(SRS.loadSRS(4326)),"geometry must be in EPSG:4326"assertgeom.GetEnvelope()[0]>=-180andgeom.GetEnvelope()[1]<=180,("Envelope must be between -180° and +180° longitude")assertgeom.GetSpatialReference().IsSame(SRS.loadSRS(4326)),"Only EPSG:4326 lat/lon supported"# first extract sub polygonssub_polys=[geom.GetGeometryRef(i)foriinrange(geom.GetGeometryCount())]# get all the bounding boxesbounds=[]forpolyinsub_polys:env=poly.GetEnvelope()# (minX, maxX, minY, maxY)bounds.append((env[0],env[1],poly))# store (minX, maxX, polygon)# sort them from left to right based on minXbounds.sort(key=lambdax:x[0])# find the largest gap iterativelymax_gap=0split_index=0curr_maxs=list()foriinrange(len(bounds)-1):curr_maxs.append(bounds[i][1])curr_max=max(curr_maxs)next_min=bounds[i+1][0]gap=next_min-curr_maxifgap>max_gap:# overwrite the max gap so far and save the index when it occursmax_gap=gapsplit_index=i# split into two sets of geoms - left and right of the gap (do only if necessary for time)ifsidein["both","main","right"]:# filter only sub polys below (or equal) split indexright_polys=[b[2]fori,binenumerate(bounds)ifi<=split_index]# merge all right polygons into a single multipolygon and assign the same spatial referenceright_multi=ogr.Geometry(ogr.wkbMultiPolygon)forpolyinright_polys:right_multi.AddGeometry(poly.Clone())right_multi.AssignSpatialReference(geom.GetSpatialReference())# same for the left sideifsidein["both","main","left"]:# only polys above split indexleft_polys=[b[2]fori,binenumerate(bounds)ifi>split_index]left_multi=ogr.Geometry(ogr.wkbMultiPolygon)forpolyinleft_polys:left_multi.AddGeometry(poly.Clone())left_multi.AssignSpatialReference(geom.GetSpatialReference())ifside=="left":returnleft_multielifside=="right":returnright_multielifside=="both":# return both left and rightreturnleft_multi,right_multielifside=="main":# return the side with the largest arealeft_area=left_multi.GetArea()right_area=right_multi.GetArea()ifleft_area>right_area:returnleft_multielse:returnright_multielse:raiseValueError("side must be 'left', 'right', 'main' or None")
Each geometry type is displayed as an appropriate plotting type
-> Points/ Multipoints are displayed as points using plt.plot(...)
-> Lines/ MultiLines are displayed as lines using plt.plot(...)
-> Polygons/ MultiPolygons are displayed as patches using the descartes
library
Each geometry can be given its own set of matplotlib plotting parameters
Notes
This function does not call plt.show() for the final display of the figure.
This must be done manually after calling this function. Otherwise
plt.savefig(...) can be called to save the output somewhere.
Sometimes geometries will disappear because of the simplification procedure.
If this happens, the procedure can be avoided by setting simplificationFactor
to None. This will take much more memory and will take longer to plot, however
The geometries to be drawn
* If a DataFrame is given, the function looks for geometries under a
columns named 'geom'
* plotting arguments can be given by adding a column named 'MPL:*'
where '*' stands in for the argument to be added
- For geometries that should ignore this argument, set it as None
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs in which to draw each geometry
* If not given, longitude/latitude is assumed
* Although geometries can be given in any SRS, it is very helpful if
they are already provided in the correct SRS
The level to which geometries should be simplified. It can be thought of
as the number of vertices allowed in either the X or Y dimension across
the figure
* A higher value means a more detailed plot, but may take longer to draw
The figure size to create when generating a new axis
* If resultign figure looks weird, altering the figure size is your best
bet to make it look nicer
The title to give to the generated colorbar
* If not given, but 'colorBy' is given, the same string for 'colorBy'
is used
* Only useful when 'colorBy' is given
All other keyword arguments are passed on to the plotting functions called
for each geometry
* Will be applied to ALL geometries. Be careful since this can cause
errors when plotting geometries of different types
Returns:
A namedtuple containing:
–
'ax' -> The map axis
'handles' -> All geometry handles which were created in the order they were
drawn
'cbar' -> The colorbar handle if it was drawn
defdrawGeoms(geoms:ogr.Geometry|list[ogr.Geometry]|pd.DataFrame|np.ndarray,# srs: srs_input = 4326,srs:srs_input|None=None,ax:None|matplotlib.axes._axes.Axes|AxHands=None,simplificationFactor:numeric|None=5000,colorBy:str|None=None,figsize:tuple[numeric,numeric]=(12,12),xlim:tuple[numeric,numeric]|None=None,ylim:tuple[numeric,numeric]|None=None,fontsize:int=16,hideAxis:bool=False,cbarTitle=None,vmin=None,vmax=None,cmap="viridis",cbargs:dict|None=None,**mplArgs,)->AxHands:"""Draw geometries onto a matplotlib figure. * Each geometry type is displayed as an appropriate plotting type -> Points/ Multipoints are displayed as points using plt.plot(...) -> Lines/ MultiLines are displayed as lines using plt.plot(...) -> Polygons/ MultiPolygons are displayed as patches using the descartes library * Each geometry can be given its own set of matplotlib plotting parameters Notes ----- This function does not call plt.show() for the final display of the figure. This must be done manually after calling this function. Otherwise plt.savefig(...) can be called to save the output somewhere. Sometimes geometries will disappear because of the simplification procedure. If this happens, the procedure can be avoided by setting simplificationFactor to None. This will take much more memory and will take longer to plot, however Parameters ---------- geoms : ogr.Geometry or [ogr.Geometry, ] or pd.DataFrame The geometries to be drawn * If a DataFrame is given, the function looks for geometries under a columns named 'geom' * plotting arguments can be given by adding a column named 'MPL:****' where '****' stands in for the argument to be added - For geometries that should ignore this argument, set it as None srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs in which to draw each geometry * If not given, longitude/latitude is assumed * Although geometries can be given in any SRS, it is very helpful if they are already provided in the correct SRS ax : matplotlib axis; optional The axis to draw the geometries on * If not given, a new axis is generated and returned simplificationFactor : float; optional The level to which geometries should be simplified. It can be thought of as the number of vertices allowed in either the X or Y dimension across the figure * A higher value means a more detailed plot, but may take longer to draw colorBy : str; optional The column in the geoms DataFrame to color by * Only useful when geoms is given as a DataFrame figsize : (int, int); optional The figure size to create when generating a new axis * If resultign figure looks weird, altering the figure size is your best bet to make it look nicer xlim : (float, float); optional The x-axis limits ylim : (float, float); optional The y-axis limits fontsize : int; optional A base font size to apply to tick marks which appear * Titles and labels are given a size of 'fontsize' + 2 hideAxis : bool; optional Instructs the created axis to hide its boundary * Only useful when generating a new axis cbarTitle : str; optional The title to give to the generated colorbar * If not given, but 'colorBy' is given, the same string for 'colorBy' is used * Only useful when 'colorBy' is given vmin : float; optional The minimum value to color * Only useful when 'colorBy' is given vmax : float; optional The maximum value to color * Only useful when 'colorBy' is given cmap : str or matplotlib ColorMap; optional The colormap to use when coloring * Only useful when 'colorBy' is given cbax : matplotlib axis; optional An explicitly given axis to use for drawing the colorbar * If not given, but 'colorBy' is given, an axis for the colorbar is automatically generated cbargs : dict; optional keyword arguments to pass on when creating the colorbar **mplArgs All other keyword arguments are passed on to the plotting functions called for each geometry * Will be applied to ALL geometries. Be careful since this can cause errors when plotting geometries of different types Returns ------- A namedtuple containing: 'ax' -> The map axis 'handles' -> All geometry handles which were created in the order they were drawn 'cbar' -> The colorbar handle if it was drawn """ifisinstance(ax,AxHands):ax=ax.axifaxisNone:passfig,ax=plt.subplots(figsize=figsize)elifisinstance(ax,matplotlib.axes._axes.Axes):passelse:raiseException("Expected None or matplotlib.axes._axes.Axes object for the 'ax' argument. However, an object of type: "+str(type(ax))+" has been provided.")ifhideAxis:ax.axis("off")ifcolorByisNone:cbax=Noneelse:divider=make_axes_locatable(ax)cbax=divider.append_axes("right",size="2.5%",pad=0.05)ax.tick_params(labelsize=fontsize)# # Be sure we have a listpargs=NoneisFrame=Falseifisinstance(geoms,ogr.Geometry):geoms=[geoms,]elifisinstance(geoms,pd.DataFrame):# We have a DataFrame with plotting argumentsisFrame=Truedata=geoms.drop("geom",axis=1)geoms=geoms["geom"].valuespargs=pd.DataFrame(index=data.index)forcindata.columns:ifnotc[:4]=="MPL:":continuepargs[c[4:]]=data[c]ifpargs.size==0:pargs=Noneelse:# Assume its an iterablegeoms=list(geoms)# Check Geometry SRSifnotsrsisNone:srs=SRS.loadSRS(srs)transformed_geoms=[]forgi,geometryinenumerate(geoms):gsrs=geometry.GetSpatialReference()ifgsrsisNone:continue# Skip it if we don't know it...ifnotgsrs.IsSame(srs):transformed_geoms.append(transform(geoms[gi],srs))else:transformed_geoms.append(geoms[gi])geoms=np.asarray(transformed_geoms)# Apply simplifications if requiredifnotsimplificationFactorisNone:ifxlimisNoneorylimisNone:xMin,yMin,xMax,yMax=1e100,1e100,-1e100,-1e100forgeometryingeoms:_xMin,_xMax,_yMin,_yMax=geometry.GetEnvelope()xMin=min(_xMin,xMin)xMax=max(_xMax,xMax)yMin=min(_yMin,yMin)yMax=max(_yMax,yMax)ifnotxlimisNone:xMin,xMax=xlimifnotylimisNone:yMin,yMax=ylimsimplificationValue=max(xMax-xMin,yMax-yMin)/simplificationFactoroGeoms=geomsgeoms=[]defdoSimplify(g):ng=g.Simplify(simplificationValue)returnngforgeometryinoGeoms:# carefulSimplification=False# if carefulSimplification and "MULTI" in g.GetGeometryName():if"MULTI"ingeometry.GetGeometryName():subgeoms=[]forgiinrange(geometry.GetGeometryCount()):ng=doSimplify(geometry.GetGeometryRef(gi))subgeoms.append(ng)geoms.append(flatten(subgeoms))else:geoms.append(doSimplify(geometry))# Handle color valueifnotcolorByisNone:color_values=data[colorBy].valuesifisinstance(cmap,str):frommatplotlibimportcmcmap=getattr(cm,cmap)ifis_numeric_dtype(arr_or_dtype=color_values.dtype):color_values_without_nan=color_values[~np.isnan(color_values)]ifcolor_values_without_nan.size==0:# All values are NaN: avoid max/min on empty array and division by zero.# Use a default normalized value (e.g. 0.0) for all entries.norm_vals=np.zeros_like(color_values,dtype=float)else:cValMax=color_values_without_nan.max()ifvmaxisNoneelsevmaxcValMin=color_values_without_nan.min()ifvminisNoneelsevmindenom=cValMax-cValMinifdenom==0:# All (non-NaN) values are identical: avoid division by zero.# Map everything to the middle of the colormap.norm_vals=np.full_like(color_values,0.5,dtype=float)else:norm_vals=(color_values-cValMin)/denom_colorVals=[cmap(v)forvinnorm_vals]else:# categorical dataunique_values=pd.unique(color_values)n_values=unique_values.sizecolor_map={val:cmap(i/n_values)fori,valinenumerate(unique_values)}_colorVals=[color_map[v]ifvincolor_mapelse(0,0,0,0)forvincolor_values]# Do Plotting# make patcheshandles=[]ifnotpargsisNone:s=[notvisNoneforvinpargs.iloc[gi]]plotargs=pargs.iloc[gi,s].to_dict()else:plotargs=dict()plotargs.update(mplArgs)forgi,geometryinenumerate(geoms):ifnotcolorByisNone:colorVal=_colorVals[gi]else:colorVal=None# Determine typeifgeometry.GetGeometryName()=="POINT":handles.append(drawPoint(g=geometry,plotargs=plotargs,ax=ax,colorVal=colorVal))elifgeometry.GetGeometryName()=="MULTIPOINT":handles.append(drawMultiPoint(geometry,plotargs,ax,colorVal))elifgeometry.GetGeometryName()=="LINESTRING":handles.append(drawLine(geometry,plotargs,ax,colorVal))elifgeometry.GetGeometryName()=="MULTILINESTRING":handles.append(drawMultiLine(geometry,plotargs,ax,colorVal))elifgeometry.GetGeometryName()=="LINEARRING":handles.append(drawLinearRing(geometry,plotargs,ax,colorVal))elifgeometry.GetGeometryName()=="POLYGON":handles.append(drawPolygon(g=geometry,plotargs=plotargs,ax=ax,colorVal=colorVal))elifgeometry.GetGeometryName()=="MULTIPOLYGON":handles.append(drawMultiPolygon(g=geometry,plotargs=plotargs,ax=ax,colorVal=colorVal))else:msg=("Could not draw geometry of type:",pargs.index[gi],"->",geometry.GetGeometryName(),)warnings.warn(msg,UserWarning)# Add the colorbar, maybeifcolorByisnotNone:norm=Normalize(vmin=cValMin,vmax=cValMax)tmp=dict(cmap=cmap,norm=norm,orientation="vertical")ifnotcbargsisNone:tmp.update(cbargs)color_bar=ColorbarBase(cbax,**tmp)color_bar.ax.tick_params(labelsize=fontsize)color_bar.set_label(colorByifcbarTitleisNoneelsecbarTitle,fontsize=fontsize+2)else:color_bar=None# Do some formattingax.set_aspect("equal")ax.autoscale(enable=True)ifnotxlimisNone:ax.set_xlim(*xlim)ifnotylimisNone:ax.set_ylim(*ylim)# Organize returnifisFrame:returnAxHands(ax,pd.Series(handles,index=data.index),color_bar)else:returnAxHands(ax,handles,color_bar)
The figure size to create when generating a new axis
* If resulting figure looks weird, altering the figure size is your best
bet to make it look nicer
The spacing padding to add between the generated axis and the generated
colorbar axis
* Only useful when generating a new axis
* Only useful when 'colorBy' is given
The title to give to the generated colorbar
* If not given, but 'colorBy' is given, the same string for 'colorBy'
is used
* Only useful when 'colorBy' is given
defdrawImage(matrix:np.ndarray,ax:matplotlib.axis.Axis|None=None,xlim:tuple[float|float]|None=None,ylim:tuple[float|float]|None=None,yAtTop:bool=True,scaling:float=1,fontsize:int=16,hideAxis=False,figsize:tuple[float,float]=(12,12),cbar:bool=True,cbarPadding=0.01,cbarTitle=None,vmin=None,vmax=None,cmap="viridis",cbax=None,cbargs=None,leftMargin=0,rightMargin=0,topMargin=0,bottomMargin=0,**kwargs,)->AxHands:"""Draw a matrix as an image on a matplotlib canvas. Parameters ---------- matrix : numpy.ndarray The matrix data to draw ax : matplotlib.axis.Axis; optional The axis to draw the geometries on * If not given, a new axis is generated and returned xlim : tuple[float, float]; optional The x-axis limits to draw the matrix on ylim : tuple[float, float]; optional The y-axis limits to draw the matrix on yAtTop : bool; optional If True, the first row of data should be plotted at the top of the image scaling : numeric; optional An integer factor by which to scale the matrix before plotting figsize : tuple[float, float]; optional The figure size to create when generating a new axis * If resulting figure looks weird, altering the figure size is your best bet to make it look nicer fontsize : int; optional A base font size to apply to tick marks which appear * Titles and labels are given a size of 'fontsize' + 2 cbar : bool; optional If True, a color bar will be drawn cbarPadding : float; optional The spacing padding to add between the generated axis and the generated colorbar axis * Only useful when generating a new axis * Only useful when 'colorBy' is given cbarTitle : str; optional The title to give to the generated colorbar * If not given, but 'colorBy' is given, the same string for 'colorBy' is used * Only useful when 'colorBy' is given vmin : float; optional The minimum value to color vmax : float; optional The maximum value to color cmap : str or matplotlib ColorMap; optional The colormap to use when coloring cbax : matplotlib axis; optional An explicitly given axis to use for drawing the colorbar cbargs : dict; optional leftMargin : float; optional Additional margin to add to the left of the figure * Before using this, try adjusting the 'figsize' rightMargin : float; optional Additional margin to add to the left of the figure * Before using this, try adjusting the 'figsize' topMargin : float; optional Additional margin to add to the left of the figure * Before using this, try adjusting the 'figsize' bottomMargin : float; optional Additional margin to add to the left of the figure * Before using this, try adjusting the 'figsize' Returns ------- A namedtuple containing: 'ax' -> The map axis 'handles' -> All geometry handles which were created in the order they were drawn 'cbar' -> The colorbar handle if it was drawn """# Create an axis, if neededifisinstance(ax,AxHands):ax=ax.axifaxisNone:newAxis=Trueplt.figure(figsize=figsize)ifnotcbar:# We don't need a colorbarifnothideAxis:leftMargin+=0.07ax=plt.axes([leftMargin,bottomMargin,1-(rightMargin+leftMargin),1-(topMargin+bottomMargin),])cbax=Noneelse:# We need a colorbarrightMargin+=0.08# Add area on the right for colorbar textifnothideAxis:leftMargin+=0.07cbarExtraPad=0.05cbarWidth=0.04ax=plt.axes([leftMargin,bottomMargin,1-(rightMargin+leftMargin+cbarWidth+cbarPadding),1-(topMargin+bottomMargin),])cbax=plt.axes([1-(rightMargin+cbarWidth),bottomMargin+cbarExtraPad,cbarWidth,1-(topMargin+bottomMargin+2*cbarExtraPad),])ifhideAxis:ax.axis("off")else:ax.tick_params(labelsize=fontsize)else:newAxis=False# handle flipped matrixifnotyAtTop:matrix=matrix[::-1,:]# Draw imageifscaling:matrix=scaleMatrix(matrix,scaling,strict=False)ifnot(xlimisNoneandylimisNone):extent=xlim[0],xlim[1],ylim[0],ylim[1]else:extent=Noneh=ax.imshow(matrix,extent=extent,cmap=cmap,vmin=vmin,vmax=vmax,interpolation="none",**kwargs,)# Draw Colorbarifcbar:tmp=dict(cmap=cmap,orientation="vertical")ifnotcbargsisNone:tmp.update(cbargs)ifcbaxisNone:cbar=plt.colorbar(h,ax=ax,**tmp)else:cbar=plt.colorbar(h,cax=cbax)cbar.ax.tick_params(labelsize=fontsize)ifnotcbarTitleisNone:cbar.set_label(cbarTitle,fontsize=fontsize+2)# Do some formattingifnewAxis:ax.set_aspect("equal")ax.autoscale(enable=True)# Done!returnAxHands(ax,h,cbar)
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs of the drawn raster data
* If not given, the raster's internal srs is assumed
* If the drawing resolution does not match the source's inherent
resolution, the source will be warped to the correct format
The resolution of the plotted raster data
* Lower resolution means more pixels to draw and can be a burden on
memory
* If a tuple is given, resolutions in the X and Y direction are expected
* Changing the resolution from the inherent resolution requires a warp
The cutline to limit the drawn data too
* If a string is given, it must be a path to a vector file
* Values outside of the cutline are given the noData value of the raster
* Requires a warp
The figure size to create when generating a new axis
* If resultign figure looks weird, altering the figure size is your best
bet to make it look nicer
The title to give to the generated colorbar
* If not given, but 'colorBy' is given, the same string for 'colorBy'
is used
* Only useful when 'colorBy' is given
defdrawRaster(source:load_raster_input|pd.DataFrame,srs:srs_input|None=None,ax:matplotlib.axes._axes.Axes|None|AxHands=None,resolution:numeric|None=None,cutline=None,figsize:tuple[numeric,numeric]=(12,12),xlim:tuple[numeric,numeric]|None=None,ylim:tuple[numeric,numeric]|None=None,fontsize:int=16,hideAxis:bool=False,cbar:bool=True,cbarTitle:str|None=None,vmin:float|None=None,vmax:float|None=None,cmap="viridis",cbargs=None,noData:numeric|None=None,zorder=0,resampleAlg:Literal["near","bilinear","cubic","cubicspline","lanczos","average","rms","mode","max","min","med","Q1","Q3","sum",]="med",**warp_kwargs,)->AxHands:"""Draw a raster as an image on a matplotlib canvas. Parameters ---------- source : Anything acceptable by loadRaster() The raster datasource to draw srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the drawn raster data * If not given, the raster's internal srs is assumed * If the drawing resolution does not match the source's inherent resolution, the source will be warped to the correct format ax : matplotlib axis; optional The axis to draw the geometries on * If not given, a new axis is generated and returned resolution : numeric or tuple; optional The resolution of the plotted raster data * Lower resolution means more pixels to draw and can be a burden on memory * If a tuple is given, resolutions in the X and Y direction are expected * Changing the resolution from the inherent resolution requires a warp cutline : str or ogr.Geometry; optional The cutline to limit the drawn data too * If a string is given, it must be a path to a vector file * Values outside of the cutline are given the noData value of the raster * Requires a warp figsize : (int, int); optional The figure size to create when generating a new axis * If resultign figure looks weird, altering the figure size is your best bet to make it look nicer xlim : (float, float); optional The x-axis limits ylim : (float, float); optional The y-axis limits fontsize : int; optional A base font size to apply to tick marks which appear * Titles and labels are given a size of 'fontsize' + 2 cbarTitle : str; optional The title to give to the generated colorbar * If not given, but 'colorBy' is given, the same string for 'colorBy' is used * Only useful when 'colorBy' is given vmin : float; optional The minimum value to color * Only useful when 'colorBy' is given vmax : float; optional The maximum value to color * Only useful when 'colorBy' is given cmap : str or matplotlib ColorMap; optional The colormap to use when coloring * Only useful when 'colorBy' is given cbargs : dict; optional noData : numeric; optional Replaces all previous noData values with this value in the output raster. resampleAlg : str, optional The resampleAlg passed on to a call of warp() if needed, by default "med" **kwargs : Passed on to a call to warp() * Determines how the warping is carried out * Consider using 'resampleAlg' or 'workingType' for finer control Returns ------- A namedtuple containing: 'ax' -> The map axis 'handles' -> All geometry handles which were created in the order they were drawn 'cbar' -> The colorbar handle if it was drawn """ifisinstance(ax,AxHands):ax=ax.axifaxisNone:new_main_axis=Truefig,ax=plt.subplots(figsize=figsize)elifisinstance(ax,matplotlib.axes._axes.Axes):new_main_axis=Falseelse:raiseException("Expected None or matplotlib.axes._axes.Axes object for the 'ax' argument. However, an object of type: "+str(type(ax))+" has been provided.")ifhideAxis:ax.axis("off")ifcbarisFalse:cbax=Noneelse:divider=make_axes_locatable(ax)cbax=divider.append_axes("right",size="2.5%",pad=0.05)ax.tick_params(labelsize=fontsize)# Load the raster datasource and check for transformationsource=loadRaster(source)info=rasterInfo(source)ifnot(srsisNoneandresolutionisNoneandcutlineisNoneandxlimisNoneandylimisNone):ifnot(xlimisNoneandylimisNone):bounds=(xlim[0],ylim[0],xlim[1],ylim[1],)else:bounds=NoneifresolutionisNone:xres,yres=None,Noneelse:try:xres,yres=resolutionexcept:xres,yres=resolution,resolutionsource=warp(source,cutline=cutline,pixelHeight=yres,pixelWidth=xres,srs=srs,bounds=bounds,noData=noData,resampleAlg=resampleAlg,**warp_kwargs,)info=rasterInfo(source)# Read the Datadata=extractMatrix(source).astype(float)data[data==info.noData]=np.nan# Draw imageext=(info.xMin,info.xMax,info.yMin,info.yMax,)h=ax.imshow(data,extent=ext,vmin=vmin,vmax=vmax,cmap=cmap,zorder=zorder,interpolation="none",)# Draw ColorbarifcbarisTrue:tmp=dict(cmap=cmap,orientation="vertical")ifcbargsisnotNone:tmp.update(cbargs)cbar_output=plt.colorbar(h,cax=cbax,**tmp)cbar_output.ax.tick_params(labelsize=fontsize)ifcbarTitleisnotNone:cbar_output.set_label(cbarTitle,fontsize=fontsize+2)else:cbar_output=None# Do some formattingax.set_aspect("equal")ax.autoscale(enable=True)ifxlimisnotNone:ax.set_xlim(*xlim)ifylimisnotNone:ax.set_ylim(*ylim)# Done!returnUTIL.AxHands(ax,h,cbar_output)
NOTE:
* The basemap is drawn using the Smopy python package. See here: https://github.com/rossant/smopy
* Be careful to adhere to the usage guidelines of the chosen tile source
- By default, this source is OSM. See here: https://wiki.openstreetmap.org/wiki/Tile_servers
defempty(gtype,srs=None):""" Make a generic OGR geometry of a desired type. *Not for the feint of heart* Parameters ---------- gtpe : str The geometry type to make * Point, MultiPoint, Line, MultiLine, Polygon, MultiPolygon, etc... srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the geometry to create Returns ------- ogr.Geometry """ifnothasattr(ogr,"wkb"+gtype):raiseGeoKitGeomError("Could not find geometry type: "+gtype)geom=ogr.Geometry(getattr(ogr,"wkb"+gtype))ifnotsrsisNone:geom.AssignSpatialReference(SRS.loadSRS(srs))returngeom
Extracts features from a source and clips them to the boundaries of a given geom.
Optionally scales numeric attribute values linearly to the overlapping area share.
An SQL-like where statement to apply to the source
* Feature attribute name do not need quotes
* String values should be wrapped in 'single quotes'
Example: If the source vector has a string attribute called "ISO" and
a integer attribute called "POP", you could use....
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs of the geometries to extract
* If not given, the geom's inherent srs is used
* If srs does not match the inherent srs, all geometries will be
transformed
The min. area share of a polygon that falls either inside or
outside the clipping geom. Allows to deal with imperfect boundary
alignments. 0 means that all clipped geoms, however small they
may be, are considered. Example: If minShare=0.001 (0.1%), a
polygon that overlaps with the clipping geom by 99.92% of its
area will NOT be clipped. If another polygon also at
minShare=0.001 overlaps by only 0.06% of its area, it will be
disregarded. By default 0.001.
defextractAndClipFeatures(source:load_vector_input|pd.DataFrame,geom,where=None,srs=None,onlyGeom=False,indexCol=None,skipMissingGeoms=True,layerName=None,scaleAttrs=None,minShare=0.001,**kwargs,)->pd.DataFrame|pd.Series|Generator:""" Extracts features from a source and clips them to the boundaries of a given geom. Optionally scales numeric attribute values linearly to the overlapping area share. Parameters ---------- source : Anything acceptable by loadVector() The vector data source to read from geom : ogr.Geometry The geometry to search with * All features touching this geometry are extracted and clipped to the geometry boundaries. where : str; optional An SQL-like where statement to apply to the source * Feature attribute name do not need quotes * String values should be wrapped in 'single quotes' Example: If the source vector has a string attribute called "ISO" and a integer attribute called "POP", you could use.... where = "ISO='DEU' AND POP>1000" srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the geometries to extract * If not given, the geom's inherent srs is used * If srs does not match the inherent srs, all geometries will be transformed onlyGeom : bool; optional If True, only feature geometries will be returned indexCol : str; optional The feature identifier to use as the DataFrams's index * Only useful when as DataFrame is True skipMissingGeoms : bool; optional If True, then the parser will not read a feature which are missing a geometry layerName : str; optional The name of the layer to extract from the source vector dataset (only applicable in case of a geopackage). scaleAttrs : str or list, optional attribute names of the source with numeric values. The values will be scaled linearly with the area share of the feature overlapping the geom. minShare : float The min. area share of a polygon that falls either inside or outside the clipping geom. Allows to deal with imperfect boundary alignments. 0 means that all clipped geoms, however small they may be, are considered. Example: If minShare=0.001 (0.1%), a polygon that overlaps with the clipping geom by 99.92% of its area will NOT be clipped. If another polygon also at minShare=0.001 overlaps by only 0.06% of its area, it will be disregarded. By default 0.001. Returns ------- * pandas.DataFrame or pandas.Series """assert0<=minShare<=1,f"minShare must be greater or equal to 0 and less or equal to 1.0"# assert and preprocess input sourceifisinstance(source,pd.DataFrame):# check validity of input dataframeifnot"geom"insource.columns:raiseAttributeError(f"source is given as a pd.DataFrame but has not 'geom' column.")ifnotisinstance(source.geom.iloc[0],ogr.Geometry):raiseTypeError(f"source is given as a pd.DataFrame but value in 'geom' column is not of type osgeo.ogr.Geometry.")# return empty dataframe with empty expected "areaShare" column if no geometries contained since vector cannot be created without geometriesiflen(source)==0:source["areaShare"]=Nonereturnsource# generate a vector from source dataframesource=createVector(source)elifisinstance(source,str)orisinstance(source,pathlib.Path):ifnotos.path.isfile(source):raiseFileNotFoundError(f"source is given as a string but is not an existing filepath: {source}")# load as vector filesource=loadVector(source)elifnotisinstance(source,gdal.Dataset):raiseTypeError(f"source must either be a pd.DataFrame, a gdal.Dataset vector instance or a str formatted shapefile path.")# extract only the overlapping geoms, first define srsifsrsisNone:srs=geom.GetSpatialReference()else:geom=GEOM.transform(geom,toSRS=srs)df=extractFeatures(source=source,geom=geom,where=where,srs=srs,onlyGeom=onlyGeom,indexCol=indexCol,skipMissingGeoms=skipMissingGeoms,layerName=layerName,**kwargs,)ifscaleAttrsisNone:scaleAttrs=[]elifisinstance(scaleAttrs,str):scaleAttrs=[scaleAttrs]else:assertisinstance(scaleAttrs,list),f"scaleAttrs must be a str or a list thereof if not None."iflen(df)>0:for_attrinscaleAttrs:ifnot_attrinlist(df.columns):raiseAttributeError(f"'{_attr}' was given as scaleAttrs but is not an attribute of the source dataframe.")ifnotall([isinstance(_val,numbers.Number)for_valindf[_attr]]):raiseTypeError(f"All values in column '{_attr}' in scaleAttrs must be numeric.")# check if we have any features to clip at alliflen(df)==0:# if not, add the mandatory areaShare column in case that it is not there already and return empty dataframedf["areaShare"]=Nonereturndf# else add the expected areaShare columnassertnot"areaShare"inlist(df.columns),f"source data must not contain a 'areaShare' attribute."df["areaShare"]=1.0# check if we need to clip the geometries at allifdf.geom.iloc[0].GetGeometryName()[:5]=="POINT":# we have only points and no further clipping is neededreturndf# else we need to add an ID column and generate a new vectorassertnot"clippingID"inlist(df.columns),f"source data must not contain a 'clippingID' attribute."df["clippingID"]=range(len(df))_vec=createVector(df)# extract only these features intersected by the outer geom boundaryouter_df=extractFeatures(source=_vec,geom=geom.Boundary(),where=where,srs=srs,indexCol=indexCol,skipMissingGeoms=skipMissingGeoms,layerName=layerName,**kwargs,)del_veciflen(outer_df)==0:# we have no features intersecting with the geom boundary, return all included featuresreturndf.drop(columns="clippingID")# else clip these features that are intersected by the geom_clippedIDs=list()_clippedGeoms=list()_areaShares=list()fori,featinzip(outer_df.clippingID,outer_df.geom):_clipped=feat.Intersection(geom)_areaShare=_clipped.Area()/feat.Area()if_areaShare>=1.0-minShare:# the feature is only touched by the boundary (or protrudes minimally outside the geom edge) -> will not be reducedcontinueelif_areaShare<=0.0+minShare:# the feature is fully outside the geom and only touches the geom boundary (or overlaps only minimally with geom)# -> set clipped feature geometry to np.nan to filter out later_clipped=np.nan_clippedGeoms.append(_clipped)_areaShares.append(_areaShare)_clippedIDs.append(i)iflen(_clippedIDs)==0:# we have not clipped any feature at all, return dfreturndf.drop(columns="clippingID")# else replace the original feature geometries with the clipped ones where needed and add area sharesdf.loc[df.clippingID.isin(_clippedIDs),"geom"]=_clippedGeomsdf.loc[df.clippingID.isin(_clippedIDs),"areaShare"]=_areaSharesfor_attrinscaleAttrs:df[_attr]=df.apply(lambdax:x[_attr]*x.areaShare,axis=1)# drop nan geometries that do not (sufficiently) overlap filter geomdf=df[~df.geom.isna()].reset_index(drop=True)# return the adapted dataframereturndf.drop(columns="clippingID")
An SQL-like where statement to apply to the source
* Feature attribute name do not need quotes
* String values should be wrapped in 'single quotes'
Example: If the source vector has a string attribute called "ISO" and
a integer attribute called "POP", you could use....
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs of the geometries to extract
* If not given, the source's inherent srs is used
* If srs does not match the inherent srs, all geometries will be
transformed
defextractAsDataFrame(source:load_vector_input,indexCol:str|None=None,geom:None|ogr.Geometry=None,where=None,srs:srs_input|None=None,**kwargs,):"""Convenience function calling extractFeatures and structuring the output as a pandas DataFrame. * Geometries are written to a column called 'geom' * Attributes are written to a column of the same name Parameters ---------- source : Anything acceptable by loadVector() The vector datasource to read from indexCol : str; optional The feature identifier to use as the DataFrams's index geom : ogr.Geometry; optional The geometry to search within * All features are extracted which touch this geometry where : str; optional An SQL-like where statement to apply to the source * Feature attribute name do not need quotes * String values should be wrapped in 'single quotes' Example: If the source vector has a string attribute called "ISO" and a integer attribute called "POP", you could use.... where = "ISO='DEU' AND POP>1000" srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the geometries to extract * If not given, the source's inherent srs is used * If srs does not match the inherent srs, all geometries will be transformed Returns ------- pandas.DataFrame """warnings.warn("This function will be removed in favor of geokit.vector.extractFeatures",DeprecationWarning,)returnextractFeatures(source=source,indexCol=indexCol,geom=geom,where=where,srs=srs,**kwargs)
If string -> An SQL-like where statement to apply to the source
If int -> The feature's ID within the dataset
* Feature attribute name do not need quotes
* String values should be wrapped in 'single quotes'
Example: If the source vector has a string attribute called "ISO" and
a integer attribute called "POP", you could use....
(Anything acceptable to geokit.srs.loadSRS(); optional)
–
The srs of the geometries to extract
* If not given, the source's inherent srs is used
* If srs does not match the inherent srs, all geometries will be
transformed
defextractFeature(source:load_vector_input,where:None|str|int=None,geom=None,srs=None,onlyGeom=False,onlyAttr=False,)->UTIL.Feature|ogr.Geometry|dict:"""Convenience function calling extractFeatures which assumes there is only one feature to extract. * Will raise an error if multiple features are found Parameters ---------- source : Anything acceptable by loadVector() The vector datasource to read from where : str, int; optional If string -> An SQL-like where statement to apply to the source If int -> The feature's ID within the dataset * Feature attribute name do not need quotes * String values should be wrapped in 'single quotes' Example: If the source vector has a string attribute called "ISO" and a integer attribute called "POP", you could use.... where = "ISO='DEU' AND POP>1000" geom : ogr.Geometry; optional The geometry to search with * All features are extracted which touch this geometry outputSRS : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the geometries to extract * If not given, the source's inherent srs is used * If srs does not match the inherent srs, all geometries will be transformed onlyGeom : bool; optional If True, only feature geometries will be returned onlyAttr : bool; optional If True, only feature attributes will be returned Returns ------- * If onlyGeom and onlyAttr are False: namedtuple -> (geom, attr) * If onlyGeom is True: ogr.Geometry * If onlyAttr is True: dict """warnings.warn(message="extractFeature is deprecated use extractFeatures instead.",category=DeprecationWarning)ifisinstance(where,int):ds=loadVector(source)lyr=ds.GetLayer()ftr=lyr.GetFeature(where)fGeom=ftr.GetGeometryRef().Clone()fItems=ftr.items().copy()else:getter=_extractFeatures(source=source,geom=geom,where=where,srs=srs,onlyGeom=onlyGeom,onlyAttr=onlyAttr,skipMissingGeoms=False,)# Get first resultres=next(getter)ifonlyGeom:fGeom=reselifonlyAttr:fItems=reselse:(fGeom,fItems)=res# try to get a second resulttry:next(getter)exceptStopIteration:passelse:raiseGeoKitVectorError("More than one feature found")ifnotsrsisNoneandnotonlyAttr:srs=SRS.loadSRS(srs)ifnotfGeom.GetSpatialReference().IsSame(srs):fGeom.TransformTo(srs)# Done!ifonlyGeom:returnfGeomelifonlyAttr:returnfItemselse:returnUTIL.Feature(fGeom,fItems)
Be careful when filtering by a geometry as the extracted features may not
necessarily be IN the given shape
* Sometimes they may only overlap
* Sometimes they are only in the geometry's envelope
* To be sure an extracted geometry fits the selection criteria, you may
still need to do further processing or use extractAndClipFeatures()
An SQL-like where statement to apply to the source
* Feature attribute name do not need quotes
* String values should be wrapped in 'single quotes'
Example: If the source vector has a string attribute called "ISO" and
a integer attribute called "POP", you could use....
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs of the geometries to extract
* If not given, the source's inherent srs is used
* If srs does not match the inherent srs, all geometries will be
transformed
Applies only in combination with given 'geom' filter. If "Touches",
all geometries will be extracted that simply touch the filter
geom. If "Overlaps", geometries to be extracted must overlap (for
lines, this represents an "Intersect") either partially or
completely (i.e. it includes "Within"), and if "CentroidWithin"
the centroid of the extracted geom must be within or on the
filter geom. By default "Touches".
NOTE: When filter geom is a polygon, centroids exactly on the
filter geom boundary will NOT be extracted.
Returns:
* If asPandas is True: pandas.DataFrame or pandas.Series
–
defextractFeatures(source,where=None,geom=None,srs=None,onlyGeom=False,onlyAttr=False,asPandas=True,indexCol=None,skipMissingGeoms=True,layerName=None,spatialPredicate:Literal["Touches","Overlaps","CentroidWithin"]="Touches",)->pd.DataFrame|pd.Series|Generator:"""Creates a generator which extract the features contained within the source. * Iteratively returns (feature-geometry, feature-fields) Note: ----- Be careful when filtering by a geometry as the extracted features may not necessarily be IN the given shape * Sometimes they may only overlap * Sometimes they are only in the geometry's envelope * To be sure an extracted geometry fits the selection criteria, you may still need to do further processing or use extractAndClipFeatures() Parameters ---------- source : Anything acceptable by loadVector() The vector data source to read from geom : ogr.Geometry; optional The geometry to search with * All features are extracted which touch this geometry where : str; optional An SQL-like where statement to apply to the source * Feature attribute name do not need quotes * String values should be wrapped in 'single quotes' Example: If the source vector has a string attribute called "ISO" and a integer attribute called "POP", you could use.... where = "ISO='DEU' AND POP>1000" srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the geometries to extract * If not given, the source's inherent srs is used * If srs does not match the inherent srs, all geometries will be transformed onlyGeom : bool; optional If True, only feature geometries will be returned onlyAttr : bool; optional If True, only feature attributes will be returned asPandas : bool; optional Whether or not the result should be returned as a pandas.DataFrame (when onlyGeom is False) or pandas.Series (when onlyGeom is True) indexCol : str; optional The feature identifier to use as the DataFrams's index * Only useful when as DataFrame is True skipMissingGeoms : bool; optional If True, then the parser will not read a feature which are missing a geometry layerName : str; optional The name of the layer to extract from the source vector dataset (only applicable in case of a geopackage). spatialPredicate : str, optional Applies only in combination with given 'geom' filter. If "Touches", all geometries will be extracted that simply touch the filter geom. If "Overlaps", geometries to be extracted must overlap (for lines, this represents an "Intersect") either partially or completely (i.e. it includes "Within"), and if "CentroidWithin" the centroid of the extracted geom must be within or on the filter geom. By default "Touches". NOTE: When filter geom is a polygon, centroids exactly on the filter geom boundary will NOT be extracted. Returns ------- * If asPandas is True: pandas.DataFrame or pandas.Series * If asPandas is False: generator """# arrange outputifnotasPandas:return_extractFeatures(source=source,geom=geom,where=where,srs=srs,onlyGeom=onlyGeom,onlyAttr=onlyAttr,skipMissingGeoms=skipMissingGeoms,layerName=layerName,spatialPredicate=spatialPredicate,)else:fields=defaultdict(list)fields["geom"]=[]forg,ain_extractFeatures(source=source,geom=geom,where=where,srs=srs,onlyGeom=False,onlyAttr=False,skipMissingGeoms=skipMissingGeoms,layerName=layerName,spatialPredicate=spatialPredicate,):fields["geom"].append(g.Clone())fork,vina.items():fields[k].append(v)df=pd.DataFrame(fields)ifnotindexColisNone:df.set_index(indexCol,inplace=True,drop=False)ifonlyGeom:assertnotonlyAttr,f"onlyGeom cannot be combined with onlyAttr."returndf["geom"]elifonlyAttr:returndf.drop("geom",axis=1)else:returndf
Extract all or part of a raster's band as a numpy matrix.
Note:
Unless one is trying to get the entire matrix from the raster dataset, usage
of this function requires intimate knowledge of the raster's characteristics.
In such a case it is probably easier to use Extent.extractMatrix
The boundary to clip the raster to before mutating
* If given as an Extent, the extent is always cast to the source's
- native srs before mutating
* If given as a tuple, (xMin, yMin, xMax, yMax) is expected
- Units must be in the srs specified by 'boundsSRS'
* This boundary must fit within the boundary of the rasters source
* The boundary is always fitted to the source's grid, so the returned
values do not necessarily match to the boundary which is provided
If True, the matrix will search for no data values and change them to
numpy.nan
* Data type will always result in a float, so be careful with large
matrices
defextractMatrix(source:load_raster_input,bounds=None,boundsSRS:srs_input="latlon",maskBand:bool=False,autocorrect:bool=False,returnBounds:bool=False,)->np.ndarray|tuple[np.ndarray,tuple[float,float,float,float]|None]:"""Extract all or part of a raster's band as a numpy matrix. Note: ----- Unless one is trying to get the entire matrix from the raster dataset, usage of this function requires intimate knowledge of the raster's characteristics. In such a case it is probably easier to use Extent.extractMatrix Parameters ---------- source : Anything acceptable by loadRaster() The raster datasource bounds: tuple or Extent The boundary to clip the raster to before mutating * If given as an Extent, the extent is always cast to the source's - native srs before mutating * If given as a tuple, (xMin, yMin, xMax, yMax) is expected - Units must be in the srs specified by 'boundsSRS' * This boundary must fit within the boundary of the rasters source * The boundary is always fitted to the source's grid, so the returned values do not necessarily match to the boundary which is provided boundsSRS: Anything acceptable to geokit.srs.loadSRS(); optional The srs of the 'bounds' argument * This is ignored if the 'bounds' argument is an Extent object or is None autocorrect : bool; optional If True, the matrix will search for no data values and change them to numpy.nan * Data type will always result in a float, so be careful with large matrices returnBounds : bool; optional If True, return the computed bounds along with the matrix data Returns ------- * If returnBounds is False: numpy.ndarray -> Two dimensional matrix * If returnBounds is True: (numpy.ndarray, tuple) - ndarray is matrix data - tuple is the (xMin, yMin, xMax, yMax) of the computed bounds """sourceDS=loadRaster(source)# BE sure we have a rasterdsInfo=rasterInfo(sourceDS)ifmaskBand:mb=sourceDS.GetMaskBand()else:sourceBand=sourceDS.GetRasterBand(1)# get band# Handle the boundariesifboundsisnotNone:# check for extenttry:isExtent=bounds._whatami=="Extent"except:isExtent=False# Ensure srs is okayifisExtent:bounds=bounds.castTo(dsInfo.srs).fit((dsInfo.dx,dsInfo.dy)).xyXYelse:boundsSRS=SRS.loadSRS(boundsSRS)ifnotdsInfo.srs.IsSame(boundsSRS):bounds=GEOM.boundsToBounds(bounds,boundsSRS,dsInfo.srs)bounds=UTIL.fitBoundsTo(bounds,dsInfo.dx,dsInfo.dy,expand=True,startAtZero=True)# Find offsetsxoff=int(np.round((bounds[0]-dsInfo.xMin)/dsInfo.dx))ifxoff<0:raiseGeoKitRasterError("The given boundary exceeds the raster's xMin value")xwin=int(np.round((bounds[2]-dsInfo.xMin)/dsInfo.dx))-xoffifxwin>dsInfo.xWinSize:raiseGeoKitRasterError("The given boundary exceeds the raster's xMax value")ifdsInfo.yAtTop:yoff=int(np.round((dsInfo.yMax-bounds[3])/dsInfo.dy))ifyoff<0:raiseGeoKitRasterError("The given boundary exceeds the raster's yMax value")ywin=int(np.round((dsInfo.yMax-bounds[1])/dsInfo.dy))-yoffifywin>dsInfo.yWinSize:raiseGeoKitRasterError("The given boundary exceeds the raster's yMin value")else:yoff=int(np.round((bounds[1]-dsInfo.yMin)/dsInfo.dy))ifyoff<0:raiseGeoKitRasterError("The given boundary exceeds the raster's yMin value")ywin=int(np.round((bounds[3]-dsInfo.yMin)/dsInfo.dy))-yoffifywin>dsInfo.yWinSize:raiseGeoKitRasterError("The given boundary exceeds the raster's yMax value")else:xoff=0yoff=0xwin=Noneywin=None# get DataifmaskBand:data=mb.ReadAsArray(xoff=xoff,yoff=yoff,win_xsize=xwin,win_ysize=ywin)else:data=sourceBand.ReadAsArray(xoff=xoff,yoff=yoff,win_xsize=xwin,win_ysize=ywin)ifdsInfo.scaleisnotNoneanddsInfo.scale!=1.0:data=data*dsInfo.scaleifdsInfo.offsetisnotNoneanddsInfo.offset!=0.0:data=data+dsInfo.offset# Correct 'nodata' valuesifautocorrect:noData=sourceBand.GetNoDataValue()data=data.astype(np.float64)data[data==noData]=np.nan# make sure we are returning data in the 'flipped-y' orientationifnotisFlipped(source):data=data[::-1,:]# DoneifreturnBounds:returndata,boundselse:returndata
The raster datasource, can be a filepath, a raster dataset etc., see
RASTER.loadRaster() for details. Alternatively, a list of multiple
such raster datasources.
Coordinates for the points to extract
* All points must be in the same SRS
* !REMEMBER! For lat and lon coordinates, X is lon and Y is lat
(opposite of what you may think...)
The window range (in pixels) to extract the values centered around the
closest raster index to the indicated locations.
* A winRange of 0 will only extract the closest raster value
* A winRange of 1 will extract a window of shape (3,3)
* A winRange of 3 will extract a window of shape (7,7)
Return only the extracted data and omit xOffeset, yOffset, inBounds
Returns:
* If only a single location is given:
–
namedtuple -> (data : The extracted data at the location
xOffset : The X index distance from the location to the
center of the closest raster pixel
yOffset : The Y index distance from the location to the
center of the closest raster pixel
inBounds: Flag for whether or not the location is within
The raster's bounds
)
* If Multiple locations are given:
–
pandas.DataFrame
* Columns are (data, xOffset, yOffset, inBounds)
- See above for column descriptions
* Index is 0...N if 'points' input is not a LocationSet
or returns numpy array if _onlyValues=True
defextractValues(source:load_raster_input|list[load_raster_input],points:tuple[float,float]|ogr.Geometry|list[tuple[float,float]]|list[ogr.Geometry]|Location|LocationSet,pointSRS:srs_input|None=None,winRange:int=0,noDataOkay:bool=True,_onlyValues:bool=False,)->ptValue|pd.DataFrame|numeric|np.ndarray:"""Extracts the value of a raster at a given point or collection of points. Can also extract a window of values if desired. * If the given raster is not in the 'flipped-y' orientation, the result will be automatically flipped Notes ----- Generally speaking, interpolateValues() should be used instead of this function Parameters ---------- source : Anything acceptable by loadRaster() or list The raster datasource, can be a filepath, a raster dataset etc., see RASTER.loadRaster() for details. Alternatively, a list of multiple such raster datasources. points : (X,Y) or [(X1,Y1), (X2,Y2), ...] or Location or LocationSet() Coordinates for the points to extract * All points must be in the same SRS * !REMEMBER! For lat and lon coordinates, X is lon and Y is lat (opposite of what you may think...) pointSRS : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the point to create * If not given, longitude/latitude is assumed * Only useful when 'points' is not a LocationSet winRange : int The window range (in pixels) to extract the values centered around the closest raster index to the indicated locations. * A winRange of 0 will only extract the closest raster value * A winRange of 1 will extract a window of shape (3,3) * A winRange of 3 will extract a window of shape (7,7) noDataOkay: bool If True, an error is raised if a 'noData' value is extracted If False, numpy.nan is inserted whenever a 'noData' value is extracted _onlyValues: bool Return only the extracted data and omit xOffeset, yOffset, inBounds Returns ------- * If only a single location is given: namedtuple -> (data : The extracted data at the location xOffset : The X index distance from the location to the center of the closest raster pixel yOffset : The Y index distance from the location to the center of the closest raster pixel inBounds: Flag for whether or not the location is within The raster's bounds ) * If Multiple locations are given: pandas.DataFrame * Columns are (data, xOffset, yOffset, inBounds) - See above for column descriptions * Index is 0...N if 'points' input is not a LocationSet or returns numpy array if _onlyValues=True """# Be sure we have a raster and srsifnotisinstance(source,list):source=[source]# make sure all source entries can be interpreted by loadRastertry:source=[loadRaster(s)forsinsource]except:raiseTypeError("At least one source cannot be loaded by geokit.raster.loadRaster().")raster_srs=Noneforsinsource:# load file to make sure it can be interpreted by loadRasterifraster_srsisNone:raster_srs=rasterInfo(s).srselse:assertraster_srs.IsSame(rasterInfo(s).srs),"All source entries must have the same SRS."# Ensure we have a list of point geometriespoints_as_list_of_geom=_convertPointsToListOfOGRPoints(points=points,point_input_SRS=pointSRS,output_srs=raster_srs)iflen(points_as_list_of_geom)>1:asSingle=Falseeliflen(points_as_list_of_geom)==1:asSingle=Trueelse:raiseGEOM.GeoKitGeomError("No points have been passed to the function extractValues()")# get the srcs that are actually overlapping with our pointssrc_mapper={}iflen(source)>1:forsinsource:# get bounds and filter indices of points that fall into bounds_bounds=rasterInfo(s).bounds_indices=[ifori,pinenumerate(points_as_list_of_geom)if(_bounds[0]<=p.GetX()<=_bounds[2])and(_bounds[1]<=p.GetY()<=_bounds[3])]# add source as key plus list of overlapped point indicessrc_mapper[s]=_indiceselse:# we have only one source which must be applied to all pointssrc_mapper[source[0]]=list(range(len(points_as_list_of_geom)))# srcs = source * len(points)# iterate over all unique srcs and extract the raster values for the affected points _src by _srcvalues=[np.nan]*len(points_as_list_of_geom)# initialize values with nan for every pointinBounds=[False]*len(points_as_list_of_geom)# initialize inbounds with False for every pointxOffset=[np.nan]*len(points_as_list_of_geom)# initialize offsets with nan for every pointyOffset=[np.nan]*len(points_as_list_of_geom)# initialize offsets with nan for every pointfor_src,_indsinsrc_mapper.items():# get the indices of the points for which data has been extracted already_xtrct=[ifori,vinenumerate(values)ifnp.all(pd.notnull(v))]# deduct them from the inds here (they may have been extracted already from an overlapping _src tile)_inds=list(set(_inds)-set(_xtrct))iflen(_inds)==0:# no indices left to extract from this _srccontinue# get the points with this _src via list indices_points=[points_as_list_of_geom[i]foriin_inds]# get the raster info for this _src_info=rasterInfo(_src)# Get x/y values as numpy arraysx=np.array([pt.GetX()forptin_points])y=np.array([pt.GetY()forptin_points])# Calculate x/y indexesxValues=(x-(_info.xMin+0.5*_info.pixelWidth))/_info.pixelWidthxIndexes=np.round(xValues)_xOffset=xValues-xIndexesif_info.yAtTop:yValues=((_info.yMax-0.5*_info.pixelWidth)-y)/abs(_info.pixelHeight)yIndexes=np.round(yValues)_yOffset=yValues-yIndexeselse:yValues=(y-(_info.yMin+0.5*_info.pixelWidth))/_info.pixelHeightyIndexes=np.round(yValues)_yOffset=-1*(yValues-yIndexes)# Calculate the starts and window sizexStarts=xIndexes-winRangeyStarts=yIndexes-winRangewindow=2*winRange+1_inBounds=xStarts>=0_inBounds=_inBounds&(yStarts>=0)_inBounds=_inBounds&(xStarts+window<=_info.xWinSize)_inBounds=_inBounds&(yStarts+window<=_info.yWinSize)# Read values_values=[]band=_src.GetRasterBand(1)forxi,yi,ibinzip(xStarts,yStarts,_inBounds):ifnotib:data=np.ones((window,window))*np.nanelse:# Open and read from rasterdata=band.ReadAsArray(xoff=xi,yoff=yi,win_xsize=window,win_ysize=window)if(_info.scale!=None)and(_info.scale!=1.0):data=data*_info.scaleif(_info.offset!=None)and(_info.offset!=0.0):data=data+_info.offset# Look for nodataif_info.noDataisnotNone:nodata=data==_info.noDataifnodata.any():ifnoDataOkay:# data will neaed to be a float type to represent a nodata valuedata=data.astype(np.float64)data[nodata]=np.nanelse:raiseGeoKitRasterError("No data values found in extractValues with 'noDataOkay' set to False")# flip if not in the 'flipped-y' orientationifnot_info.yAtTop:data=data[::-1,:]ifwinRange==0:# If winRange is 0, there's no need to return a 2D matrixdata=data[0][0]# Append to values_values.append(data)# check if not _inbounds, then replace values with nanforiinrange(len(_values)):ifnot_inBounds[i]:_values[i]=np.nan*np.ones_like(_values[i])# write _values, _inbounds and offsets into the overall container at the respective indicesfori,v,ib,x,yinzip(_inds,_values,_inBounds,_xOffset,_yOffset):values[i]=vinBounds[i]=ibxOffset[i]=xyOffset[i]=yifnp.isnan(values).any():# we still have NaN values leftmsg=f"WARNING: {np.isnan(values).sum()} of the given points (or extraction windows) exceed/s the source's limits. Values were replaced with nan."warnings.warn(msg,UserWarning)# Done!ifasSingle:# A single point was given, so return a single resultif_onlyValues:returnvalues[0]else:returnptValue(values[0],xOffset[0],yOffset[0],inBounds[0])else:if_onlyValues:returnnp.array(values)else:returnpd.DataFrame(dict(data=values,xOffset=xOffset,yOffset=yOffset,inBounds=inBounds),index=None,)
defextractVerticies(geom):"""Get all vertices found on the geometry as a Nx2 numpy.ndarray."""isMulti="MULTI"ingeom.GetGeometryName()# Check geometry typeif"LINE"ingeom.GetGeometryName():ifisMulti:pts=[]forgiinrange(geom.GetGeometryCount()):pts.append(geom.GetGeometryRef(gi).GetPoints())else:pts=geom.GetPoints()elif"POLYGON"ingeom.GetGeometryName():ifisMulti:pts=[]forgiinrange(geom.GetGeometryCount()):newGeom=geom.GetGeometryRef(gi).GetBoundary()pts.append(extractVerticies(newGeom))else:newGeom=geom.GetBoundary()pts=extractVerticies(newGeom)elif"POINT"ingeom.GetGeometryName():ifisMulti:pts=[]forgiinrange(geom.GetGeometryCount()):pts.append(geom.GetGeometryRef(gi).GetPoints())else:pts=geom.GetPoints()else:raiseGeoKitGeomError("Cannot extract points from geometry ")ifisMulti:out=np.concatenate(pts)else:out=np.array(pts)ifout.shape[1]==3:# This can happen when POINTs are extractedout=out[:,:2]returnout
This function allows to deal with polygons that protrude over the SRS bounds
at +/-180° longitude respectively +/-90° latitude. Polygon areas that exceed
those bounds are either clipped or shifted to the "opposite end of the map"
with shapes at the poles being inverted and shifted by 180° to create a
"fold-over" effect.
geom : osgeo-ogr.Geometry
Geometry to fix.
how : str, optional
The way how to deal with sub shapes extending over the bounds:
'shift' : split off and shift to the "opposite end of the map"
'clip' : clip and remove extending shapes completely
deffixOutOfBoundsGeoms(geom,how="shift"):""" This function allows to deal with polygons that protrude over the SRS bounds at +/-180° longitude respectively +/-90° latitude. Polygon areas that exceed those bounds are either clipped or shifted to the "opposite end of the map" with shapes at the poles being inverted and shifted by 180° to create a "fold-over" effect. geom : osgeo-ogr.Geometry Geometry to fix. how : str, optional The way how to deal with sub shapes extending over the bounds: 'shift' : split off and shift to the "opposite end of the map" 'clip' : clip and remove extending shapes completely """assertisinstance(geom,ogr.Geometry),f"geom must be an osgeo.ogr.Geometry"asserthowin["clip","shift"],f"how must be in 'clip', 'shift'"# get the envelope and srs of original geomenv=geom.GetEnvelope()_srs=geom.GetSpatialReference()assert_srs.IsSame(SRS.loadSRS(4326)),f"SRS must be EPSG:4326"ifenv[0]>=-180andenv[1]<=180andenv[2]>=-90andenv[3]<=90:# the polygon is completely within bounds already, return as isreturngeom# else we need to clip and shift/rotategeom_fixed=copy(geom)# HORIZONTZAL BOUNDSifenv[0]<-180orenv[1]>180:# we need to clip & shift in HORIZONTAL directionbasebox_tripleheight=polygon([(-180,-90*3),(-180,90*3),(180,90*3),(180,-90*3)],srs=4326)geom_fixed=geom_fixed.Intersection(basebox_tripleheight)# clip off outer partsifenv[0]<-180andhow=="shift":# extends over left edge, get the left part, shift and mergeleft_part=shift(basebox_tripleheight,lonShift=-360).Intersection(geom)ifgeom_fixed.IsEmpty():geom_fixed=shift(left_part,lonShift=360)# overwrite when no center partelse:geom_fixed=geom_fixed.Union(shift(left_part,lonShift=360))ifenv[1]>180andhow=="shift":# extends over right edge, get the right part, shift and mergeright_part=shift(basebox_tripleheight,lonShift=+360).Intersection(geom)ifgeom_fixed.IsEmpty():geom_fixed=shift(right_part,lonShift=-360)# overwrite when no center partelse:geom_fixed=geom_fixed.Union(shift(right_part,lonShift=-360))# VERTICAL BOUNDSifenv[2]<-90orenv[3]>90:# we need to clip in VERTICAL direction and rotate (horizontal issue are fixed already)deffold_over_pole(geom):"""Aux function to fold a geometry over the +/-90° latitude line and rotate it."""env=geom.GetEnvelope()center_lon=(env[0]+env[1])/2# x value of center axis of whole geomdeffold_polygon(polygon):"""Function that folds a polygon geometry over 90° lat line."""def_fold_ring(ring):"""Core function for every linear ring."""new_ring=ogr.Geometry(ogr.wkbLinearRing)foriinrange(ring.GetPointCount()):x,y,*z=ring.GetPoint(i)# if needed, mirror y at +/-90° line and flip x value around center axis + shift by 180° (fold over)if-90<=y<=90:# all goody_new=yx_new=xelse:# rotate and flip x values, fold y values_x_new=x+2*(center_lon-x)x_new=(_x_new+180)%360ify>90:y_new=180-yy_new=min(y_new,90.0-1e-6)# avoid that geoms touch each other at the poleelify<-90:y_new=-180-yy_new=max(y_new,-90.0+1e-6)# avoid that geoms touch each other at the polenew_ring.AddPoint(x_new,y_new)# create new ring pointreturnnew_ringnew_geom=ogr.Geometry(ogr.wkbPolygon)outer_ring=_fold_ring(polygon.GetGeometryRef(0))new_geom.AddGeometry(outer_ring)# now do this for every potential other (inner) ringforiinrange(1,polygon.GetGeometryCount()):inner_ring=_fold_ring(polygon.GetGeometryRef(i))new_geom.AddGeometry(inner_ring)returnnew_geomifgeom.GetGeometryName()=="POLYGON":# apply function directlynew_geom=fold_polygon(geom)elifgeom.GetGeometryName()=="MULTIPOLYGON":new_geom=ogr.Geometry(ogr.wkbMultiPolygon)# fold iteratively every single sub polyforiinrange(geom.GetGeometryCount()):sub_geom=geom.GetGeometryRef(i)rotated=fold_polygon(sub_geom)new_geom.AddGeometry(rotated)else:raiseNotImplementedError(f"Geometry type '{geom.GetGeometryName()}' not supported")# assign SRS and returnnew_geom.AssignSpatialReference(geom.GetSpatialReference())returnnew_geombasebox=polygon([(-180,-90),(-180,90),(180,90),(180,-90)],srs=4326)geom_fixed_horizontally=copy(geom_fixed)geom_fixed=geom_fixed.Intersection(basebox)ifenv[2]<-90andhow=="shift":# clip off the bottom part and rotate to the other side (by 180°)bottom_part=shift(basebox,latShift=-180).Intersection(geom_fixed_horizontally)bottom_part_shifted=fold_over_pole(bottom_part)geom_fixed=geom_fixed.Union(bottom_part_shifted)ifenv[3]>90andhow=="shift":# clip off the top part and rotate to the other side (by 180°)top_part=shift(basebox,latShift=+180).Intersection(geom_fixed_horizontally)top_part_shifted=fold_over_pole(top_part)geom_fixed=geom_fixed.Union(top_part_shifted)# assign srs and returngeom_fixed.AssignSpatialReference(_srs)returngeom_fixed
Flatten a list of geometries into a single geometry object.
Combine geometries by iteratively union-ing neighbors (according to index)
* example, given a list of geometries (A,B,C,D,E,F,G,H,I,J):
[ A B C D E F G H I J ]
[ AB CD EF GH IJ ]
[ ABCD EFGH IJ ]
[ ABCDEFGH IJ ]
[ ABCDEFGHIJ ] <- This becomes the resulting geometry
Example:
* A list of Polygons/Multipolygons will become a single Multipolygon
* A list of Linestrings/MultiLinestrings will become a single MultiLinestring
defflatten(geoms):"""Flatten a list of geometries into a single geometry object. Combine geometries by iteratively union-ing neighbors (according to index) * example, given a list of geometries (A,B,C,D,E,F,G,H,I,J): [ A B C D E F G H I J ] [ AB CD EF GH IJ ] [ ABCD EFGH IJ ] [ ABCDEFGH IJ ] [ ABCDEFGHIJ ] <- This becomes the resulting geometry Example: -------- * A list of Polygons/Multipolygons will become a single Multipolygon * A list of Linestrings/MultiLinestrings will become a single MultiLinestring """ifnotisinstance(geoms,list):geoms=list(geoms)try:# geoms is not a list, but it might be iterablegeoms=list(geoms)except:raiseValueError("argument must be a list of geometries")iflen(geoms)==0:returnNone# Begin flatteningwhilelen(geoms)>1:newGeoms=[]forgiinrange(0,len(geoms),2):try:ifnotgeoms[gi].IsValid():warnings.warn("WARNING: Invalid Geometry encountered",UserWarning)ifnotgeoms[gi+1].IsValid():warnings.warn("WARNING: Invalid Geometry encountered",UserWarning)newGeoms.append(geoms[gi].Union(geoms[gi+1]))exceptIndexError:# should only occur when length of geoms is oddnewGeoms.append(geoms[gi])geoms=newGeomsreturngeoms[0]
The scaling factor relating the units of the x & y dimensions to the z
dimension
* If factor is 'latlonToM', the x & y units are assumed to be degrees
(lat & lon) and the z units are assumed to be meters. A factor is then
computed for coordinates at the source's center.
* Example: If x,y units are meters and z units are feet, factor should
be 0.3048
defgradient(source,mode="total",factor=1,asMatrix=False,**kwargs):"""Calculate a raster's gradient and return as a new dataset or simply a matrix. Parameters ---------- source : Anything acceptable by loadRaster() The raster datasource mode : str; optional Determines the type of gradient to compute * Options are.... "total" : Calculates the absolute gradient as a ratio "slope" : Same as 'total' "north-south" : Calculates the "north-facing" gradient as a ratio where negative numbers indicate a south facing gradient "east-west" : Calculates the "east-facing" gradient as a ratio where negative numbers indicate a west facing gradient "aspect" : calculates the gradient's direction in radians (0 is east) "dir" : same as 'aspect' factor : numeric or 'latlonToM' The scaling factor relating the units of the x & y dimensions to the z dimension * If factor is 'latlonToM', the x & y units are assumed to be degrees (lat & lon) and the z units are assumed to be meters. A factor is then computed for coordinates at the source's center. * Example: If x,y units are meters and z units are feet, factor should be 0.3048 asMatrix : bool If True, makes the returned value a matrix If False, makes the returned value a raster dataset **kwargs : All extra key word arguments are passed on to a final call to 'createRaster' * Only useful when 'asMatrix' is True Returns ------- * If 'asMatrix' is True: numpy.ndarray * If 'asMatrix' is False: gdal.Dataset """# Make sure source is a sourcesource=loadRaster(source)# Check modeacceptable=["total","slope","north-south","east-west","dir","ew","ns","aspect",]ifmodenotinacceptable:raiseValueError("'mode' not understood. Must be one of: ",acceptable)# Get the factorsourceInfo=rasterInfo(source)iffactor=="latlonToM":latMid=(sourceInfo.yMax+sourceInfo.yMin)/2R_EARTH=6371000DEGtoRAD=np.pi/180yFactor=R_EARTH*DEGtoRAD# Get arc length in meters/DegreexFactor=R_EARTH*DEGtoRAD*np.cos(latMid*DEGtoRAD)# ditto...else:try:xFactor,yFactor=factorexcept:yFactor=factorxFactor=factor# Calculate gradientarr=extractMatrix(source)ifmodein["north-south","ns","total","slope","dir","aspect"]:ns=np.zeros(arr.shape)ns[1:-1,:]=(arr[2:,:]-arr[:-2,:])/(2*sourceInfo.dy*yFactor)ifmodein["north-south","ns"]:output=nsifmodein["east-west","total","slope","dir","aspect"]:ew=np.zeros(arr.shape)ew[:,1:-1]=(arr[:,:-2]-arr[:,2:])/(2*sourceInfo.dx*xFactor)ifmodein["east-west","ew"]:output=ewifmode=="total"ormode=="slope":output=np.sqrt(ns*ns+ew*ew)ifmode=="dir"ormode=="aspect":output=np.arctan2(ns,ew)# Done!ifasMatrix:returnoutputelse:returncreateRaster(bounds=sourceInfo.bounds,pixelWidth=sourceInfo.dx,pixelHeight=sourceInfo.dy,srs=sourceInfo.srs,data=output,**kwargs,)
defindexToCoord(yi:int|np.ndarray,xi:int|np.ndarray,source=None,asPoint:bool=False,bounds=None,dx=None,dy=None,yAtTop=True,srs=None,)->ogr.Geometry|tuple[int|int]|np.ndarray:"""Convert the index of a raster to coordinate values. Parameters ---------- xi : int The x index * a numpy array of ints is also acceptable yi : int The y index * a numpy array of ints is also acceptable source : Anything acceptable by loadRaster() The contentual raster datasource asPoint : bool Instruct program to return point geometries instead of x,y coordinates Returns ------- * If 'asPoint' is True: ogr.Geometry * If 'asPoint' is False: tuple -> (x,y) coordinates """ifsourceisnotNone:# Get source infoifnotisinstance(source,RasterInfo):source=rasterInfo(source)xMin=source.xMin# xMax = source.xMaxyMin=source.yMinyMax=source.yMaxdx=source.dxdy=source.dyyAtTop=source.yAtTopelse:try:xMin,yMin,xMax,yMax=boundsexcept:xMin,yMin,xMax,yMax=bounds.xyXY# Calculate coordinatesifyAtTop:x=xMin+dx*(xi+0.5)y=yMax-dy*(yi+0.5)else:x=xMin+dx*(xi+0.5)y=yMin+dy*(yi+0.5)# make the outputifasPoint:try:# maybe x and y are iterableoutput=[GEOM.point((xx,yy),srs=srs)forxx,yyinzip(x,y)]exceptTypeError:# x and y should be a single pointoutput=GEOM.point((x,y),srs=srs)else:output=np.column_stack([x,y])# Done!returnoutput
The raster datasource, can be a filepath, a raster dataset etc., see
RASTER.loadRaster() for details. Alternatively, a list of multiple
such raster datasources.
Coordinates for the points to extract
* All points must be in the same SRS
* !REMEMBER! For lat and lon coordinates, X is lon and Y is lat
(opposite of what you may think...)
The interpolation scheme to use
* options are...
"near" - Just gets the nearest value (this is default)
"linear-spline" - calculates a linear spline in between points
"cubic-spline" - calculates a cubic spline in between points
"average" - calculates average across a window
"func" - uses user-provided calculator
A user defined interpolation function
* Only utilized when 'mode' equals "func"
* The function must take three arguments in this order...
- A 2 dimensional data matrix
- A x-index-offset
- A y-index-offset
* See the example below for more information
The window range (in pixels) to extract the values centered around the
closest raster index to the indicated locations.
* A winRange of 0 will only extract the closest raster value
* A winRange of 1 will extract a window of shape (3,3)
* A winRange of 3 will extract a window of shape (7,7)
* Only utilized when 'mode' equals "func"
* All interpolation schemes have a predefined window range which is
appropriate to their use
- near -> 0
- linear-spline -> 2
- cubic-spline -> 4
- average -> 3
- func -> 3
Returns:
* If only a single location is given: numeric
–
namedtuple -> (data : The extracted data at the location
xOffset : The X index distance from the location to the
center of the closest raster pixel
yOffset : The Y index distance from the location to the
center of the closest raster pixel
inBounds: Flag for whether or not the location is within
The raster's bounds
)
* If Multiple locations are given: numpy.ndrray -> (N,)
–
where N is the number of locations
Example:
"Interpolate" according to the median value in a 5x5 window
definterpolateValues(source:load_raster_input,points:tuple[float,float]|ogr.Geometry|list[tuple[float,float]]|list[ogr.Geometry]|Location|LocationSet,pointSRS:srs_input="latlon",mode:Literal["near","linear-spline","cubic-spline","average","func"]="near",func:Callable|None=None,winRange:int|None=None,**kwargs,):"""Interpolates the value of a raster at a given point or collection of points. Supports various interpolation schemes: 'near', 'linear-spline', 'cubic-spline', 'average', or user-defined Parameters ---------- source : Anything acceptable by loadRaster() or list The raster datasource, can be a filepath, a raster dataset etc., see RASTER.loadRaster() for details. Alternatively, a list of multiple such raster datasources. points : (X,Y) or [(X1,Y1), (X2,Y2), ...] or Location or LocationSet() Coordinates for the points to extract * All points must be in the same SRS * !REMEMBER! For lat and lon coordinates, X is lon and Y is lat (opposite of what you may think...) pointSRS : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the point to create * If not given, longitude/latitude is assumed * Only useful when 'points' is not a LocationSet mode : str; optional The interpolation scheme to use * options are... "near" - Just gets the nearest value (this is default) "linear-spline" - calculates a linear spline in between points "cubic-spline" - calculates a cubic spline in between points "average" - calculates average across a window "func" - uses user-provided calculator func - function A user defined interpolation function * Only utilized when 'mode' equals "func" * The function must take three arguments in this order... - A 2 dimensional data matrix - A x-index-offset - A y-index-offset * See the example below for more information winRange : int The window range (in pixels) to extract the values centered around the closest raster index to the indicated locations. * A winRange of 0 will only extract the closest raster value * A winRange of 1 will extract a window of shape (3,3) * A winRange of 3 will extract a window of shape (7,7) * Only utilized when 'mode' equals "func" * All interpolation schemes have a predefined window range which is appropriate to their use - near -> 0 - linear-spline -> 2 - cubic-spline -> 4 - average -> 3 - func -> 3 Returns ------- * If only a single location is given: numeric namedtuple -> (data : The extracted data at the location xOffset : The X index distance from the location to the center of the closest raster pixel yOffset : The Y index distance from the location to the center of the closest raster pixel inBounds: Flag for whether or not the location is within The raster's bounds ) * If Multiple locations are given: numpy.ndrray -> (N,) - where N is the number of locations Example: -------- "Interpolate" according to the median value in a 5x5 window >>> def medianFinder( data, xOff, yOff ): >>> return numpy.median(data) >>> >>> result = interpolateValues( <source>, <points>, mode='func', >>> func=medianFinder, winRange=2) """# Determine what the user probably wants as an outputifisinstance(points,(tuple,ogr.Geometry,Location)):asSingle=True# make points a list of length 1 so that the rest works (will be unpacked later)elifisinstance(points,list):iflen(points)==1:asSingle=Trueelse:asSingle=Falseelifisinstance(points,LocationSet):ifpoints.count==1:asSingle=Trueelse:asSingle=Falseelse:# Assume points is already an iterable of some sortraiseGeoKitRasterError("The following datatype has been provided as point but something different was expected: "+str(type(points)))# Do interpolationifmode=="near":# Simple get the nearest valuewin=0ifwinRangeisNoneelsewinRangevalues=extractValues(source,points,pointSRS=pointSRS,winRange=win,_onlyValues=True,**kwargs)ifasSingleisTrue:result=[values]else:result=valueselifmode=="linear-spline":# use a spline interpolation scheme# setup inputswin=2ifwinRangeisNoneelsewinRangex=np.linspace(-1*win,win,2*win+1)y=np.linspace(-1*win,win,2*win+1)# get raw datavalues=extractValues(source,points,pointSRS=pointSRS,winRange=win,**kwargs)ifasSingleisTrue:assertisinstance(values,ptValue)data_frame_values=pd.DataFrame(dict(dict(data=[values.data],xOffset=values.xOffset,yOffset=values.yOffset,inBounds=values.inBounds),),index=None,)else:data_frame_values=values# Calculate interpolated valuesresult=[]forvindata_frame_values.itertuples(index=False):rbs=RectBivariateSpline(y,x,v.data,kx=1,ky=1)result.append(rbs(v.yOffset,v.xOffset)[0][0])elifmode=="cubic-spline":# use a spline interpolation scheme# setup inputswin=4ifwinRangeisNoneelsewinRangex=np.linspace(-1*win,win,2*win+1)y=np.linspace(-1*win,win,2*win+1)# Get raw datavalues=extractValues(source,points,pointSRS=pointSRS,winRange=win,**kwargs)ifasSingleisTrue:assertisinstance(values,ptValue)data_frame_values=pd.DataFrame(dict(dict(data=[values.data],xOffset=values.xOffset,yOffset=values.yOffset,inBounds=values.inBounds),),index=None,)else:data_frame_values=values# Calculate interpolated valuesresult=[]forvindata_frame_values.itertuples(index=False):rbs=RectBivariateSpline(y,x,v.data)result.append(rbs(v.yOffset,v.xOffset)[0][0])elifmode=="average":# Get the average in a windowwin=3ifwinRangeisNoneelsewinRangevalues=extractValues(source,points,pointSRS=pointSRS,winRange=win,**kwargs)ifasSingleisTrue:assertisinstance(values,ptValue)data_frame_values=pd.DataFrame(dict(dict(data=[values.data],xOffset=values.xOffset,yOffset=values.yOffset,inBounds=values.inBounds),),index=None,)else:data_frame_values=valuesresult=[]forvindata_frame_values.itertuples(index=False):result.append(v.data.mean())elifmode=="func":# Use a general function processoriffuncisNone:raiseGeoKitRasterError("'func' mode chosen, but no func kwargs was given")win=3ifwinRangeisNoneelsewinRangevalues=extractValues(source,points,pointSRS=pointSRS,winRange=win,**kwargs)ifasSingleisTrue:assertisinstance(values,ptValue)data_frame_values=pd.DataFrame(dict(dict(data=[values.data],xOffset=values.xOffset,yOffset=values.yOffset,inBounds=values.inBounds),),index=None,)else:data_frame_values=valuesresult=[]forvindata_frame_values.itertuples(index=False):result.append(func(v.data,v.xOffset,v.yOffset))else:raiseGeoKitRasterError("Interpolation mode not understood: ",mode)# Done!ifasSingle:returnresult[0]else:returnnp.array(result)
defisRaster(source):""" Test if loadRaster fails for the given input. Parameters ---------- source : str The path to the raster file to load Returns ------- bool -> True if the given input is a raster """ifisinstance(source,gdal.Dataset):try:ifsource.GetLayerCount()==0:returnTrue# This should always be true?else:returnFalseexcept:returnFalseelifisinstance(source,str):d=gdal.IdentifyDriver(source)meta=d.GetMetadata()ifmeta.get("DCAP_RASTER",False)=="YES":returnTrueelse:returnFalseelse:returnFalse
defisVector(source):""" Test if loadVector fails for the given input. Parameters ---------- source : str The path to the vector file to load Returns ------- bool -> True if the given input is a vector """ifisinstance(source,gdal.Dataset):ifsource.GetLayerCount()>0:returnTrueelse:returnFalseelifisinstance(source,str):d=gdal.IdentifyDriver(source)meta=d.GetMetadata()ifmeta.get("DCAP_VECTOR",False)=="YES":returnTrueelse:returnFalseelse:returnFalse
defline(points,srs=4326):"""Creates an OGR Line object from a given set of points. Parameters ---------- Points : [(x,y), ], Nx2 numpy.ndarray or list of osgeo.ogr.Geometry points. The points defining the line srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the line to create Returns ------- ogr.Geometry """# Make the complete geometryg=ogr.Geometry(ogr.wkbLineString)ifnotsrsisNone:g.AssignSpatialReference(SRS.loadSRS(srs))# Make the lineifall([isinstance(p,ogr.Geometry)forpinpoints]):# convert points into a list of coordinate tuples in correct srspoints=[(transform(p,toSRS=srs).GetX(),transform(p,toSRS=srs).GetY())forpinpoints][g.AddPoint(x,y)forx,yinpoints]# g.AddGeometry(otr)# Ensure validifnotg.IsValid():raiseGeoKitGeomError("Polygon is invalid")# Done!returng
deflistLayers(source,):"""Returns the layer names for each layer that is stored in a geopackage. Parameters ---------- source : Anything acceptable by loadVector() The vector datasource to read from Returns ------- list A list of layer names for the source geopackage. """layer_names=[]ds=loadVector(source)# Loop over the layers to get their names.foriinrange(ds.GetLayerCount()):name=ds.GetLayer(i).GetName()layer_names.append(name)returnlayer_names
defloadRaster(source:load_raster_input,mode=0)->gdal.Dataset:""" Load a raster dataset from a path to a file on disc. Parameters ---------- source : str or gdal.Dataset * If a string is given, it is assumed as a path to a raster file on disc * If a gdal.Dataset is given, it is assumed to already be an open raster and is returned immediately Returns ------- gdal.Dataset """ifisinstance(source,pathlib.Path):source=str(source)ifisinstance(source,str):ds=gdal.Open(source,mode)else:ds=sourceifdsisNone:raiseGeoKitRasterError("Could not load input dataSource: ",str(source))returnds
The SRS to load. Can be:
* osr.SpatialReference object
* EPSG integer ID
* Standardized SRS string like 'EPSG:4326'
* Common SRS name in SRSCOMMON (e.g. 'europe_laea')
* WKT string
* 'laea' — creates a Lambert Azimuthal Equal Area projection
(optionally centered on a geometry)
defloadSRS(source:srs_input,geom:ogr.Geometry|None=None,**kwargs)->osr.SpatialReference:""" Load a spatial reference system (SRS) from various sources. Parameters ---------- source : str | int | osr.SpatialReference | None The SRS to load. Can be: * osr.SpatialReference object * EPSG integer ID * Standardized SRS string like 'EPSG:4326' * Common SRS name in SRSCOMMON (e.g. 'europe_laea') * WKT string * 'laea' — creates a Lambert Azimuthal Equal Area projection (optionally centered on a geometry) geom : osgeo.ogr.Geometry, optional Geometry to center the projection on if 'laea' is used. **kwargs : Passed directly to `SRSCOMMON.laea(**kwargs)`, allowing customization such as `lon`, `lat`, or `name`. Returns ------- osr.SpatialReference """# Do initial check of sourceifisinstance(source,osr.SpatialReference):returnsourceelifsourceisNone:returnNone# Create an empty SRSsrs=osr.SpatialReference()# Check if source is a stringifisinstance(source,str):src_upper=source.strip().upper()# Handle special LAEA caseifsrc_upper=="LAEA":# merge geom into kwargs if providedifgeomisnotNoneand"geom"notinkwargs:kwargs["geom"]=geomreturnSRSCOMMON.laea(**kwargs)elifhasattr(SRSCOMMON,source):# assume a name for one of the common SRS's was givensrs=SRSCOMMON[source]else:try:# try handling as a standardized epsg or esri etc. codesrs=osr.SpatialReference()_val=srs.SetFromUserInput(source)assert_val==0except:srs.ImportFromWkt(source)# assume a Wkt string was inputelifisinstance(source,int):srs.ImportFromEPSG(source)else:raiseGeoKitSRSError("Unknown srs source type: ",type(source))ifgdal.__version__>="3.0.0":srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)# assert that the srs is valid (may be invalid if e.g. wrong integer codes were passed)ifsrs.Validate()!=0:raiseRuntimeError(f"Created srs is invalid.")returnsrs
defloadVector(x:load_vector_input)->gdal.Dataset:""" Load a vector dataset from a path to a file on disc. Parameters ---------- source : str or gdal.Dataset * If a string is given, it is assumed as a path to a vector file on disc * If a gdal.Dataset is given, it is assumed to already be an open vector and is returned immediately Returns ------- gdal.Dataset """ifisinstance(x,str)orisinstance(x,pathlib.Path):ifnotos.path.exists(x):raiseFileNotFoundError(f"Vector file, directory, or resource not found: {x}")# Since we know the path exists, try to open it.ds=gdal.OpenEx(x)ifdsisNone:# This error now clearly means path exists, but GDAL failed to open it or the file is corrupted.raiseGeoKitVectorError(f"Could not load input dataSource: {x}")elifxisNone:# Raise and error if 'None' is passed as the objectraiseGeoKitVectorError("Input dataSource cannot be None.")elifisinstance(x,gdal.Dataset):# If it is an already-opened GDAL Dataset object.ds=xelse:# Handle any other invalid typeraiseTypeError(f"Invalid input type: Expected str, or gdal.Dataset, got {type(x)}")returnds
defmakeBox(*args,**kwargs):"""Alias for geokit.geom.box(...)."""msg="makeBox will be removed soon. Switch to 'box'"warnings.warn(msg,Warning)returnbox(*args,**kwargs)
defmakeEmpty(*args,**kwargs):"""Alias for geokit.geom.empty(...)."""msg="makeEmpty will be removed soon. Switch to 'empty'"warnings.warn(msg,Warning)returnempty(*args,**kwargs)
defmakeLine(*args,**kwargs):"""Alias for geokit.geom.line(...)."""msg="makeLine will be removed soon. Switch to 'line'"warnings.warn(msg,Warning)returnline(*args,**kwargs)
defmakePoint(*args,**kwargs):"""Alias for geokit.geom.point(...)."""msg="makePoint will be removed soon. Switch to 'point'"warnings.warn(msg,Warning)returnpoint(*args,**kwargs)
defmakePolygon(*args,**kwargs):"""Alias for geokit.geom.polygon(...)."""msg="makePolygon will be removed soon. Switch to 'polygon'"warnings.warn(msg,Warning)returnpolygon(*args,**kwargs)
Process all pixels in a raster according to a given function. The boundaries
of the resulting raster can be changed as long as the new boundaries are within
the scope of the original raster, but the resolution cannot.
The function performing the mutation of the raster's data
* The function will take single argument (a 2D numpy.ndarray)
* The function must return a numpy.ndarray of the same size as the input
* The return type must also be containable within a Float32 (int and
boolean is okay)
* See example below for more info
The boundary to clip the raster to before mutating
* If given as an Extent, the extent is always cast to the source's native
srs before mutating
* If given as a tuple, (xMin, yMin, xMax, yMax) is expected
- Units must be in the srs specified by 'boundsSRS'
* This boundary must fit within the boundary of the rasters source
* The boundary is always fitted to the source's grid, so the returned
values do not necessarily match to the boundary which is provided
If True, then before mutating the matrix extracted from the source will have
pixels equal to its 'noData' value converted to numpy.nan
* Data type will always result in a float, so be careful with large
matrices
A path to an output file
* If output is None, the raster will be created in memory and a dataset
handle will be returned
* If output is given, the raster will be written to disk and nothing will
be returned
(Type, str, or numpy-dtype; optional, default:
None
)
–
If given, forces the processed data to be a particular datatype
* Example
- A python numeric type such as bool, int, or float
- A Numpy datatype such as numpy.uint8 or numpy.float64
- a String such as "Byte", "UInt16", or "Double"
defmutateRaster(source:load_raster_input,processor:Callable|None=None,bounds:tuple[numeric,numeric,numeric,numeric]|None=None,boundsSRS:srs_input="latlon",autocorrect:bool=False,output:str|None=None,dtype:geokit_c_data_types_literal|None=None,**create_raster_kwargs,):"""Process all pixels in a raster according to a given function. The boundaries of the resulting raster can be changed as long as the new boundaries are within the scope of the original raster, but the resolution cannot. Parameters ---------- source : Anything acceptable by loadRaster() The raster datasource processor: function; optional The function performing the mutation of the raster's data * The function will take single argument (a 2D numpy.ndarray) * The function must return a numpy.ndarray of the same size as the input * The return type must also be containable within a Float32 (int and boolean is okay) * See example below for more info bounds: tuple or Extent The boundary to clip the raster to before mutating * If given as an Extent, the extent is always cast to the source's native srs before mutating * If given as a tuple, (xMin, yMin, xMax, yMax) is expected - Units must be in the srs specified by 'boundsSRS' * This boundary must fit within the boundary of the rasters source * The boundary is always fitted to the source's grid, so the returned values do not necessarily match to the boundary which is provided boundsSRS: Anything acceptable to geokit.srs.loadSRS(); optional The srs of the 'bounds' argument * This is ignored if the 'bounds' argument is an Extent object or is None autocorrect : bool; optional If True, then before mutating the matrix extracted from the source will have pixels equal to its 'noData' value converted to numpy.nan * Data type will always result in a float, so be careful with large matrices output : str; optional A path to an output file * If output is None, the raster will be created in memory and a dataset handle will be returned * If output is given, the raster will be written to disk and nothing will be returned dtype : Type, str, or numpy-dtype; optional If given, forces the processed data to be a particular datatype * Example - A python numeric type such as bool, int, or float - A Numpy datatype such as numpy.uint8 or numpy.float64 - a String such as "Byte", "UInt16", or "Double" **kwargs: * All kwargs are passed on to a call to createRaster() Example: -------- If you wanted to assign suitability factors based on a raster containing integer identifiers >>> def calcSuitability( data ): >>> # create an output matrix >>> outputMatrix = numpy.zeros( data.shape ) >>> >>> # do the processing >>> outputMatrix[ data == 1 ] = 0.1 >>> outputMatrix[ data == 2 ] = 0.2 >>> outputMatrix[ data == 10] = 0.4 >>> outputMatrix[ np.logical_and(data > 15, data < 20) ] = 0.5 >>> >>> # return the output matrix >>> return outputMatrix >>> >>> result = processRaster( <source-path>, processor=calcSuitability ) """# open the dataset and get SRSworkingDS=loadRaster(source)# Get ds infodsInfo=rasterInfo(workingDS)# Read data into arraysourceData,bounds=extractMatrix(source,bounds=bounds,boundsSRS=boundsSRS,autocorrect=autocorrect,returnBounds=True,)workingExtent=dsInfo.boundsif(boundsisNone)elsebounds# Perform processingifprocessor:processedData=processor(sourceData)else:processedData=sourceData# Ensure returned matrix is okayifprocessedData.shape!=sourceData.shape:raiseGeoKitRasterError("Processed matrix does not have the correct shape \nIs {0}\nShould be {1}",format(processedData.shape,sourceData.shape),)delsourceDatalist_of_numbers=[processedData.min(),processedData.max()]minimum_gdal_type_list=[str(processedData.dtype)]ifisinstance(dtype,str):minimum_gdal_type_list.append(dtype)gdal_data_string=MinimumCDataTypeHandler.get_valid_gdal_data_type_as_string(list_of_numbers=list_of_numbers,minimum_gdal_type_list=minimum_gdal_type_list)# Create an output rasterifoutputisNone:returnUTIL.quickRaster(dy=dsInfo.dy,dx=dsInfo.dx,bounds=workingExtent,dtype=gdal_data_string,srs=dsInfo.srs,data=processedData,**create_raster_kwargs,)else:createRaster(pixelHeight=dsInfo.dy,pixelWidth=dsInfo.dx,bounds=workingExtent,srs=dsInfo.srs,data=processedData,output=output,dtype=gdal_data_string,**create_raster_kwargs,)returnoutput
Process a vector dataset according to an arbitrary function.
Note:
If this is called without a processor, it simply clips the vector source to
the selection criteria given by 'geom' and 'where' as well as translates
all geometries to 'srs'
The function which mutates each feature
* If no function is given, this stage is skipped
* The function will take 1 arguments: a pandas.Series object containing
one 'geom' key indicating the geometry and the other keys indicating
attributes
* The function must return something understandable by pd.Series
containing the geometry under the index 'geom' and any other keys
- These will be used to update the old geometries and attributes
* The attribute dictionary's values should only be numerics and strings
* See example below for more info
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs of the vector to create
* If not given, the source's inherent srs is used
* If the given SRS is different from the source's SRS, all feature
geometries will be cast to the given SRS before processing
An SQL-like where statement to apply to the source
* Feature attribute name do not need quotes
* String values should be wrapped in 'single quotes'
Example: If the source vector has a string attribute called "ISO" and
a integer attribute called "POP", you could use....
A dictionary specifying the datatype of each attribute when written into
the final dataset
* Options are defined from ogr.OFT[...]
- ex. Integer, Real, String
* The ogrType() function can be used to map typical python and numpy types
to appropriate ogr types
A path on disk to create the output vector
* If output is None, the vector dataset will be created in memory
* Assumed to be of "ESRI Shapefile" format
* Will create a number of files with different extensions
If True, the old attributes will be kept in the output vector
* Unless they are over written by the processor
If False, only the newly specified attributes are kept
Say you had a vector source which contains point geometries, and where each
feature also had an float-attribute called "value". You want to create a new
vector set wherein you have circle geometries at the same locations as the
original points and whose radius is equal to the original features' "value"
attribute. Furthermore, let's say you only want to do this for feature's who's
"value" attribute is greater than zero and less than 10. Do as such...
def growPoints( row ):
# Create a new geom
newGeom = row.geom.Buffer(row.radius)
# Return the new geometry/attribute set
return { 'geom':newGeom }
result = processVector( , where="value>0 AND value<10",
processor=growPoints )
defmutateVector(source,processor=None,srs=None,geom=None,where=None,fieldDef=None,output=None,keepAttributes=True,_slim:bool=False,**create_vector_kwargs,)->None|gdal.Dataset|str:"""Process a vector dataset according to an arbitrary function. Note: ----- If this is called without a processor, it simply clips the vector source to the selection criteria given by 'geom' and 'where' as well as translates all geometries to 'srs' Parameters ---------- source : Anything acceptable by loadVector() The vector datasource to read from processor : function; optional The function which mutates each feature * If no function is given, this stage is skipped * The function will take 1 arguments: a pandas.Series object containing one 'geom' key indicating the geometry and the other keys indicating attributes * The function must return something understandable by pd.Series containing the geometry under the index 'geom' and any other keys - These will be used to update the old geometries and attributes * The attribute dictionary's values should only be numerics and strings * See example below for more info srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the vector to create * If not given, the source's inherent srs is used * If the given SRS is different from the source's SRS, all feature geometries will be cast to the given SRS before processing geom : ogr.Geometry; optional The geometry to search within * All features are extracted which touch this geometry where : str; optional An SQL-like where statement to apply to the source * Feature attribute name do not need quotes * String values should be wrapped in 'single quotes' Example: If the source vector has a string attribute called "ISO" and a integer attribute called "POP", you could use.... where = "ISO='DEU' AND POP>1000" fieldDef : dict; optional A dictionary specifying the datatype of each attribute when written into the final dataset * Options are defined from ogr.OFT[...] - ex. Integer, Real, String * The ogrType() function can be used to map typical python and numpy types to appropriate ogr types output : str; optional A path on disk to create the output vector * If output is None, the vector dataset will be created in memory * Assumed to be of "ESRI Shapefile" format * Will create a number of files with different extensions keepAttributes : bool; optional If True, the old attributes will be kept in the output vector * Unless they are over written by the processor If False, only the newly specified attributes are kept _slim: bool Returns ------- * If 'output' is None: gdal.Dataset * If 'output' is given: None Example: -------- Say you had a vector source which contains point geometries, and where each feature also had an float-attribute called "value". You want to create a new vector set wherein you have circle geometries at the same locations as the original points and whose radius is equal to the original features' "value" attribute. Furthermore, let's say you only want to do this for feature's who's "value" attribute is greater than zero and less than 10. Do as such... >>> def growPoints( row ): >>> # Create a new geom >>> newGeom = row.geom.Buffer(row.radius) >>> >>> # Return the new geometry/attribute set >>> return { 'geom':newGeom } >>> >>> result = processVector( <source-path>, where="value>0 AND value<10", >>> processor=growPoints ) >>> """# Extract filtered featuresgeoms=extractFeatures(source,geom=geom,where=where,srs=srs)ifgeoms.size==0:returnNone# Hold on to the SRS in case we need itifsrsisNone:vecds=loadVector(source)veclyr=vecds.GetLayer()srs=veclyr.GetSpatialRef()# Do processingifnotprocessorisNone:result=geoms.apply(lambdax:pd.Series(processor(x)),axis=1)ifkeepAttributes:forcinresult.columns:geoms[c]=result[c].valueselse:geoms=resultifnot"geom"ingeoms:raiseGeoKitVectorError("There is no 'geom' field in the resulting vector table")# make sure the geometries have an srsifnotgeoms.geom[0].GetSpatialReference():srs=SRS.loadSRS(srs)geoms.geom.apply(lambdax:x.AssignSpatialReference(srs))# Create a new shapefile from the resultsif_slimisTrue:returncreateVector(geoms.geom)else:returncreateVector(geoms,srs=srs,output=output,fieldDef=fieldDef,**create_vector_kwargs)
defogrType(s):"""Tries to determine the corresponding OGR type according to the input."""ifisinstance(s,str):ifhasattr(ogr,s):returnselifs.lower()in_ogrStrToType:return_ogrStrToType[s.lower()]elifhasattr(ogr,"OFT%s"%s):return"OFT%s"%sreturn"OFTString"elifsisstr:return"OFTString"elifisinstance(s,np.dtype):returnogrType(str(s))elifisinstance(s,np.generic):returnogrType(s.dtype)elifsisbool:return"OFTInteger"elifsisint:return"OFTInteger64"elifisinstance(s,int):return_ogrIntToType[s]elifsisfloat:return"OFTReal"elifisinstance(s,Iterable):returnogrType(s[0])raiseValueError("OGR type could not be determined")
Partition a Polygon into some number of pieces whose areas should be close
to the targetArea.
WARNING: Not tested for several version. Will probably be removed later
Inputs:
geom : The geometry to partition
- a single ogr Geometry object of POLYGON type
targetArea - float : The ideal area of each partition
* Most of the geometries will be around this area, but they can also be anywhere in the range 0 and 2x
growStep - float : The incremental buffer to add while searching for a suitable partition
* Choose carefully!
- A large growStep will make the algorithm run faster
- A small growStep will produce a more accurate result
* If no growStep is given, a decent one will be calculated
defpartition(geom,targetArea,growStep=None,_startPoint=0):"""Partition a Polygon into some number of pieces whose areas should be close to the targetArea. WARNING: Not tested for several version. Will probably be removed later Inputs: geom : The geometry to partition - a single ogr Geometry object of POLYGON type targetArea - float : The ideal area of each partition * Most of the geometries will be around this area, but they can also be anywhere in the range 0 and 2x growStep - float : The incremental buffer to add while searching for a suitable partition * Choose carefully! - A large growStep will make the algorithm run faster - A small growStep will produce a more accurate result * If no growStep is given, a decent one will be calculated """ifgrowStepisNone:growStep=np.sqrt(targetArea/np.pi)/2# Be sure we are working with a polygonifgeom.GetGeometryName()=="POLYGON":passelifgeom.GetGeometryName()=="MULTIPOLYGON":results=[]forgiinrange(geom.GetGeometryCount()):tmpResults=partition(geom.GetGeometryRef(gi),targetArea,growStep)results.extend(tmpResults)returnresultselse:raiseGeoKitGeomError("Geometry is not a polygon or multipolygon object")# Check the geometry's sizegArea=geom.Area()ifgArea<1.5*targetArea:return[geom.Clone(),]# Find the most starting boundary coordinateboundary=geom.Boundary()ifboundary.GetGeometryCount()==0:coords=boundary.GetPoints()else:coords=boundary.GetGeometryRef(0).GetPoints()y=np.array([c[1]forcincoords])x=np.array([c[0]forcincoords])if_startPoint==0:# start from the TOP-LEFTyStart=y.max()xStart=x[y==yStart].min()elif_startPoint==1:# Start from the LEFT-TOPxStart=x.min()yStart=y[x==xStart].max()elif_startPoint==2:# Start from the RIGHT-TOPxStart=x.max()yStart=y[x==xStart].max()elif_startPoint==3:# Start from the BOT-RIGHTyStart=y.min()xStart=x[y==yStart].max()elif_startPoint==4:# Start from the BOT-LEFTyStart=y.min()xStart=x[y==yStart].min()elif_startPoint==5:# Start from the LEFT-BOTxStart=x.min()yStart=y[x==xStart].min()else:raiseGeoKitGeomError("Start point failure. There may be an infinite loop in one of the geometries")start=point(xStart,yStart,srs=geom.GetSpatialReference())# start searchingtmp=start.Buffer(float(growStep))tmp.Simplify(growStep)searchGeom=tmp.Intersection(geom)sgArea=searchGeom.Area()ifgArea<2*targetArea:# use a slightly smalled target area when the whole geometry# is close to twice the target area in order to increase the# likelihood of a usable leftoverworkingTarget=0.9*targetAreaelse:workingTarget=targetAreawhilesgArea<workingTarget:tmp=searchGeom.Buffer(float(growStep))tmp.Simplify(growStep)newGeom=tmp.Intersection(geom)newArea=newGeom.Area()ifnewArea>workingTarget*1.1:dA=(newArea-sgArea)/growStepweightedGrowStep=(workingTarget-sgArea)/dAtmp=start.Buffer(float(weightedGrowStep))tmp.Simplify(growStep)searchGeom=tmp.Intersection(geom)breaksearchGeom=newGeomsgArea=newArea# fix the search geometry# - For some reason the searchGeometry will sometime create a geometry and a linestring,# in these cases the real geometry was always the second item...ifsearchGeom.GetGeometryName()=="GEOMETRYCOLLECTION":forgiinrange(searchGeom.GetGeometryCount()):g=searchGeom.GetGeometryRef(gi)ifg.GetGeometryName()=="POLYGON":searchGeom=g.Clone()break# Check the left over geometry, maybe some poops have been created and we need to glue them togetheroutputGeom=searchGeom.Simplify(growStep/20)geomToDo=[]leftOvers=geom.Difference(searchGeom)ifleftOvers.GetGeometryName()=="MULTIPOLYGON":forgiinrange(leftOvers.GetGeometryCount()):leftOver=leftOvers.GetGeometryRef(gi)ifleftOver.Area()<targetArea*0.5:outputGeom=outputGeom.Union(leftOver)else:geomToDo.append(leftOver)elifleftOvers.GetGeometryName()=="POLYGON":geomToDo.append(leftOvers)else:raiseGeoKitGeomError("FATAL ERROR: Difference did not result in a polygon")# make an output arrayifoutputGeom.Area()<targetArea*1.5:output=[outputGeom]else:# the search geom plus some (or maybe all) of the left over poops total an area which is too large,# so it will need recomputing. But in order to decrease the liklihhod of an infinite loop,# use a difference starting point than the one used before# - This will loop a maximum of 6 times before an exception is raisedoutput=partition(outputGeom,targetArea,growStep,_startPoint+1)forgingeomToDo:tmpOutput=partition(g,targetArea,growStep)output.extend(tmpOutput)# Done!returnoutput
defpoint(*args,srs:srs_input="latlon"):"""Make a simple point geometry. Parameters ---------- *args : numeric, numeric or (numeric, numeric) The X and Y coordinate of the point to create srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the point to create * If not given, longitude/latitude is assumed * srs MUST be given as a keyword argument Returns ------- ogr.Geometry Example: -------- point(x, y [,srs]) point( (x, y) [,srs] ) """iflen(args)==1:x,y=args[0]eliflen(args)==2:x=args[0]y=args[1]else:raiseGeoKitGeomError('Too many positional inputs. Did you mean to specify "srs="?')"""make a point geometry from given coordinates (x,y) and srs"""pt=ogr.Geometry(ogr.wkbPoint)pt.AddPoint(float(x),float(y))ifnotsrsisNone:pt.AssignSpatialReference(SRS.loadSRS(srs))returnpt
([(x,y), ] or [ogr.Geometry, ] or Nx2 numpy.ndarray, default:
()
)
–
The inner edges of the polygon
* Each input forms a single edge
* Inner rings cannot interset the outer ring or one another
* NOTE! For proper drawing in matplotlib, inner rings must be given in
the opposite orientation as the outer ring (clockwise vs
counterclockwise)
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
'default'
)
–
The srs of the polygon to create. By default "default", i.e. if
point geometries are passed, srs will be extracted from first point of outer ring,
if points are passed as (x, y) tuples, EPSG:4326 will be assigned
by default unless given otherwise. If given as None, no srs will be assigned
defpolygon(outerRing,*args,srs:srs_input|None="default"):"""Creates an OGR Polygon object from a given set of points. Parameters ---------- outerRing : [(x,y), ] or [ogr.Geometry, ] or Nx2 numpy.ndarray The polygon's outer edge *args : [(x,y), ] or [ogr.Geometry, ] or Nx2 numpy.ndarray The inner edges of the polygon * Each input forms a single edge * Inner rings cannot interset the outer ring or one another * NOTE! For proper drawing in matplotlib, inner rings must be given in the opposite orientation as the outer ring (clockwise vs counterclockwise) srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the polygon to create. By default "default", i.e. if point geometries are passed, srs will be extracted from first point of outer ring, if points are passed as (x, y) tuples, EPSG:4326 will be assigned by default unless given otherwise. If given as None, no srs will be assigned Returns ------- ogr.Geometry Example: -------- Make a diamond cut out of a box... box = [(-2,-2), (-2,2), (2,2), (2,-2), (-2,-2)] diamond = [(0,1), (-0.5,0), (0,-1), (0.5,0), (0,1)] geom = polygon( box, diamond ) """# check if we have all point geometriespointGeometries=all([isinstance(_p,ogr.Geometry)for_pinouterRing])ifsrs=="default":ifpointGeometries:# we have geometries, set srs to the srs of the first outer ring pointsrs=outerRing[0].GetSpatialReference()else:# set srs to EPSG:4326 as standardsrs=SRS.loadSRS(4326)elifsrsisNone:passelse:srs=SRS.loadSRS(srs)# Make the complete geometryg=ogr.Geometry(ogr.wkbPolygon)ifnotsrsisNone:g.AssignSpatialReference(srs)# Make the outer ringotr=ogr.Geometry(ogr.wkbLinearRing)ifnotsrsisNone:otr.AssignSpatialReference(srs)# convert to tuples if we have point geometries at handifpointGeometries:outerRing=[(_p.GetX(),_p.GetY())for_pinouterRing]else:assertall([isinstance(x,tuple)forxinouterRing]),(f"All outerRing entries must be (x,y) tuples in given (or default) srs.")[otr.AddPoint(float(x),float(y))forx,yinouterRing]g.AddGeometry(otr)# Make the inner rings (maybe)forinnerRinginargs:# extract (x, y) tuples first if neededifall([isinstance(_p,ogr.Geometry)for_pininnerRing]):innerRing=[(_p.GetX(),_p.GetY())for_pininnerRing]tmp=ogr.Geometry(ogr.wkbLinearRing)ifnotsrsisNone:tmp.AssignSpatialReference(srs)[tmp.AddPoint(float(x),float(y))forx,yininnerRing]g.AddGeometry(tmp)# Make sure geom is validg.CloseRings()ifnotg.IsValid():raiseGeoKitGeomError("Polygon is invalid")# Done!returng
((xMin, yMin, xMax, yMax) or Extent, default:
None
)
–
Determines the boundary context for the given mask and will scale
the resulting geometry's coordinates accordingly
* If a boundary is not given, the geometry coordinates will
correspond to the mask's indices
* If the boundary is given as an Extent object, an srs input is not
required
If True, shrink all geoms by a tiny amount in order to avoid geometry
overlapping issues
* The total amount shrunk should be very very small
* Generally this should be left as True unless it is ABSOLUTELY
necessary to maintain the same area
defpolygonizeMask(mask,bounds=None,srs=None,flat=True,shrink=True):"""Create a geometry set from a matrix mask. Each True-valued group of pixels will be converted to a geometry Parameters ---------- mask : matrix_like The mask which will be turned into a geometry set * Must be 2 dimensional * Must be boolean type * True values are interpreted as 'in the geometry' bounds : (xMin, yMin, xMax, yMax) or geokit.Extent Determines the boundary context for the given mask and will scale the resulting geometry's coordinates accordingly * If a boundary is not given, the geometry coordinates will correspond to the mask's indices * If the boundary is given as an Extent object, an srs input is not required srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the geometries to create flat : bool If True, flattens the resulting geometries into a single geometry shrink : bool If True, shrink all geoms by a tiny amount in order to avoid geometry overlapping issues * The total amount shrunk should be very very small * Generally this should be left as True unless it is ABSOLUTELY necessary to maintain the same area Returns ------- If 'flat' is True: ogr.Geometry else: [ogr.Geometry, ] """# Make sure we have a boolean numpy matrixifnotisinstance(mask,np.ndarray):mask=np.array(mask)ifnot(mask.dtype==boolormask.dtype==np.uint8):raiseGeoKitGeomError("Mask must be a 2D boolean numpy ndarray")# Do vectorizationresult=polygonizeMatrix(matrix=mask,bounds=bounds,srs=srs,flat=flat,shrink=shrink,_raw=True)[0]ifflat:result=result[0]# Done!returnresult
((xMin, yMin, xMax, yMax) or Extent, default:
None
)
–
Determines the boundary context for the given matrix and will scale
the resulting geometry's coordinates accordingly
* If a boundary is not given, the geometry coordinates will
correspond to the mask's indices
* If the boundary is given as an Extent object, an srs input is not
required
If True, shrink all geoms by a tiny amount in order to avoid geometry
overlapping issues
* The total amount shrunk should be very very small
* Generally this should be left as True unless it is ABSOLUTELY
necessary to maintain the same area
return a tuple with with two lists instead of a data frame. The first
list contains the The contiguous-valued geometries and the second list the value
the value for each geometry
Returns:
pandas.DataFrame -> With columns:
–
'geom' -> The contiguous-valued geometries
'value' -> The value for each geometry
| tuple[ list | list ]
–
The first list contains the contiguous-valued geometries. The seconds
–
defpolygonizeMatrix(matrix:np.ndarray,bounds:tuple[numeric,numeric,numeric,numeric,]|None=None,srs:srs_input|None=None,flat:bool=False,shrink:bool=True,_raw:bool=False,)->pd.DataFrame|tuple[list,list]:"""Create a geometry set from a matrix of integer values. Each unique-valued group of pixels will be converted to a geometry Parameters ---------- matrix : matrix_like The matrix which will be turned into a geometry set * Must be 2 dimensional * Must be integer or boolean type bounds : (xMin, yMin, xMax, yMax) or geokit.Extent Determines the boundary context for the given matrix and will scale the resulting geometry's coordinates accordingly * If a boundary is not given, the geometry coordinates will correspond to the mask's indices * If the boundary is given as an Extent object, an srs input is not required srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs context for the given matrix and of the geometries to create flat : bool If True, flattens the resulting geometries which share a contiguous matrix value into a single geometry object shrink : bool If True, shrink all geoms by a tiny amount in order to avoid geometry overlapping issues * The total amount shrunk should be very very small * Generally this should be left as True unless it is ABSOLUTELY necessary to maintain the same area _raw: bool return a tuple with with two lists instead of a data frame. The first list contains the The contiguous-valued geometries and the second list the value the value for each geometry Returns ------- pandas.DataFrame -> With columns: 'geom' -> The contiguous-valued geometries 'value' -> The value for each geometry | tuple[ list | list ] The first list contains the contiguous-valued geometries. The seconds list contains the value for each geometry. """# Make sure we have a boolean numpy matrixifnotisinstance(matrix,np.ndarray):matrix=np.array(matrix)ifmatrix.dtype==boolormatrix.dtype==np.uint8:dtype="GDT_Byte"elifnp.issubdtype(matrix.dtype,np.integer):dtype="GDT_Int32"else:raiseGeoKitGeomError("matrix must be a 2D boolean or integer numpy ndarray")# Make boundaries if not givenifboundsisNone:# bounds in xMin, yMin, xMax, yMaxbounds=(0,0,matrix.shape[1],matrix.shape[0])pixelHeight=1pixelWidth=1try:# first try for a tuplexMin,yMin,xMax,yMax=boundsexcept:# next assume the user gave an extent objecttry:xMin,yMin,xMax,yMax=bounds.xyXYsrs=bounds.srsexcept:raiseGeoKitGeomError("Could not understand 'bounds' input")pixelHeight=(yMax-yMin)/matrix.shape[0]pixelWidth=(xMax-xMin)/matrix.shape[1]ifnotsrsisNone:srs=SRS.loadSRS(srs)# Make a raster dataset and pull the band/maskBand objects# used 'round' instead of 'int' because this matched GDAL behavior bettercols=int(round((xMax-xMin)/pixelWidth))rows=int(round((yMax-yMin)/abs(pixelHeight)))originX=xMinoriginY=yMax# Always use the "Y-at-Top" orientation# Open the driverdriver=gdal.GetDriverByName("Mem")# create a raster in memoryraster=driver.Create("",cols,rows,1,getattr(gdal,dtype))ifrasterisNone:raiseGeoKitGeomError("Failed to create temporary raster")raster.SetGeoTransform((originX,abs(pixelWidth),0,originY,0,-1*abs(pixelHeight)))# Set the SRSifnotsrsisNone:rasterSRS=SRS.loadSRS(srs)raster.SetProjection(rasterSRS.ExportToWkt())# Set data into bandband=raster.GetRasterBand(1)band.SetNoDataValue(0)band.WriteArray(matrix)band.FlushCache()raster.FlushCache()# rasDS = createRaster(bounds=bounds, data=matrix, noDataValue=0, pixelWidth=pixelWidth, pixelHeight=pixelHeight, srs=srs)# Do a polygonizerasBand=raster.GetRasterBand(1)maskBand=rasBand.GetMaskBand()vecDS=gdal.GetDriverByName("Memory").Create("",0,0,0,gdal.GDT_Unknown)vector_layer=vecDS.CreateLayer("mem",srs=srs)field=ogr.FieldDefn("DN",ogr.OFTInteger)vector_layer.CreateField(field)# Polygonize geometryresult=gdal.Polygonize(rasBand,maskBand,vector_layer,0)ifresult!=0:raiseGeoKitGeomError("Failed to polygonize geometry")# Check how many features were createdfeature_count=vector_layer.GetFeatureCount()iffeature_count==0:# raise GlaesError("No features in created in temporary layer")message="No features created in temporary layer"warnings.warn(message,UserWarning)if_rawisTrue:return([],[])elif_rawisFalse:returnpd.DataFrame(dict(geom=[],value=[]))else:raiseException("_raw is suped to be True or False but is: "+str(type(_raw)))# Extract geometries and valuesgeoms=[]rid=[]foriinrange(feature_count):ftr=vector_layer.GetFeature(i)geoms.append(ftr.GetGeometryRef().Clone())rid.append(ftr.items()["DN"])# Do shrink, maybeifshrink:# Compute shrink factorshrinkFactor=-0.00001*(xMax-xMin)/matrix.shape[1]geoms=[g.Buffer(float(shrinkFactor))forgingeoms]# Do flatten, maybeifflat:geoms=np.array(geoms)rid=np.array(rid)finalGeoms=[]finalRID=[]for_ridinset(rid):smallGeomSet=geoms[rid==_rid]finalGeoms.append(flatten(smallGeomSet)iflen(smallGeomSet)>1elsesmallGeomSet[0])finalRID.append(_rid)else:finalGeoms=geomsfinalRID=rid# Cleanupvector_layer=NonevecDS=NonemaskBand=NonerasBand=Noneraster=None# Done!if_raw:returnfinalGeoms,finalRIDelse:returnpd.DataFrame(dict(geom=finalGeoms,value=finalRID))
If True, shrink all geoms by a tiny amount in order to avoid geometry
overlapping issues
* The total amount shrunk should be very very small
* Generally this should be left as True unless it is ABSOLUTELY
necessary to maintain the same area
Returns:
pandas.DataFrame -> With columns:
–
'geom' -> The contiguous-valued geometries
'value' -> The value for each geometry
defpolygonizeRaster(source,srs=None,flat=False,shrink=True):"""Polygonize a raster or an integer-valued data matrix. Parameters ---------- source : Anything acceptable by loadRaster() The raster datasource to polygonize * The Datatype MUST be of boolean of integer type srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the polygons to create * If not given, the raster's internal srs is assumed flat : bool If True, flattens the resulting geometries which share a contiguous value into a single geometry object shrink : bool If True, shrink all geoms by a tiny amount in order to avoid geometry overlapping issues * The total amount shrunk should be very very small * Generally this should be left as True unless it is ABSOLUTELY necessary to maintain the same area Returns ------- pandas.DataFrame -> With columns: 'geom' -> The contiguous-valued geometries 'value' -> The value for each geometry """# Load the information we will needsource=loadRaster(source)band=source.GetRasterBand(1)maskBand=band.GetMaskBand()ifsrsisNone:srs=SRS.loadSRS(source.GetProjectionRef())# Do polygonizevecDS=gdal.GetDriverByName("Memory").Create("",0,0,0,gdal.GDT_Unknown)vecLyr=vecDS.CreateLayer("mem",srs=srs)# vecDS = gdal.GetDriverByName("ESRI Shapefile").Create("deleteme.tif", 0, 0, 0, gdal.GDT_Unknown )# vecLyr = vecDS.CreateLayer("layer",srs=srs)vecField=ogr.FieldDefn("DN",ogr.OFTInteger)vecLyr.CreateField(vecField)# Polygonize geometryresult=gdal.Polygonize(band,maskBand,vecLyr,0)ifresult!=0:raiseGeoKitGeomError("Failed to polygonize geometry")# Check the geomsftrN=vecLyr.GetFeatureCount()ifftrN==0:# raise GlaesError("No features in created in temporary layer")msg="No features in created in temporary layer"warnings.warn(msg,UserWarning)return# Extract geometries and valuesgeoms=[]rid=[]foriinrange(ftrN):ftr=vecLyr.GetFeature(i)geoms.append(ftr.GetGeometryRef().Clone())rid.append(ftr.items()["DN"])# Do shrink, maybeifshrink:# Compute shrink factorshrinkFactor=-0.00001geoms=[g.Buffer(float(shrinkFactor))forgingeoms]# Do flatten, maybeifflat:geoms=np.array(geoms)rid=np.array(rid)finalGeoms=[]finalRID=[]for_ridinset(rid):smallGeomSet=geoms[rid==_rid]finalGeoms.append(GEOM.flatten(smallGeomSet)iflen(smallGeomSet)>1elsesmallGeomSet[0])finalRID.append(_rid)else:finalGeoms=geomsfinalRID=rid# CleanupvecLyr=NonevecDS=NonemaskBand=None# Done!returnpd.DataFrame(dict(geom=finalGeoms,value=finalRID))
Returns the raster cell number for one or multiple points defined by geometry or lon/lat. Cell numeration
starting in the top left corner cell of the raster with (0,0). Cells with (-1,-1) are out of bounds.
Args:
points (osgeo.ogr.Geometry, tuple, iterable): Can be an osgeo.ogr.Geometry point or an iterable thereof, else a (lon, lat) tuple (in EPSG:4326) or an iterable thereof.
source (gdal.Dataset, str, optional): A gdal.Dataset type raster or a str formatted path to a raster file. Defaults to None.
bounds (tuple, optional): Raster boundaries in EPSG:4326 in the form of (minX, minY, maxX, maxY). Defaults to None.
cellWidth (int, float, optional): The cell width in EPSG:4326 units. Defaults to None.
cellHeight (int, float, optional): The cell height in EPSG:4326 units. Defaults to None.
NOTE: If source is given, all of the others must be None, else they must be given.
Returns:
tuple or iterable: tuple with (X, Y) cell No or an iterable thereof if multiple points were given.
–
defrasterCellNo(points,source=None,bounds=None,cellWidth=None,cellHeight=None):""" Returns the raster cell number for one or multiple points defined by geometry or lon/lat. Cell numeration starting in the top left corner cell of the raster with (0,0). Cells with (-1,-1) are out of bounds. Args: points (osgeo.ogr.Geometry, tuple, iterable): Can be an osgeo.ogr.Geometry point or an iterable thereof, else a (lon, lat) tuple (in EPSG:4326) or an iterable thereof. source (gdal.Dataset, str, optional): A gdal.Dataset type raster or a str formatted path to a raster file. Defaults to None. bounds (tuple, optional): Raster boundaries in EPSG:4326 in the form of (minX, minY, maxX, maxY). Defaults to None. cellWidth (int, float, optional): The cell width in EPSG:4326 units. Defaults to None. cellHeight (int, float, optional): The cell height in EPSG:4326 units. Defaults to None. NOTE: If source is given, all of the others must be None, else they must be given. Returns ------- tuple or iterable: tuple with (X, Y) cell No or an iterable thereof if multiple points were given. """# check and preprocess points inputsifisinstance(points,ogr.Geometry):points=[points]elifisinstance(points,tuple)andlen(points)==2:points=[points]asserthasattr(points,"__iter__"),("points must be an osgeo.ogr.Geometry POINT object, a tuple of (lon, lat) or an iterable of any of them.")ifisinstance(points[0],ogr.Geometry):ifnotall([p.GetGeometryName()=="POINT"forpinpoints]):raiseTypeError("Only POINT geometries allowed")ifnotall([p.GetSpatialReference().IsSame(SRS.loadSRS(4326))forpinpoints]):raiseValueError("SRS of all points must be EPSG:4326")# convert to lons and latspoints=[(p.GetX(),p.GetY())forpinpoints]else:assertall([isinstance(x,tuple)andlen(x)==2andall([isinstance(_x,(int,float))for_xinx])forxinpoints]),("All entries in points must be (lon, lat) tuples with length of 2 and int or float coordinates if given as tuples or iterable thereof.")# get bounds, cellWidth and cellHeight from the inputsifsourceisnotNone:ifnot(boundsisNoneandcellWidthisNoneandcellHeightisNone):raiseValueError("bounds, cellWidth and cellHeight must be None when source raster is given.")ifisinstance(source,str)andnotos.path.isfile(source):raiseFileNotFoundError("source must be an existing raster file if given as string.")try:sourceRasterInfo=rasterInfo(source)bounds=sourceRasterInfo.boundscellWidth=sourceRasterInfo.pixelWidthcellHeight=sourceRasterInfo.pixelHeightrasterWidth=sourceRasterInfo.xWinSizerasterHeight=sourceRasterInfo.yWinSizeexcept:raiseTypeError("source must be path to a gdal.Dataset raster or a gdal.Dataset raster itself if not None.")ifnotsourceRasterInfo.srs.IsSame(SRS.loadSRS(4326)):raiseValueError("raster source must be in epsg:4326")else:ifboundsisNoneorcellWidthisNoneorcellHeightisNone:raiseValueError("If source is None, bounds, cellWidth and cellHeight must be given.")assert(isinstance(bounds,tuple)andlen(bounds)==4andall([isinstance(x,(int,float))forxinbounds])andbounds[0]<bounds[2]andbounds[1]<bounds[3]),"bounds must be a tuple of length = 4 with int or float entries like such: (minX, minY, maxX, maxY)"assertisinstance(cellHeight,(int,float)),"cellHeight must be an int or float."assertisinstance(cellWidth,(int,float)),"cellWidth must be an int or float."# calculate the raster width and height in Nos of cellsrasterWidth=(bounds[2]-bounds[0])/cellWidthrasterHeight=(bounds[3]-bounds[1])/cellHeightassertnp.isclose(rasterWidth%cellWidth,0)ornp.isclose(rasterWidth%cellWidth,cellWidth),(f"rasterWidth {rasterWidth} is not a multiple of cellWidth {cellWidth}")assertnp.isclose(rasterHeight%cellHeight,0)ornp.isclose(rasterHeight%cellHeight,cellHeight),(f"rasterHeight {rasterHeight} is not a multiple of cellHeight {cellHeight}")# calculate the cell NoscellNos=list()outOfBounds=Falseforpointinpoints:X=int(np.floor((point[0]-bounds[0])/cellWidth))Y=int(np.floor((bounds[3]-point[1])/cellHeight))# if any dimension out of boundsif(X<0orX>=rasterWidth)or(Y<0orY>=rasterHeight):X=Y=-1outOfBounds=True# append to cell No listcellNos.append((X,Y))ifoutOfBounds:warnings.warn("points contain out-of-bounds locations!")# return cellNos as tuple for single point, else as list of tuplesiflen(cellNos)==1:returncellNos[0]else:returncellNos
If True, the maximum and minimum value of the raster data are computed.
This is computationally expensive for large rasters and repeated calls of
this function on the same raster should be avoided.
Returns:
namedtuple -> ( srs: The spatial reference system (as an OGR object)
–
dtype: The datatype
flipY: A flag which indicates that the raster starts at the
'bottom' as opposed to at the 'top'
bounds: The (xMin, yMin, xMax, and yMax) values as a tuple
xMin: The minimal X boundary
yMin:The minimal Y boundary
xMax:The maximal X boundary
yMax: The maximal Y boundary
pixelWidth: The raster's pixelWidth
pixelHeight: The raster's pixelHeight
dx:The raster's pixelWidth
dy: The raster's pixelHeight
noData: The noData value used by the raster
scale: The scale value used by the raster
offset: The offset value used by the raster
xWinSize: The width of the raster is pixels
yWinSize: The height of the raster is pixels
meta: The raster's meta data )
defrasterInfo(sourceDS:load_raster_input,compute_statistics:bool=False)->RasterInfo:"""Returns a named tuple containing information relating to the input raster. Parameters ---------- sourceDS : Anything acceptable by loadRaster() The raster datasource compute_statistics : bool; optional If True, the maximum and minimum value of the raster data are computed. This is computationally expensive for large rasters and repeated calls of this function on the same raster should be avoided. Returns ------- namedtuple -> ( srs: The spatial reference system (as an OGR object) dtype: The datatype flipY: A flag which indicates that the raster starts at the 'bottom' as opposed to at the 'top' bounds: The (xMin, yMin, xMax, and yMax) values as a tuple xMin: The minimal X boundary yMin:The minimal Y boundary xMax:The maximal X boundary yMax: The maximal Y boundary pixelWidth: The raster's pixelWidth pixelHeight: The raster's pixelHeight dx:The raster's pixelWidth dy: The raster's pixelHeight noData: The noData value used by the raster scale: The scale value used by the raster offset: The offset value used by the raster xWinSize: The width of the raster is pixels yWinSize: The height of the raster is pixels meta: The raster's meta data ) """output={}sourceDS=loadRaster(sourceDS)# get srsifsourceDS.GetProjectionRef()=="":# return None directly if raster has no srssrs=Noneelse:# generate an srs object if we have srs informationsrs=SRS.loadSRS(sourceDS.GetProjectionRef())output["srs"]=srs# get extent and resolutionsourceBand=sourceDS.GetRasterBand(1)# Is required to get the maximum and minimum value of the raster dataoutput["dtype"]=sourceBand.DataTypeoutput["noData"]=sourceBand.GetNoDataValue()output["scale"]=sourceBand.GetScale()output["offset"]=sourceBand.GetOffset()ifcompute_statisticsisTrue:try:sourceBand.ComputeStatistics(0)maximum_value=sourceBand.GetMaximum()minimum_value=sourceBand.GetMinimum()output["maximum_value"]=maximum_valueoutput["minimum_value"]=minimum_valueexcept:output["maximum_value"]=output["noData"]output["minimum_value"]=output["noData"]else:output["maximum_value"]=Noneoutput["minimum_value"]=NonexSize=sourceBand.XSizeySize=sourceBand.YSizexOrigin,dx,_,yOrigin,_,dy=sourceDS.GetGeoTransform()xMin=xOriginxMax=xOrigin+dx*xSizeifdy<0:yMax=yOriginyMin=yMax+dy*ySizedy=-1*dyoutput["flipY"]=Trueoutput["yAtTop"]=Trueelse:yMin=yOriginyMax=yOrigin+dy*ySizeoutput["flipY"]=Falseoutput["yAtTop"]=Falseoutput["pixelWidth"]=dxoutput["pixelHeight"]=dyoutput["dx"]=dxoutput["dy"]=dyoutput["xMin"]=xMinoutput["xMax"]=xMaxoutput["yMin"]=yMinoutput["yMax"]=yMaxoutput["xWinSize"]=xSizeoutput["yWinSize"]=ySizeoutput["bounds"]=(xMin,yMin,xMax,yMax)output["meta"]=sourceDS.GetMetadata_Dict()output["source"]=sourceDS.GetDescription()output["data_type_name_str"]=gdal.GetDataTypeName(output["dtype"])# clean updelsourceBand,sourceDSraster_info=RasterInfo(**output)returnraster_info
defrasterStats(source,cutline=None,ignoreValue=None,**kwargs):"""Compute basic statistics of the values contained in a raster dataset. Parameters ---------- source : Anything acceptable by loadRaster() The raster datasource cutline : ogr.Geometry; optional The geometry over which to cut out the raster's data * Must be a Polygon or MultiPolygon ignoreValue : numeric A value to ignore when computing the statistics * If the raster source has a 'no Data' value, it is automatically ignored **kwargs * All kwargs are passed on to warp() when 'geom' is given * See gdal.WarpOptions for more details * For example, 'allTouched' may be useful Returns ------- Results from a call to scipy.stats.describe """fromscipy.statsimportdescribesource=loadRaster(source)# Get the matrix to calculate overifcutlineisnotNone:source=warp(source,cutline=cutline,noData=ignoreValue,**kwargs)rawData=extractMatrix(source)dataInfo=rasterInfo(source)# exclude nodata and ignore valuessel=np.ones(rawData.shape,dtype="bool")ifignoreValueisnotNone:np.logical_and(rawData!=ignoreValue,sel,sel)ifdataInfo.noDataisnotNone:np.logical_and(rawData!=dataInfo.noData,sel,sel)# compute statisticsdata=rawData[sel].flatten()returndescribe(data)
Rasterize a vector datasource onto a raster context.
Note:
When creating an 'in memory' raster vs one which is saved to disk, a slightly
different algorithm is used which can sometimes add an extra row of pixels. Be
aware of this if you intend to compare value-matricies directly from rasters
generated with this function.
The values to burn into the raster
* If a numeric is given, all pixels are burned with the specified value
* If a string is given, then one the feature attribute names is expected
A path to an output file
* If output is None, the raster will be created in memory and a dataset
handle will be returned
* If output is given, the raster will be written to disk and nothing will
be returned
The datatype of the represented by the created raster's band
* Options are: Byte, Int16, Int32, Int64, Float32, Float64
* If dtype is None and data is None, the assumed datatype is a 'Byte'
* If dtype is None and data is not None, the datatype will be inferred
from the given data
A flag instructing the output raster to use a compression algorithm
* only useful if 'output' has been defined
* "DEFLATE" used for Linux/Mac, "LZW" used for Windows
Specifies which value should be considered as 'no data' in the created
raster
* Must be the same datatype as the 'dtype' input (or that which is derived)
defrasterize(source:str|pathlib.Path|ogr.Geometry|gdal.Dataset,pixelWidth:numeric,pixelHeight:numeric,srs:srs_input|None=None,bounds:tuple[numeric,numeric,numeric,numeric]|Extent|None=None,where:str|None=None,value:numeric|str=1,output:str|None=None,dtype:geokit_c_data_types_literal|None=None,compress=True,noData=None,overwrite:bool=True,**kwargs,)->gdal.Dataset|str:"""Rasterize a vector datasource onto a raster context. Note: ----- When creating an 'in memory' raster vs one which is saved to disk, a slightly different algorithm is used which can sometimes add an extra row of pixels. Be aware of this if you intend to compare value-matricies directly from rasters generated with this function. Parameters ---------- source : str or ogr.Geometry If str, the path to the vector file to load If ogr.Geometry, an Polygon geometry - Will be immediately turned into a vector pixelWidth : numeric The pixel width of the raster in the working srs * Is 'srs' is not given, these are the units of the source's inherent srs pixelHeight : numeric The pixel height of the raster in the working srs * Is 'srs' is not given, these are the units of the source's inherent srs srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the point to create * If 'bounds' is an Extent object, the bounds' internal srs will override this input bounds : (xMin, yMix, xMax, yMax) or Extent; optional The geographic extents spanned by the raster * If not given, the whole bounds spanned by the input is used where : str; optional An SQL-like where statement to use to filter the vector before rasterizing value : numeric, str The values to burn into the raster * If a numeric is given, all pixels are burned with the specified value * If a string is given, then one the feature attribute names is expected output : str; optional A path to an output file * If output is None, the raster will be created in memory and a dataset handle will be returned * If output is given, the raster will be written to disk and nothing will be returned dtype : str; optional The datatype of the represented by the created raster's band * Options are: Byte, Int16, Int32, Int64, Float32, Float64 * If dtype is None and data is None, the assumed datatype is a 'Byte' * If dtype is None and data is not None, the datatype will be inferred from the given data compress : bool A flag instructing the output raster to use a compression algorithm * only useful if 'output' has been defined * "DEFLATE" used for Linux/Mac, "LZW" used for Windows noData : numeric; optional Specifies which value should be considered as 'no data' in the created raster * Must be the same datatype as the 'dtype' input (or that which is derived) overwrite : bool A flag to overwrite a pre-existing output file * If set to False and an 'output' is specified which already exists, an error will be raised Returns ------- * If 'output' is None: gdal.Dataset * If 'output' is a string: The path to the output is returned (for easy opening) """# Normalize some inputsifisinstance(source,ogr.Geometry):source=createVector(source)else:source=loadVector(source)# Get the vector's infovecinfo=vectorInfo(source)ifsrsisNone:srs=vecinfo.srssrsOkay=Trueelse:srs=SRS.loadSRS(srs)ifsrs.IsSame(vecinfo.srs):srsOkay=Trueelse:srsOkay=False# Look for bounds inputifboundsisNone:bounds=vecinfo.boundsifnotsrsOkay:bounds=GEOM.boundsToBounds(bounds,vecinfo.srs,srs)else:try:bounds=bounds.xyXY# Get a tuple from an Extent objectexcept:pass# Bounds should already be a tuplebounds=UTIL.fitBoundsTo(bounds,pixelWidth,pixelHeight)# Collect rasterization optionsifoutputisNoneandnot"bands"inkwargs:kwargs["bands"]=[1]list_of_data_types=[]ifisinstance(dtype,str):list_of_data_types.append(dtype)list_of_numbers=[]ifisinstance(value,str):kwargs["attribute"]=valuedata_type_of_field_as_string=vecinfo.attribute_data_types_strlist_of_data_types.append(data_type_of_field_as_string[value])else:kwargs["burnValues"]=[value,]list_of_numbers.append(value)ifisinstance(noData,(numeric,bool)):list_of_numbers.append(noData)# minimum_data_type = dtype# just to raise error early if invalid# Do 'in memory' rasterization# We need to follow this path in both cases since the below fails when simultaneously rasterizing and writing to disk (I couldn't figure out why...)ifoutputisNoneornotsrsOkay:minimum_data_type_string=MinimumCDataTypeHandler.get_valid_gdal_data_type_as_string(list_of_numbers=list_of_numbers,minimum_gdal_type_list=list_of_data_types)# Create temporary output fileoutputDS=UTIL.quickRaster(bounds=bounds,srs=srs,dx=pixelWidth,dy=pixelHeight,dtype=minimum_data_type_string,noData=noData,)# Do rasterizetmp=gdal.Rasterize(outputDS,source,where=where,**kwargs)iftmp==0:raiseGeoKitRasterError("Rasterization failed!")outputDS.FlushCache()ifoutputisNone:returnoutputDSelse:ri=RASTER.rasterInfo(outputDS)RASTER.createRasterLike(ri,output=output,data=RASTER.extractMatrix(outputDS))returnoutput# Do a rasterization to a file on diskelse:# Check for existing fileifos.path.isfile(output):ifoverwrite==True:os.remove(output)ifos.path.isfile(output+".aux.xml"):# Because QGIS....os.remove(output+".aux.xml")else:raiseGeoKitRasterError("Output file already exists: %s"%output)# Arrange some inputsaligned=kwargs.pop("targetAlignedPixels",True)ifnot"creationOptions"inkwargs:ifcompress:creation_options=RASTER.COMPRESSION_OPTIONelse:creation_options=[]else:creation_options=kwargs.pop("creationOptions")# Fix the bounds issue by making them just a little bit smaller, which should be fixed by gdalwarpbounds=(bounds[0]+0.001*pixelWidth,bounds[1]+0.001*pixelHeight,bounds[2]-0.001*pixelWidth,bounds[3]-0.001*pixelHeight,)minimum_data_type_constant=MinimumCDataTypeHandler.get_valid_gdal_data_type_as_constant(list_of_numbers=list_of_numbers,minimum_gdal_type_list=list_of_data_types)# Do rasterizetmp=gdal.Rasterize(output,source,outputBounds=bounds,xRes=pixelWidth,yRes=pixelHeight,outputSRS=srs,noData=noData,where=where,creationOptions=creation_options,targetAlignedPixels=aligned,outputType=minimum_data_type_constant,**kwargs,)ifnotUTIL.isRaster(tmp):raiseGeoKitRasterError("Rasterization failed!")returnoutput
defsaveRasterAsTif(source:gdal.Dataset,output:str,**kwargs):"""Write a osgeo.gdal.Dataset in memory to a GeoTiff file to disk. Parameters ---------- source : osgeo.gdal.Dataset output : str A path to an output file Returns ------- str Path to the saved file on disk. """# assert os.path.isdir(os.path.dirname(output)), 'Output folder does not exist!'assertoutput.split(".")[-1]in["tif","tiff",],"Wrong type specified, use *.tif or *.tiff"sourceInfo=rasterInfo(source)data=extractMatrix(source)returncreateRaster(bounds=sourceInfo.bounds,pixelWidth=sourceInfo.dx,pixelHeight=sourceInfo.dy,noData=sourceInfo.noData,dtype=sourceInfo.data_type_name_str,srs=sourceInfo.srs,data=data,output=output,**kwargs,)
Scale a 2-dimensional matrix. For example, a 2x2 matrix, with a scale of 2,
will become a 4x4 matrix. Or scaling a 24x24 matrix with a scale of -3 will
produce an 8x8 matrix.
Scaling UP (positive) results in a dimensionally larger matrix where each
value is repeated scale^2 times
scaling DOWN (negative) results in a dimensionally smaller matrix where each
value is the average of the associated 'up-scaled' block
defscaleMatrix(mat,scale,strict=True):"""Scale a 2-dimensional matrix. For example, a 2x2 matrix, with a scale of 2, will become a 4x4 matrix. Or scaling a 24x24 matrix with a scale of -3 will produce an 8x8 matrix. * Scaling UP (positive) results in a dimensionally larger matrix where each value is repeated scale^2 times * scaling DOWN (negative) results in a dimensionally smaller matrix where each value is the average of the associated 'up-scaled' block Parameters ---------- mat : numpy.ndarray or [[numeric,],] The data to be scaled * Must be two dimensional scale : int or (int, int) The dimensional scaling factor for either both, or independently for the Y and X dimensions * If scaling down, the scaling factors must be a factor of the their associated dimension in the input matrix (unless 'strict' is set to False) strict : bool Whether or not to force a fail when scaling-down by a scaling factor which is not a dimensional factor * Really intended for internal use... * When scaling down by a non-dimensional factor, the matrix will be padded with zeros such that the new matrix has dimensional sizes which are divisible by the scaling factor. The points which are not at the right or bottom boundary are averaged, same as before. The points which lie on the edge however, are also averaged across all the values which lie in those pixels, but they are corrected so that the averaging does NOT take into account the padded zeros. Returns ------- numpy.ndarray Examples -------- INPUT Scaling Factor Output -------------------------------------- | 1 2 | 2 | 1 1 2 2 | | 3 4 | | 1 1 2 2 | | 3 3 4 4 | | 3 3 4 4 | | 1 1 1 1 | -2 | 1.5 2.0 | | 2 2 3 3 | | 5.25 6.75| | 4 4 5 5 | | 6 7 8 9 | | 1 1 1 1 | -3 | 2.55 3.0 | | 2 2 3 3 | * strict=False | 7.0 9 | | 4 4 5 5 | | 6 7 8 9 | *padded* -------------------------- | 1 1 1 1 0 0 | | 2 2 3 3 0 0 | | 4 4 5 5 0 0 | | 6 7 8 9 0 0 | | 0 0 0 0 0 0 | | 0 0 0 0 0 0 | """# unpack scaletry:yScale,xScale=scaleexcept:yScale,xScale=scale,scale# check for intsifnot(isinstance(xScale,int)andisinstance(yScale,int)):raiseValueError("scale must be integer types")ifxScale==0andyScale==0:returnmat# no scaling (it would just be silly to call this)elifxScale>0andyScale>0:# scale upout=np.zeros((mat.shape[0]*yScale,mat.shape[1]*xScale),dtype=mat.dtype)foryoinrange(yScale):forxoinrange(xScale):out[yo::yScale,xo::xScale]=matelifxScale<0andyScale<0:# scale downxScale=-1*xScaleyScale=-1*yScale# ensure scale is a factor of both xSize and ySizeifstrict:ifnot(mat.shape[0]%yScale==0andmat.shape[1]%xScale==0):raiseGeoKitError("Matrix can only be scaled down by a factor of it's dimensions")yPad=0xPad=0else:# get the amount to pad in the y directionyPad=yScale-mat.shape[0]%yScale# get the amount to pad in the x directionxPad=xScale-mat.shape[1]%xScaleifyPad==yScale:yPad=0ifxPad==xScale:xPad=0# Do y-paddingifyPad>0:mat=np.concatenate((mat,np.zeros((yPad,mat.shape[1]))),0)ifxPad>0:mat=np.concatenate((mat,np.zeros((mat.shape[0],xPad))),1)out=np.zeros((mat.shape[0]//yScale,mat.shape[1]//xScale),dtype="float")foryoinrange(yScale):forxoinrange(xScale):out+=mat[yo::yScale,xo::xScale]out=out/(xScale*yScale)# Correct the edges if a padding was providedifyPad>0:# fix the right edge EXCLUDING the bot-left pointout[:-1,-1]*=yScale/(yScale-yPad)ifxPad>0:# fix the bottom edge EXCLUDING the bot-left pointout[-1,:-1]*=xScale/(xScale-xPad)ifyPad>0:out[-1,-1]*=yScale*xScale/(yScale-yPad)/(xScale-xPad)# fix the bot-left pointelse:# we have both a scaleup and a scale downraiseGeoKitError("Dimensions must be scaled in the same direction")returnout
Shift a polygon in longitudinal and/or latitudinal direction.
Inputs:
geom : The geometry to be shifted
- a single ogr Geometry object of POINT, LINESTRING, POLYGON or MULTIPOLYGON type
- NOTE: Accepts only 2D geometries, z value must be zero.
lonShift - (int, float) : The shift in longitudinal direction in units of the geom srs, may be positive or negative
latShift - (int, float) : The shift in latitudinal direction in units of the geom srs, may be positive or negative
Returns :
osgeo.ogr.Geometry object of the input type with shifted coordinates
defshift(geom,lonShift=0,latShift=0):"""Shift a polygon in longitudinal and/or latitudinal direction. Inputs: geom : The geometry to be shifted - a single ogr Geometry object of POINT, LINESTRING, POLYGON or MULTIPOLYGON type - NOTE: Accepts only 2D geometries, z value must be zero. lonShift - (int, float) : The shift in longitudinal direction in units of the geom srs, may be positive or negative latShift - (int, float) : The shift in latitudinal direction in units of the geom srs, may be positive or negative Returns : --------- osgeo.ogr.Geometry object of the input type with shifted coordinates """ifnotisinstance(geom,ogr.Geometry):raiseTypeError(f"geom must be of type osgeo.ogr.Geometry")geom=geom.Clone()# do not modify input geom# first get srs of input geom_srs=geom.GetSpatialReference()# get the dimension of the geometrydims=geom.GetCoordinateDimension()assertdimsin[2,3],f"Only 2D and 3D points are supported, but got {dims}D"# define sub method to shift collection of single pointsdef_movePoints(pointCollection,lonShift,latShift):"""Auxiliary function shifting individual points."""points=list()foriinrange(len(str(pointCollection).split(","))):points.append(pointCollection.GetPoint(i))# shift the points individuallypoints_shifted=list()forpinpoints:assertp[2]==0,f"All z-values must be zero!"points_shifted.append((p[0]+lonShift,p[1]+latShift))returnpoints_shifted# first check if geom is a point and shiftif"POINT"ingeom.GetGeometryName():p=geom.GetPoint()point_shifted=point((p[0]+lonShift,p[1]+latShift),srs=_srs)ifdims==2:point_shifted.FlattenTo2D()returnpoint_shifted# else check if line and adaptelif"LINESTRING"ingeom.GetGeometryName()andnot"MULTILINE"ingeom.GetGeometryName():assertnotgeom.IsEmpty(),f"Line is empty"line_shifted=line(_movePoints(pointCollection=geom,lonShift=lonShift,latShift=latShift),srs=_srs,)ifdims==2:line_shifted.FlattenTo2D()returnline_shifted# else check if we have a (multi)polygonelif"POLYGON"ingeom.GetGeometryName():assertnotgeom.IsEmpty(),f"Polygon is empty"ifnot"MULTIPOLYGON"ingeom.GetGeometryName():# only a simple polygon, generate single entry list to allow iterationgeom=[geom]# iterate over individual polygonsforip,polyinenumerate(geom):assert"POLYGON"inpoly.GetGeometryName(),f"MULTIPOLYGON is not composed of only POLYGONS"# iterate over sub linear ringsforir,ringinenumerate(poly):assert"LINEARRING"inring.GetGeometryName(),(f"POLYGON (or sub polygon of MULTIPOLYGON) is not composed of only LINEARRINGS")poly_shifted=polygon(_movePoints(pointCollection=ring,lonShift=lonShift,latShift=latShift),srs=_srs,)ifip==0andir==0:multi_shifted=poly_shiftedelse:multi_shifted=multi_shifted.Union(poly_shifted)# the shifted polygon should have the same dimensions as the inputifdims==2:multi_shifted.FlattenTo2D()returnmulti_shiftedelse:raiseTypeError(f"geom must be a 'POINT', 'LINESTRING', 'POLYGON' or 'MULTIPOLYGON' osgeo.ogr.Geometry")
Removes raster polygons smaller than a provided threshold size (in pixels) and
replaces them with the pixel value of the largest neighbour polygon.
It is useful if you have a large amount of small areas on your raster map.
Returns:
* If 'output' is None : gdal.Dataset
–
* If 'output' is a string : The path to the output is returned (for easy opening)
–
defsieve(source,threshold:int=100,connectedness:Literal[4,8]=8,mask="none",quiet_flag:bool=False,output:str|None=None,**kwargs,)->gdal.Dataset|str:""" Removes raster polygons smaller than a provided threshold size (in pixels) and replaces them with the pixel value of the largest neighbour polygon. It is useful if you have a large amount of small areas on your raster map. Returns ------- * If 'output' is None : gdal.Dataset * If 'output' is a string : The path to the output is returned (for easy opening) Parameters ---------- source : Anything acceptable by loadRaster() threshold (int): minimum polygon size (number of pixels) to retain. connectedness (int): either 4 indicating that diagonal pixels are not considered directly adjacent for polygon membership purposes or 8 indicating they are. mask (str): 'none' or 'default'. An optional mask band. All pixels in the mask band with a value other than zero will be considered suitable for inclusion in polygons. quiet_flag (bool): 0 or 1. Callback for reporting algorithm progress output : str; optional The path on disk where the new raster should be created **kwargs: * All kwargs are passed on to SieveFilter() * See gdal.SieveFilter for more details """gdal.AllRegister()try:gdal.SieveFilterexcept:print("")print('gdal.SieveFilter() not available. You are likely using "old gen"')print("bindings or an older version of the next gen bindings.")print("")### Open source file# if output is None:# src_ds = gdal.Open( source, gdal.GA_Update )# else:# src_ds = gdal.Open( source, gdal.GA_ReadOnly )src_ds=loadRaster(source)ifsrc_dsisNone:print("Unable to open %s "%source)srcband=src_ds.GetRasterBand(1)ifmask=="default":maskband=srcband.GetMaskBand()elifmask=="none":maskband=Noneelse:mask_ds=loadRaster(mask)maskband=mask_ds.GetRasterBand(1).GetMaskBand()### Create output file if one is specified.ifoutputisnotNone:format="GTiff"driver=gdal.GetDriverByName(format)dst_ds=driver.Create(output,src_ds.RasterXSize,src_ds.RasterYSize,1,srcband.DataType,COMPRESSION_OPTION,)# ["COMPRESS=LZW"]else:format="Mem"driver=gdal.GetDriverByName(format)# create a raster in memorydst_ds=driver.Create("",src_ds.RasterXSize,src_ds.RasterYSize,1,srcband.DataType,COMPRESSION_OPTION,)# ["COMPRESS=LZW"]wkt=src_ds.GetProjection()ifwkt!="":dst_ds.SetProjection(wkt)dst_ds.SetGeoTransform(src_ds.GetGeoTransform())dstband=dst_ds.GetRasterBand(1)ifquiet_flag:prog_func=Noneelse:prog_func=gdal.TermProgress_nocbresult=gdal.SieveFilter(srcband,maskband,dstband,threshold,connectedness,callback=prog_func,**kwargs,)dst_ds.FlushCache()# Return raster if in memoryifoutputisNone:returndst_ds# Return output path to raster if on diskelse:src_ds=Nonedst_ds=Nonemask_ds=Nonereturnoutput
defsubTiles(geom,zoom,checkIntersect=True,asGeom=False):""" Generate a collection of tiles which encompass the passed geometry. Parameters ---------- geom : ogr.Geometry The geometry to be analyzed zoom : int The zoom level to generate tiles on checkIntersect : bool If True, only tiles which overlap the given geometry are returned asGeom : bool If True, geometry object corresponding to each tile is yielded, instead of (xi,yi,zoom) tuples Returns ------- If asGeom is False: Generates (xi, yi, zoom) tuples If asGeom is True: Generates Geometry objects """geom4326=transform(geom,toSRS=SRS.EPSG4326)ifcheckIntersect:geom3857=transform(geom,toSRS=SRS.EPSG3857)xmin,xmax,ymin,ymax=geom4326.GetEnvelope()tl_xi,tl_yi=smopy.deg2num(ymax,xmin,zoom)br_xi,br_yi=smopy.deg2num(ymin,xmax,zoom)forxiinrange(tl_xi,br_xi+1):foryiinrange(tl_yi,br_yi+1):ifcheckIntersectorasGeom:gtile=tile(xi,yi,zoom)ifcheckIntersect:ifnotgeom3857.Intersects(gtile):continueifasGeom:yieldgtileelse:yieldTile(xi,yi,zoom)
deftile(xi,yi,zoom):"""Generates a box corresponding to a tile used for "slippery maps". Parameters ---------- xi : int The tile's X-index - Range depends on zoom value yi : int The tile's Y-index - Range depends on zoom value zoom : int The tile's zoom index - Range is between 0 and 18 Returns ------- ogr.Geometry """tl=smopy.num2deg(xi-0.0,yi+1.0,zoom)[::-1]br=smopy.num2deg(xi+1.0,yi-0.0,zoom)[::-1]o=SRS.xyTransform([tl,br],fromSRS=SRS.EPSG4326,toSRS=SRS.EPSG3857,outputFormat="xy")returnbox(o.x.min(),o.y.min(),o.x.max(),o.y.max(),srs=SRS.EPSG3857)
deftileAt(x,y,zoom,srs):"""Generates a box corresponding to a tile at the coordinates 'x' and 'y' in the given srs,. Parameters ---------- x : float The X coordinate to search for a tile around y : float The Y coordinate to search for a tile around zoom : int The tile's zoom index - Range is between 0 and 18 srs : anything acceptable to SRS.loadSRS The SRS of the given 'x' & 'y' coordinates Returns ------- ogr.Geometry """t=SRS.tileIndexAt(x=x,y=y,zoom=zoom,srs=srs)returntile(t.xi,t.yi,t.zoom)
deftileIndexAt(x,y,zoom,srs):"""Get the "slippery tile" index at the given zoom, around the coordinates ('x', 'y') within the specified 'srs'. """srs=loadSRS(srs)iterable_input=isinstance(x,Iterable)orisinstance(y,Iterable)ifnotsrs.IsSame(EPSG4326):pt=xyTransform(x,y,fromSRS=srs,toSRS=EPSG4326,outputFormat="xy")ifiterable_input:x,y=pt.x,pt.yelse:x,y=pt.x[0],pt.y[0]ifiterable_input:x=np.array(x)y=np.array(y)xi,yi=smopy.deg2num(y,x,zoom)returnTile(xi,yi,zoom)
deftileize(geom,zoom):"""Deconstruct a given geometry into a set of tiled geometries. Returns: Generator of ogr.Geometry objects """geom=transform(geom,toSRS=SRS.EPSG3857)fortile_insubTiles(geom,zoom,asGeom=True,checkIntersect=True):yieldgeom.Intersection(tile_)
(Anything acceptable to geokit.srs.loadSRS(); optional, default:
None
)
–
The srs of the input geometries
* Only needed if an SRS cannot be inferred from the geometry inputs or
is, for whatever reason, the geometry's SRS is wrong
If True, will revert the 360° shift that is applied by PROJ to points
beyond the antimeridian when transforming to EPSG:4326 to avoid invalid,
distorted geometries when points are only partially shifted. Will then
devliver a geometry that may partially be outside the -180° to 180°
longitude range. By default False.
An optional segmentation length to apply to the input geometries BEFORE
transformation occurs. The input geometries will be segmented such that
no line segment is longer than the given segment size
* Units are in the input geometry's native unit
* Use this for a more detailed transformation!
Returns:
Geometry or [Geometry]
–
Note:
When inferring the SRS from the given geometries, only the FIRST geometry
is checked for an existing SRS
deftransform(geoms,toSRS,fromSRS=None,revert360degProj=False,segment=None)->ogr.Geometry|list[ogr.Geometry]:"""Transform a geometry, or a list of geometries, from one SRS to another. Parameters ---------- geoms : ogr.Geometry or [ogr.Geometry, ] The geometry or geometries to transform * All geometries must have the same spatial reference toSRS : Anything acceptable to geokit.srs.loadSRS() The srs of the output geometries fromSRS : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the input geometries * Only needed if an SRS cannot be inferred from the geometry inputs or is, for whatever reason, the geometry's SRS is wrong revert360degProj : bool; optional If True, will revert the 360° shift that is applied by PROJ to points beyond the antimeridian when transforming to EPSG:4326 to avoid invalid, distorted geometries when points are only partially shifted. Will then devliver a geometry that may partially be outside the -180° to 180° longitude range. By default False. segment : float; optional An optional segmentation length to apply to the input geometries BEFORE transformation occurs. The input geometries will be segmented such that no line segment is longer than the given segment size * Units are in the input geometry's native unit * Use this for a more detailed transformation! Returns ------- ogr.Geometry or [ogr.Geometry, ] Note: ----- When inferring the SRS from the given geometries, only the FIRST geometry is checked for an existing SRS """# make sure geoms is a listifisinstance(geoms,ogr.Geometry):returnSingle=Truegeoms=[geoms,]else:# assume geoms is iterablereturnSingle=Falsetry:geoms=list(geoms)exceptExceptionase:msg="Geoms is neither an osgeo.ogr.Geometry instance nor an iterable that can be converted to a list."warnings.warn(msg,UserWarning)raisee# make sure geoms is a listiffromSRSisNone:fromSRS=geoms[0].GetSpatialReference()assertall([g.GetSpatialReference().IsSame(fromSRS)forgingeoms]),f"geoms have different SRS."iffromSRSisNone:raiseGeoKitGeomError("Could not determine fromSRS from geometry")# load srs'sfromSRS=SRS.loadSRS(fromSRS)toSRS=SRS.loadSRS(toSRS)# make a transformertrx=osr.CoordinateTransformation(fromSRS,toSRS)# Do transformationgeoms=[g.Clone()forgingeoms]ifnotsegmentisNone:[g.Segmentize(segment)forgingeoms]r=[g.Transform(trx)forgingeoms]ifsum(r)>0:# check for errorsraiseGeoKitGeomError("Errors in geometry transformations")ifrevert360degProjandtoSRS.IsSame(SRS.loadSRS(4326)):def_revert_antimeridian_proj(geom,max_width=180):"""Shift back points across the antimeridian to undo the effect of the PROJ function."""def_shift_points(g):"""Shifts points back by +/-360° if PROJ shifted them due to the antimeridian."""prev_x,_,_=g.GetPoint(0)forjinrange(1,g.GetPointCount()):x,y,z=g.GetPoint(j)dx=x-prev_xifdx>max_width:x-=360elifdx<-max_width:x+=360g.SetPoint(j,x,y,z)prev_x=xreturngeom# apply at different levels (directly to lines, indirectly to rings per polygon or recursively to multi-polygons/lines)geom_type=geom.GetGeometryType()ifgeom_typein(ogr.wkbPolygon,ogr.wkbPolygon25D):foriinrange(geom.GetGeometryCount()):_shift_points(geom.GetGeometryRef(i))elifgeom_typein(ogr.wkbMultiPolygon,ogr.wkbMultiPolygon25D,ogr.wkbMultiLineString,ogr.wkbMultiLineString25D,):foriinrange(geom.GetGeometryCount()):_revert_antimeridian_proj(geom.GetGeometryRef(i),max_width)elifgeom_typein(ogr.wkbLineString,ogr.wkbLineString25D):_shift_points(geom)# points are returned unchanged since here the projection across the antimeridian is desiredreturngeom# polygons or lines crossing the antimeridian may have been distorted by PROJ operationgeoms=[_revert_antimeridian_proj(geom=g)forgingeoms]# make sure all geoms are validtry:assertall([g.IsValid()forgingeoms])# no msg, fall back to except only if needed to save timeexcept:geoms=[g.Buffer(0)forgingeoms]# trick to reset boundary of unnecessarily invalid geomsassertall([g.IsValid()forgingeoms]),f"Some geometries are invalid after transformation"# Done!ifreturnSingle:returngeoms[0]else:returngeoms
bounds : The source's boundaries (in the srs's units),
xMin : The source's xMin boundaries (in the srs's units),
yMin : The source's xMax boundaries (in the srs's units),
xMax : The source's yMin boundaries (in the srs's units),
yMax : The source's yMax boundaries (in the srs's units),
count : The number of features in the source,
attributes : The attribute titles for the source's features,)
defvectorInfo(source)->vecInfo:"""Extract general information about a vector source. Determines: Parameters ---------- source : Anything acceptable by loadVector() The vector datasource to read from Returns ------- namedtuple -> (srs : The source's SRS system, bounds : The source's boundaries (in the srs's units), xMin : The source's xMin boundaries (in the srs's units), yMin : The source's xMax boundaries (in the srs's units), xMax : The source's yMin boundaries (in the srs's units), yMax : The source's yMax boundaries (in the srs's units), count : The number of features in the source, attributes : The attribute titles for the source's features,) """info={}vecDS=loadVector(source)vecLyr:ogr.Layer=vecDS.GetLayer()info["srs"]=vecLyr.GetSpatialRef()xMin,xMax,yMin,yMax=vecLyr.GetExtent()info["bounds"]=(xMin,yMin,xMax,yMax)info["xMin"]=xMininfo["xMax"]=xMaxinfo["yMin"]=yMininfo["yMax"]=yMaxinfo["count"]=vecLyr.GetFeatureCount()info["source"]=vecDS.GetDescription()info["attributes"]=[]info["attribute_data_types_constant"]={}info["attribute_data_types_str"]={}layerDef:ogr.FeatureDefn=vecLyr.GetLayerDefn()forlayer_numberinrange(layerDef.GetFieldCount()):field_definition:ogr.FieldDefn=layerDef.GetFieldDefn(layer_number)field_name=field_definition.GetName()field_type_constant=field_definition.GetType()field_type_constant_str=gdal.GetDataTypeName(field_type_constant)info["attributes"].append(field_name)info["attribute_data_types_constant"][field_name]=field_type_constantinfo["attribute_data_types_str"][field_name]=field_type_constant_strreturnvecInfo(**info)
The resampling algorithm to use when translating pixel values
* Knowing which option to use can have significant impacts!
* Options are: near , bilinear, cubic,
cubicspline, lanczos, average, rms, mode,
max, min, med, Q1, Q3, sum
The cutline to limit the drawn data too
* If a string is given, it must be a path to a vector file
* Values outside of the cutline are given the value 'cutlineFillValue'
* Requires a warp
(Type, str, or numpy-dtype; optional, default:
None
)
–
If given, forces the processed data to be a particular datatype
* Only required if this value should be changed
* Example
- A python numeric type such as bool, int, or float
- A Numpy datatype such as numpy.uint8 or numpy.float64
- a String such as "Byte", "UInt16", or "Double"
All keyword arguments are passed on to a call to gdal.WarpOptions
Use these to fine-tune the warping procedure
Key Options are (from gdal.WarpOptions):
format --- output format ("GTiff", etc...)
targetAlignedPixels --- whether to force output bounds to be multiple
of output resolution
workingType --- working type (gdal.GDT_Byte, etc...)
warpMemoryLimit --- size of working buffer in bytes
creationOptions --- list of creation options
srcNodata --- source nodata value(s)
dstNodata --- output nodata value(s)
multithread --- whether to multithread computation and I/O operations
cutlineWhere --- cutline WHERE clause
cropToCutline --- whether to use cutline extent for output bounds
setColorInterpretation --- whether to force color interpretation of
input bands to output bands
Returns:
* If 'output' is None: gdal.Dataset
–
* If 'output' is a string: The path to the output is returned (for easy opening)
–
defwarp(source:load_raster_input,resampleAlg:Literal["near","bilinear","cubic","cubicspline","lanczos","average","rms","mode","max","min","med","Q1","Q3","sum",]="bilinear",cutline:str|ogr.Geometry|None=None,output:str|pathlib.Path|None=None,pixelHeight:numeric|None=None,pixelWidth:numeric|None=None,srs:srs_input|None=None,bounds:tuple[numeric,numeric,numeric,numeric]|None=None,dtype:None|geokit_c_data_types_literal=None,noData:numeric|None=None,overwrite:bool=True,meta:None|dict[str,str]=None,**kwargs,)->gdal.Dataset|str:"""Warps a given raster source to another context. * Can be used to 'warp' a raster in memory to a raster on disk Note: ----- Unless manually altered as keyword arguments, the gdal.Warp options 'targetAlignedPixels' and 'copyMetadata' are both set to True Parameters ---------- source : Anything acceptable by loadRaster() The raster datasource to draw srs : Anything acceptable to geokit.srs.loadSRS(); optional The srs of the resulting raster * If not given, the raster's internal srs is assumed resampleAlg : str; optional The resampling algorithm to use when translating pixel values * Knowing which option to use can have significant impacts! * Options are: near , bilinear, cubic, cubicspline, lanczos, average, rms, mode, max, min, med, Q1, Q3, sum cutline : str or ogr.Geometry; optional The cutline to limit the drawn data too * If a string is given, it must be a path to a vector file * Values outside of the cutline are given the value 'cutlineFillValue' * Requires a warp output : str; pathlib.Path optional The path on disk where the new raster should be created pixelHeight : numeric; optional The pixel height (y-resolution) of the output raster * Only required if this value should be changed pixelWidth : numeric; optional The pixel width (x-resolution) of the output raster * Only required if this value should be changed bounds : tuple; optional The (xMin, yMin, xMax, yMax) limits of the output raster * Only required if this value should be changed dtype : Type, str, or numpy-dtype; optional If given, forces the processed data to be a particular datatype * Only required if this value should be changed * Example - A python numeric type such as bool, int, or float - A Numpy datatype such as numpy.uint8 or numpy.float64 - a String such as "Byte", "UInt16", or "Double" noData : numeric; optional Replaces all previous noData values with this value in the output raster. meta: dict; optional: contains a key value pair that is passed to the output gdal.dataset using the SetMetadataItem method. **kwargs: * All keyword arguments are passed on to a call to gdal.WarpOptions * Use these to fine-tune the warping procedure * Key Options are (from gdal.WarpOptions): format --- output format ("GTiff", etc...) targetAlignedPixels --- whether to force output bounds to be multiple of output resolution workingType --- working type (gdal.GDT_Byte, etc...) warpMemoryLimit --- size of working buffer in bytes creationOptions --- list of creation options srcNodata --- source nodata value(s) dstNodata --- output nodata value(s) multithread --- whether to multithread computation and I/O operations cutlineWhere --- cutline WHERE clause cropToCutline --- whether to use cutline extent for output bounds setColorInterpretation --- whether to force color interpretation of input bands to output bands Returns ------- * If 'output' is None: gdal.Dataset * If 'output' is a string: The path to the output is returned (for easy opening) """# open source and get infosource=loadRaster(source)dsInfo=rasterInfo(sourceDS=source,compute_statistics=True)ifdsInfo.scale!=1.0ordsInfo.offset!=0.0:isAdjusted=Trueelse:isAdjusted=False# Handle potentially missing argumentsifsrsisnotNone:srs=SRS.loadSRS(srs)ifsrsisNone:srs=dsInfo.srssrsOkay=Trueelse:ifsrs.IsSame(dsInfo.srs):srsOkay=Trueelse:srsOkay=FalseifboundsisNone:ifsrsOkay:bounds=dsInfo.boundselse:bounds=GEOM.boundsToBounds(dsInfo.bounds,dsInfo.srs,srs)ifpixelHeightisNone:ifsrsOkay:pixelHeight=dsInfo.dyelse:pixelHeight=(bounds[3]-bounds[1])/(dsInfo.yWinSize*1.1)ifpixelWidthisNone:ifsrsOkay:pixelWidth=dsInfo.dxelse:pixelWidth=(bounds[2]-bounds[0])/(dsInfo.xWinSize*1.1)bounds=UTIL.fitBoundsTo(bounds,pixelWidth,pixelHeight)ifnoDataisNone:noDataRead=dsInfo.noDataelse:noDataRead=noDatalist_of_numbers=[]ifisinstance(noDataRead,(numeric,bool)):list_of_numbers.append(noDataRead)elifnoDataReadisNone:passelse:raiseGeoKitRasterError("noData must be a numeric or boolean value but got: %s"%str(type(noDataRead)))list_of_datatypes=[]ifisinstance(dtype,str):list_of_datatypes.append(dtype)elifdtypeisNone:passelse:raiseGeoKitRasterError("dtype must be a gdal data type, string or None value but got: %s"%str(type(dtype)))list_of_datatypes.append(dsInfo.data_type_name_str)list_of_numbers.append(dsInfo.minimum_value)list_of_numbers.append(dsInfo.maximum_value)# If a cutline is given, create the outputifcutlineisnotNone:ifisinstance(cutline,ogr.Geometry):tempdir=TemporaryDirectory()cutline=UTIL.quickVector(cutline,output=os.path.join(tempdir.name,"tmp.shp"))# cutline is already a path to a vectorelifUTIL.isVector(cutline):tempdir=Noneelse:raiseGeoKitRasterError("cutline must be a Geometry or a path to a shape file")# Workflow depends on whether or not we have an outputifisinstance(output,pathlib.Path):output=str(output)ifisinstance(output,str):# Simply do a translateifos.path.isfile(output):ifoverwriteisTrue:os.remove(output)ifos.path.isfile(output+".aux.xml"):# Because QGIS....os.remove(output+".aux.xml")else:raiseGeoKitRasterError("Output file already exists: %s"%output)gdal_data_type_constant=MinimumCDataTypeHandler.get_valid_gdal_data_type_as_constant(list_of_numbers=list_of_numbers,minimum_gdal_type_list=list_of_datatypes,user_defined_minimum_gdal_type=dtype,)# # Check some for bad input configurations# if not srs is None:# if (pixelHeight is None or pixelWidth is None):# raise GeoKitRasterError("When warping between srs's and writing to a file, pixelWidth and pixelHeight must be given")# Arange inputsco=kwargs.pop("creationOptions",COMPRESSION_OPTION)copyMeta=kwargs.pop("copyMetadata",True)aligned=kwargs.pop("targetAlignedPixels",True)# Fix the bounds issue by making them just a little bit smaller, which should be fixed by gdalwarpbounds=(bounds[0]+0.001*pixelWidth,bounds[1]+0.001*pixelHeight,bounds[2]-0.001*pixelWidth,bounds[3]-0.001*pixelHeight,)# Let gdalwarp do everything...gdal_warp_options=gdal.WarpOptions(outputType=gdal_data_type_constant,xRes=pixelWidth,yRes=pixelHeight,creationOptions=co,outputBounds=bounds,dstSRS=srs,dstNodata=noDataRead,resampleAlg=resampleAlg,copyMetadata=copyMeta,targetAlignedPixels=aligned,cutlineDSName=cutline,**kwargs,)result_dataset=gdal.Warp(destNameOrDestDS=output,srcDSOrSrcDSTab=source,options=gdal_warp_options)ifnotUTIL.isRaster(result_dataset):raiseGeoKitRasterError("Failed to translate raster")destination_raster=outputelse:if"cropToCutline"inkwargs:msg="The 'cropToCutline' option is not taken into account when writing to a raster in memory. Try using geokit.Extent.warp instead"warnings.warn(msg,UserWarning)gdal_data_type_string=MinimumCDataTypeHandler.get_valid_gdal_data_type_as_string(list_of_numbers=list_of_numbers,minimum_gdal_type_list=list_of_datatypes,user_defined_minimum_gdal_type=dtype,)# Warp to a raster in memorydestination_raster=UTIL.quickRaster(bounds=bounds,srs=srs,dx=pixelWidth,dy=pixelHeight,dtype=gdal_data_type_string,noData=noDataRead,)# Do a warpresult_dataset=gdal.Warp(destination_raster,source,resampleAlg=resampleAlg,cutlineDSName=cutline,**kwargs)destination_raster.FlushCache()# Do we have meta data?ifmetaisnotNone:ifisinstance(destination_raster,str):loaded_output_raster=loadRaster(destination_raster,1)else:loaded_output_raster=destination_rasterfork,vinmeta.items():loaded_output_raster.SetMetadataItem(k,v)# FlushCache writes all changes to diskloaded_output_raster.FlushCache()ifcutlineisnotNone:deltempdirreturndestination_raster
Convenience function to warp a raster to the context of another raster
as returned from a call to rasterInfo(contextSource).
dataSource : Anything acceptable by loadRaster()
The raster data source to draw
contextSource : Anything acceptable by loadRaster()
The raster context source to draw, i.e. the data source raster will be
warped to this pixelWidth, pixelHeight, bounds, extent, srs etc.
copyMetadata : bool, optional
If True, the metadata of the dataSource raster will be copied, else
metadata will be empty or as possibly provided in kwargs. Defaults to False.
**kwargs
All kwargs will be passed on to raster.warp()
NOTE: If no 'dtype' value as kwargs is given, dtype will be defined
automatically based on the value range, this can be time-consuming
depending on data size. Avoid by specifying dtype explicitly.
defwarpLike(dataSource:load_raster_input,contextSource:load_raster_input,copyMetadata:bool=False,**kwargs):""" Convenience function to warp a raster to the context of another raster as returned from a call to rasterInfo(contextSource). dataSource : Anything acceptable by loadRaster() The raster data source to draw contextSource : Anything acceptable by loadRaster() The raster context source to draw, i.e. the data source raster will be warped to this pixelWidth, pixelHeight, bounds, extent, srs etc. copyMetadata : bool, optional If True, the metadata of the dataSource raster will be copied, else metadata will be empty or as possibly provided in kwargs. Defaults to False. **kwargs All kwargs will be passed on to raster.warp() NOTE: If no 'dtype' value as kwargs is given, dtype will be defined automatically based on the value range, this can be time-consuming depending on data size. Avoid by specifying dtype explicitly. """ifUTIL.isRaster(dataSource):dataInfo=rasterInfo(dataSource)ifnotisinstance(dataInfo,RasterInfo):raiseGeoKitRasterError("Could not understand dataSource")ifUTIL.isRaster(contextSource):contextInfo=rasterInfo(contextSource)ifnotisinstance(contextInfo,RasterInfo):raiseGeoKitRasterError("Could not understand contextSource")# first get data-related parameters from DATA sourceifcopyMetadata:if"meta"inkwargs:raiseGeoKitRasterError("If metadata is given as 'meta' in kwargs, copyMetadata cannot be True!")meta=dataInfo.metaelse:meta=kwargs.pop("meta",None)# then get context related parameters from CONTEXT sourcebounds=kwargs.pop("bounds",contextInfo.bounds)pixelWidth=kwargs.pop("pixelWidth",contextInfo.pixelWidth)pixelHeight=kwargs.pop("pixelHeight",contextInfo.pixelHeight)srs=kwargs.pop("srs",contextInfo.srs)noData=kwargs.pop("noData",contextInfo.noData)returnwarp(source=dataSource,bounds=bounds,pixelWidth=pixelWidth,pixelHeight=pixelHeight,srs=srs,noData=noData,meta=meta,**kwargs,)
defxyTransform(*args,toSRS:srs_input,fromSRS:srs_input,outputFormat:Literal["raw","xy","xyz"]="raw")->list[tuple]|TransformedPointsXY|TransformedPointsXYZ:"""Transform xy points between coordinate systems. Parameters ---------- xy : A single, or an iterable of (x,y) tuples The coordinates to transform toSRS : Anything acceptable by geokit.srs.loadSRS The srs of the output points fromSRS : Anything acceptable by geokit.srs.loadSRS The srs of the input points outputFormat : str Determine return value format * if 'raw', the raw output from osr.TransformPoints is given * if 'xy', or 'xyz' the points are given as named tuples Returns ------- list of tuples, or namedtuple * See the point for the 'outputFormat' argument """# load srs'sfromSRS=loadSRS(fromSRS)toSRS=loadSRS(toSRS)# make a transformertrx=osr.CoordinateTransformation(fromSRS,toSRS)# Do transformationiflen(args)==0:raiseGeoKitSRSError("no positional inputs given")eliflen(args)==1:xy=args[0]ifisinstance(xy,tuple):x,y=xyout=[trx.TransformPoint(x,y),]else:out=trx.TransformPoints(xy)eliflen(args)==2:x=np.array(args[0])y=np.array(args[1])xy=np.column_stack([x,y])out=trx.TransformPoints(xy)else:raiseGeoKitSRSError("Too many positional inputs")# Done!ifoutputFormat=="raw":returnoutelifoutputFormat=="xy":x=np.array([o[0]foroinout])y=np.array([o[1]foroinout])returnTransformedPointsXY(x,y)elifoutputFormat=="xyz":x=out[:,0]y=out[:,1]z=out[:,2]returnTransformedPointsXYZ(x,y,z)