Python osgeo.ogr.Geometry() Examples

The following are code examples for showing how to use osgeo.ogr.Geometry(). They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don't like. You can also save this page to your account.

Example 1
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 8 votes vote down vote up
def exporttogeojson(geojsonfilename, buildinglist):
    #
    # geojsonname should end with .geojson
    # building list should be list of dictionaries
    # list of Dictionaries {'ImageId': image_id, 'BuildingId': building_id, 'polyPix': poly,
    #                       'polyGeo': poly}
    # image_id is a string,
    # BuildingId is an integer,
    # poly is a ogr.Geometry Polygon
    #
    # returns geojsonfilename

    driver = ogr.GetDriverByName('geojson')
    if os.path.exists(geojsonfilename):
        driver.DeleteDataSource(geojsonfilename)
    datasource = driver.CreateDataSource(geojsonfilename)
    layer = datasource.CreateLayer('buildings', geom_type=ogr.wkbPolygon)
    field_name = ogr.FieldDefn("ImageId", ogr.OFTString)
    field_name.SetWidth(75)
    layer.CreateField(field_name)
    layer.CreateField(ogr.FieldDefn("BuildingId", ogr.OFTInteger))

    # loop through buildings
    for building in buildinglist:
        # create feature
        feature = ogr.Feature(layer.GetLayerDefn())
        feature.SetField("ImageId", building['ImageId'])
        feature.SetField("BuildingId", building['BuildingId'])
        feature.SetGeometry(building['polyPix'])

        # Create the feature in the layer (geojson)
        layer.CreateFeature(feature)
        # Destroy the feature to free resources
        feature.Destroy()

    datasource.Destroy()

    return geojsonfilename 
Example 2
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 6 votes vote down vote up
def mergePolyList(geojsonfilename):

    multipolygon = ogr.Geometry(ogr.wkbMultiPolygon)
    datasource = ogr.Open(geojsonfilename, 0)

    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())


    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:

            multipolygon.AddGeometry(feature.GetGeometryRef().Clone())

    return multipolygon 
Example 3
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 6 votes vote down vote up
def get_envelope(poly):
    env = poly.GetEnvelope()

    # Get Envelope returns a tuple (minX, maxX, minY, maxY)
    # Create ring
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(env[0], env[2])
    ring.AddPoint(env[0], env[3])
    ring.AddPoint(env[1], env[3])
    ring.AddPoint(env[1], env[2])
    ring.AddPoint(env[0], env[2])

    # Create polygon
    poly1 = ogr.Geometry(ogr.wkbPolygon)
    poly1.AddGeometry(ring)

    return poly1 
Example 4
Project: wradlib   Author: wradlib   File: vector.py    (license) View Source Project 6 votes vote down vote up
def ogr_add_geometry(layer, geom, attrs):
    """ Copies single OGR.Geometry object to an OGR.Layer object.

    .. versionadded:: 0.7.0

    Given OGR.Geometry is copied to new OGR.Feature and
    written to given OGR.Layer by given index. Attributes are attached.

    Parameters
    ----------
    layer : OGR.Layer
        object
    geom : OGR.Geometry
        object
    attrs : list
        attributes referring to layer fields

    """
    defn = layer.GetLayerDefn()
    feat = ogr.Feature(defn)

    for i, item in enumerate(attrs):
        feat.SetField(i, item)
    feat.SetGeometry(geom)
    layer.CreateFeature(feat) 
Example 5
Project: wradlib   Author: wradlib   File: vector.py    (license) View Source Project 6 votes vote down vote up
def ogr_to_numpy(ogrobj):
    """Backconvert a gdal/ogr geometry to a numpy vertex array.

    .. versionadded:: 0.7.0

    Using JSON as a vehicle to efficiently deal with numpy arrays.

    Parameters
    ----------
    ogrobj : ogr.Geometry
        object

    Returns
    -------
    out : :class:`numpy:numpy.ndarray`
        a nested ndarray of vertices of shape (num vertices, 2)

    """
    jsonobj = eval(ogrobj.ExportToJson())

    return np.squeeze(jsonobj['coordinates']) 
Example 6
Project: wradlib   Author: wradlib   File: vector.py    (license) View Source Project 6 votes vote down vote up
def get_centroid(polyg):
    """Return centroid of a polygon

    Parameters
    ----------
    polyg : :class:`numpy:numpy.ndarray`
        of shape (num vertices, 2) or ogr.Geometry object

    Returns
    -------
    out : x and y coordinate of the centroid

    """
    if not type(polyg) == ogr.Geometry:
        polyg = numpy_to_ogr(polyg, 'Polygon')
    return polyg.Centroid().GetPoint()[0:2] 
Example 7
Project: beachfront-py   Author: venicegeo   File: mask.py    (license) View Source Project 6 votes vote down vote up
def get_features_as_geojson(layer, bbox=None, union=False):
    """ Get features in this layer and return as GeoJSON """
    if bbox is not None:
        layer.SetSpatialFilterRect(bbox[0], bbox[3], bbox[2], bbox[1])
    poly = ogr.Geometry(ogr.wkbPolygon)
    if union:
        for feature in layer:
            geom = feature.GetGeometryRef()
            #  required for ogr2
            if hasattr(geom, 'GetLinearGeometry'):
                geom = geom.GetLinearGeometry()
            poly = poly.Union(geom)
    if bbox is not None:
        wkt = "POLYGON ((%s %s, %s %s, %s %s, %s %s, %s %s))" % \
              (bbox[0], bbox[1], bbox[2], bbox[1], bbox[2], bbox[3], bbox[0], bbox[3], bbox[0], bbox[1])
        bbox_wkt = ogr.CreateGeometryFromWkt(wkt)
        poly = poly.Intersection(bbox_wkt)
    return json.loads(poly.ExportToJson()) 
Example 8
Project: policosm   Author: ComplexCity   File: addMetricDistanceToEdges.py    (license) View Source Project 6 votes vote down vote up
def addMetricDistanceToEdge(x1,y1,x2,y2, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	line = ogr.Geometry(ogr.wkbLineString)
	line.AddPoint(x1,y1)
	line.AddPoint(x2,y2)
	if epsgCode != 4326:
		line.Transform(coordTrans)
	length = line.Length()

	return length 
Example 9
Project: policosm   Author: ComplexCity   File: addMetricDistanceToEdges.py    (license) View Source Project 6 votes vote down vote up
def addMetricDistanceToEdges(graph, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	metricDistance = {}
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	for edge in graph.edges():
		node1, node2 =  edge
		line = ogr.Geometry(ogr.wkbLineString)
		line.AddPoint(graph.node[node1]['longitude'], graph.node[node1]['latitude'])
		line.AddPoint(graph.node[node2]['longitude'], graph.node[node2]['latitude'])
		if epsgCode != 4326:
			line.Transform(coordTrans)
		length = line.Length()
		metricDistance[edge] = length
	nx.set_edge_attributes(graph, 'length', metricDistance)

	return graph 
Example 10
Project: policosm   Author: ComplexCity   File: testAddMetricDistanceToEdges.py    (license) View Source Project 6 votes vote down vote up
def addMetricDistanceToEdges(graph, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	metricDistance = {}
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	for edge in graph.edges():
		node1, node2 =  edge
		line = ogr.Geometry(ogr.wkbLineString)
		line.AddPoint(graph.node[node1]['longitude'], graph.node[node1]['latitude'])
		line.AddPoint(graph.node[node2]['longitude'], graph.node[node2]['latitude'])
		if epsgCode != 4326:
			line.Transform(coordTrans)
		length = line.Length()
		metricDistance[edge] = length
	nx.set_edge_attributes(graph, 'length', metricDistance)

	return graph 
Example 11
Project: pfb-network-connectivity   Author: azavea   File: detect_utm_zone.py    (license) View Source Project 5 votes vote down vote up
def get_bbox(shapefile):
    """ Gets the bounding box of a shapefile in EPSG 4326.
        If shapefile is not in WGS84, bounds are reprojected.
    """
    driver = ogr.GetDriverByName('ESRI Shapefile')
    # 0 means read, 1 means write
    data_source = driver.Open(shapefile, 0)
    # Check to see if shapefile is found.
    if data_source is None:
        failure('Could not open %s' % (shapefile))
    else:
        layer = data_source.GetLayer()
        shape_bbox = layer.GetExtent()

        spatialRef = layer.GetSpatialRef()
        target = osr.SpatialReference()
        target.ImportFromEPSG(4326)

        # this check for non-WGS84 projections gets some false positives, but that's okay
        if target.ExportToProj4() != spatialRef.ExportToProj4():
            transform = osr.CoordinateTransformation(spatialRef, target)
            point1 = ogr.Geometry(ogr.wkbPoint)
            point1.AddPoint(shape_bbox[0], shape_bbox[2])
            point2 = ogr.Geometry(ogr.wkbPoint)
            point2.AddPoint(shape_bbox[1], shape_bbox[3])
            point1.Transform(transform)
            point2.Transform(transform)
            shape_bbox = [point1.GetX(), point2.GetX(), point1.GetY(), point2.GetY()]

        return shape_bbox 
Example 12
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 5 votes vote down vote up
def readwktcsv(csv_path,removeNoBuildings=True, groundTruthFile=True):
    #
    # csv Format Expected = ['ImageId', 'BuildingId', 'PolygonWKT_Pix', 'PolygonWKT_Geo']
    # returns list of Dictionaries {'ImageId': image_id, 'BuildingId': building_id, 'poly': poly}
    # image_id is a string,
    # BuildingId is an integer,
    # poly is a ogr.Geometry Polygon

    buildinglist = []
    with open(csv_path, 'rb') as csvfile:
        building_reader = csv.reader(csvfile, delimiter=',', quotechar='"')
        next(building_reader, None)  # skip the headers
        for row in building_reader:

            if removeNoBuildings:
                if int(row[1]) != -1:
                    polyPix = ogr.CreateGeometryFromWkt(row[2])
                    if groundTruthFile:
                        polyGeo = ogr.CreateGeometryFromWkt(row[3])
                    else:
                        polyGeo = []
                    buildinglist.append({'ImageId': row[0], 'BuildingId': int(row[1]), 'polyPix': polyPix,
                                         'polyGeo': polyGeo,
                                         })

            else:

                polyPix = ogr.CreateGeometryFromWkt(row[2])
                if groundTruthFile:
                    polyGeo = ogr.CreateGeometryFromWkt(row[3])
                else:
                    polyGeo = []
                buildinglist.append({'ImageId': row[0], 'BuildingId': int(row[1]), 'polyPix': polyPix,
                                     'polyGeo': polyGeo,
                                     })

    return buildinglist 
Example 13
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 5 votes vote down vote up
def latlon2pixel(lat, lon, input_raster='', targetsr='', geom_transform=''):
    # type: (object, object, object, object, object) -> object

    sourcesr = osr.SpatialReference()
    sourcesr.ImportFromEPSG(4326)

    geom = ogr.Geometry(ogr.wkbPoint)
    geom.AddPoint(lon, lat)

    if targetsr == '':
        src_raster = gdal.Open(input_raster)
        targetsr = osr.SpatialReference()
        targetsr.ImportFromWkt(src_raster.GetProjectionRef())
    coord_trans = osr.CoordinateTransformation(sourcesr, targetsr)
    if geom_transform == '':
        src_raster = gdal.Open(input_raster)
        transform = src_raster.GetGeoTransform()
    else:
        transform = geom_transform

    x_origin = transform[0]
    # print(x_origin)
    y_origin = transform[3]
    # print(y_origin)
    pixel_width = transform[1]
    # print(pixel_width)
    pixel_height = transform[5]
    # print(pixel_height)
    geom.Transform(coord_trans)
    # print(geom.GetPoint())
    x_pix = (geom.GetPoint()[0] - x_origin) / pixel_width
    y_pix = (geom.GetPoint()[1] - y_origin) / pixel_height

    return (x_pix, y_pix) 
Example 14
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 5 votes vote down vote up
def returnBoundBox(xOff, yOff, pixDim, inputRaster, targetSR='', pixelSpace=False):
    # Returns Polygon for a specific square defined by a center Pixel and
    # number of pixels in each dimension.
    # Leave targetSR as empty string '' or specify it as a osr.SpatialReference()
    # targetSR = osr.SpatialReference()
    # targetSR.ImportFromEPSG(4326)
    if targetSR == '':
        targetSR = osr.SpatialReference()
        targetSR.ImportFromEPSG(4326)
    xCord = [xOff - pixDim / 2, xOff - pixDim / 2, xOff + pixDim / 2,
             xOff + pixDim / 2, xOff - pixDim / 2]
    yCord = [yOff - pixDim / 2, yOff + pixDim / 2, yOff + pixDim / 2,
             yOff - pixDim / 2, yOff - pixDim / 2]

    ring = ogr.Geometry(ogr.wkbLinearRing)
    for idx in xrange(len(xCord)):
        if pixelSpace == False:
            geom = pixelToGeoCoord(xCord[idx], yCord[idx], inputRaster)
            ring.AddPoint(geom[0], geom[1], 0)
        else:
            ring.AddPoint(xCord[idx], yCord[idx], 0)

    poly = ogr.Geometry(ogr.wkbPolygon)
    if pixelSpace == False:
        poly.AssignSpatialReference(targetSR)

    poly.AddGeometry(ring)

    return poly 
Example 15
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 5 votes vote down vote up
def search_rtree(test_building, index):
    # input test poly ogr.Geometry  and rtree index
    if test_building.GetGeometryName() == 'POLYGON' or \
                    test_building.GetGeometryName() == 'MULTIPOLYGON':
        fidlist = index.intersection(test_building.GetEnvelope())
    else:
        fidlist = []

    return fidlist 
Example 16
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 5 votes vote down vote up
def getRasterExtent(srcImage):
    geoTrans = srcImage.GetGeoTransform()
    ulX = geoTrans[0]
    ulY = geoTrans[3]
    xDist = geoTrans[1]
    yDist = geoTrans[5]
    rtnX = geoTrans[2]
    rtnY = geoTrans[4]

    cols = srcImage.RasterXSize
    rows = srcImage.RasterYSize

    lrX = ulX + xDist * cols
    lrY = ulY + yDist * rows

    # Create ring
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(lrX, lrY)
    ring.AddPoint(lrX, ulY)
    ring.AddPoint(ulX, ulY)
    ring.AddPoint(ulX, lrY)
    ring.AddPoint(lrX, lrY)

    # Create polygon
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)

    return geoTrans, poly, ulX, ulY, lrX, lrY 
Example 17
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 5 votes vote down vote up
def createPolygonFromCorners(lrX,lrY,ulX, ulY):
    # Create ring
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(lrX, lrY)
    ring.AddPoint(lrX, ulY)
    ring.AddPoint(ulX, ulY)
    ring.AddPoint(ulX, lrY)
    ring.AddPoint(lrX, lrY)

    # Create polygon
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)

    return poly 
Example 18
Project: utilities   Author: SpaceNetChallenge   File: labelTools.py    (license) View Source Project 5 votes vote down vote up
def evaluateLineStringPlane(geom, label='Airplane'):
    ring = ogr.Geometry(ogr.wkbLinearRing)

    for i in range(0, geom.GetPointCount()):
        # GetPoint returns a tuple not a Geometry
        pt = geom.GetPoint(i)
        ring.AddPoint(pt[0], pt[1])
    pt = geom.GetPoint(0)
    ring.AddPoint(pt[0], pt[1])
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)

    transform_WGS84_To_UTM, transform_UTM_To_WGS84, utm_cs = gT.createUTMTransform(geom)
    geom.Transform(transform_WGS84_To_UTM)
    pt0 = geom.GetPoint(0) # Tail
    pt1 = geom.GetPoint(1) # Wing
    pt2 = geom.GetPoint(2) # Nose
    pt3 = geom.GetPoint(3) # Wing
    Length = math.sqrt((pt2[0]-pt0[0])**2 + (pt2[1]-pt0[1])**2)
    Width = math.sqrt((pt3[0] - pt1[0])**2 + (pt3[1] - pt1[1])**2)
    Aspect = Length/Width
    Direction = (math.atan2(pt2[0]-pt0[0], pt2[1]-pt0[1])*180/math.pi) % 360


    geom.Transform(transform_UTM_To_WGS84)

    return [poly, Length, Width, Aspect, Direction] 
Example 19
Project: wikiwhere   Author: mkrnr   File: countries.py    (license) View Source Project 5 votes vote down vote up
def __init__(self, lat, lng):
        """ Coordinates are in degrees """
        self.point = ogr.Geometry(ogr.wkbPoint)
        self.point.AddPoint(lng, lat) 
Example 20
Project: chorospy   Author: spyrostheodoridis   File: transFunc.py    (license) View Source Project 5 votes vote down vote up
def transPoint(point, inproj, outproj):
    sr_srs = osr.SpatialReference()
    sr_srs.ImportFromEPSG(inproj)
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromEPSG(outproj)

    coordTransform = osr.CoordinateTransformation(sr_srs, dst_srs)

    gps_point = ogr.Geometry(ogr.wkbPoint)
    gps_point.AddPoint(point[0], point[1])
    gps_point.Transform(coordTransform)
    x = float(gps_point.ExportToWkt().split()[1].split('(')[1])
    y = float(gps_point.ExportToWkt().split()[2])
    return [x,y] 
Example 21
Project: wradlib   Author: wradlib   File: zonalstats.py    (license) View Source Project 5 votes vote down vote up
def get_data_by_geom(self, geom=None):
        """ Returns DataSource geometries filtered by given OGR geometry

        Parameters
        ----------
        geom : OGR.Geometry object
        """
        lyr = self.ds.GetLayer()
        lyr.ResetReading()
        lyr.SetAttributeFilter(None)
        lyr.SetSpatialFilter(geom)
        return self._get_data() 
Example 22
Project: wradlib   Author: wradlib   File: zonalstats.py    (license) View Source Project 5 votes vote down vote up
def _get_intersection(self, trg=None, idx=None, buf=0.):
        """Just a toy function if you want to inspect the intersection
        points/polygons of an arbitrary target or an target by index.
        """
        # TODO: kwargs necessary?

        # check wether idx is given
        if idx is not None:
            if self.trg:
                try:
                    lyr = self.trg.ds.GetLayerByName('trg')
                    feat = lyr.GetFeature(idx)
                    trg = feat.GetGeometryRef()
                except Exception:
                    raise TypeError("No target polygon found at index {0}".
                                    format(idx))
            else:
                raise TypeError('No target polygons found in object!')

        # check for trg
        if trg is None:
            raise TypeError('Either *trg* or *idx* keywords must be given!')

        # check for geometry
        if not type(trg) == ogr.Geometry:
            trg = georef.numpy_to_ogr(trg, 'Polygon')

        # apply Buffer value
        trg = trg.Buffer(buf)

        if idx is None:
            intersecs = self.src.get_data_by_geom(trg)
        else:
            intersecs = self.dst.get_data_by_att('trg_index', idx)

        return intersecs 
Example 23
Project: wradlib   Author: wradlib   File: zonalstats.py    (license) View Source Project 5 votes vote down vote up
def _create_dst_features(self, dst, trg, **kwargs):
        """ Create needed OGR.Features in dst OGR.Layer

        Parameters
        ----------
        dst : OGR.Layer
            destination layer
        trg : OGR.Geometry
            target polygon
        """
        # TODO: kwargs necessary?

        # claim and reset source ogr layer
        layer = self.src.ds.GetLayerByName('src')
        layer.ResetReading()

        # if given, we apply a buffer value to the target polygon filter
        trg_index = trg.GetField('index')
        trg = trg.GetGeometryRef()
        trg = trg.Buffer(self._buffer)
        layer.SetSpatialFilter(trg)

        feat_cnt = layer.GetFeatureCount()

        if feat_cnt:
            [georef.ogr_add_geometry(dst, ogr_src.GetGeometryRef(),
                                     [ogr_src.GetField('index'), trg_index])
             for ogr_src in layer]
        else:
            layer.SetSpatialFilter(None)
            src_pts = np.array([ogr_src.GetGeometryRef().GetPoint_2D(0)
                                for ogr_src in layer])
            centroid = georef.get_centroid(trg)
            tree = cKDTree(src_pts)
            distnext, ixnext = tree.query([centroid[0], centroid[1]], k=1)
            feat = layer.GetFeature(ixnext)
            georef.ogr_add_geometry(dst, feat.GetGeometryRef(),
                                    [feat.GetField('index'), trg_index]) 
Example 24
Project: wradlib   Author: wradlib   File: vector.py    (license) View Source Project 5 votes vote down vote up
def get_shape_points(geom):
    """
    Extract coordinate points from given ogr geometry as generator object

    If geometries are nested, function recurses.

    .. versionadded:: 0.6.0

    Parameters
    ----------
    geom : ogr.Geometry

    Returns
    -------
    result : generator object
        expands to Nx2 dimensional nested point arrays
    """

    type = geom.GetGeometryType()
    if type:
        # 1D Geometries, LINESTRINGS
        if type == 2:
            result = np.array(geom.GetPoints())
            yield result
        # RINGS, POLYGONS, MULTIPOLYGONS, MULTILINESTRINGS
        elif type > 2:
            # iterate over geometries and recurse
            for item in geom:
                for result in get_shape_points(item):
                    yield result
    else:
        print("Unknown Geometry") 
Example 25
Project: wradlib   Author: wradlib   File: vector.py    (license) View Source Project 5 votes vote down vote up
def numpy_to_ogr(vert, geom_name):
    """Convert a vertex array to gdal/ogr geometry.

    .. versionadded:: 0.7.0

    Using JSON as a vehicle to efficiently deal with numpy arrays.

    Parameters
    ----------
    vert : array_like
        a numpy array of vertices of shape (num vertices, 2)
    geom_name : string
        Name of Geometry

    Returns
    -------
    out : ogr.Geometry
        object of type geom_name
    """

    if geom_name in ['Polygon', 'MultiPolygon']:
        json_str = "{{'type':{0!r},'coordinates':[{1!r}]}}".\
            format(geom_name, vert.tolist())
    else:
        json_str = "{{'type':{0!r},'coordinates':{1!r}}}".\
            format(geom_name, vert.tolist())

    return ogr.CreateGeometryFromJson(json_str) 
Example 26
Project: wradlib   Author: wradlib   File: vector.py    (license) View Source Project 5 votes vote down vote up
def ogr_geocol_to_numpy(ogrobj):
    """Backconvert a gdal/ogr geometry Collection to a numpy vertex array.

    .. versionadded:: 0.7.0

    This extracts only Polygon geometries!

    Using JSON as a vehicle to efficiently deal with numpy arrays.

    Parameters
    ----------
    ogrobj : ogr.Geometry
        Collection object

    Returns
    -------
    out : :class:`numpy:numpy.ndarray`
        a nested ndarray of vertices of shape (num vertices, 2)

    """
    jsonobj = eval(ogrobj.ExportToJson())

    mpol = []
    for item in jsonobj['geometries']:
        if item['type'] == 'Polygon':
            mpol.append(item['coordinates'])

    return np.squeeze(mpol) 
Example 27
Project: qgis-geogiglight-plugin   Author: boundlessgeo   File: testgpkg.py    (license) View Source Project 5 votes vote down vote up
def testAddingFeatureUsingOgr(self):
        dataSource =self._getOgrTestLayer()
        layer = dataSource.GetLayer()
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(123, 456)
        featureDefn = layer.GetLayerDefn()
        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(point)
        feature.SetField("fid", 5)
        feature.SetField("n", 5)
        layer.CreateFeature(feature)
        dataSource.Destroy() 
Example 28
Project: pyroSAR   Author: johntruckenbrodt   File: vector.py    (license) View Source Project 5 votes vote down vote up
def bbox(coordinates, crs, outname=None, format='ESRI Shapefile', overwrite=True):
    """
    create a bounding box vector object or shapefile from coordinates and coordinate reference system
    coordinates must be provided in a dictionary containing numerical variables with names 'xmin', 'xmax', 'ymin' and 'ymax'
    the coordinate reference system can be in either WKT, EPSG or PROJ4 format
    """
    srs = crsConvert(crs, 'osr')

    ring = ogr.Geometry(ogr.wkbLinearRing)

    ring.AddPoint(coordinates['xmin'], coordinates['ymin'])
    ring.AddPoint(coordinates['xmin'], coordinates['ymax'])
    ring.AddPoint(coordinates['xmax'], coordinates['ymax'])
    ring.AddPoint(coordinates['xmax'], coordinates['ymin'])
    ring.CloseRings()

    geom = ogr.Geometry(ogr.wkbPolygon)
    geom.AddGeometry(ring)

    geom.FlattenTo2D()

    bbox = Vector(driver='Memory')
    bbox.addlayer('bbox', srs, ogr.wkbPolygon)
    bbox.addfield('id', width=4)
    bbox.addfeature(geom, 'id', 1)
    geom.Destroy()
    if outname is None:
        return bbox
    else:
        bbox.write(outname, format, overwrite) 
Example 29
Project: pyroSAR   Author: johntruckenbrodt   File: bordernoise.py    (license) View Source Project 5 votes vote down vote up
def createPoly(xn, yn, xmax, ymax):
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(0, 0)
    for item in zip(xn, yn):
        item = map(int, item)
        if item != [0, 0] and item != [xmax, ymax]:
            ring.AddPoint(item[0], item[1])
    ring.AddPoint(xmax, ymax)
    ring.AddPoint(xmax, 0)
    ring.CloseRings()
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)
    return poly 
Example 30
Project: pyroSAR   Author: johntruckenbrodt   File: linesimplify.py    (license) View Source Project 5 votes vote down vote up
def createPoly(xn, yn, xmax, ymax):
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(0, 0)
    for item in zip(xn, yn):
        item = map(int, item)
        if item != [0, 0] and item != [xmax, ymax]:
            ring.AddPoint(item[0], item[1])
    ring.AddPoint(xmax, ymax)
    ring.AddPoint(xmax, 0)
    ring.CloseRings()
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)
    return poly 
Example 31
Project: policosm   Author: ComplexCity   File: convertNodesCoordinates.py    (license) View Source Project 5 votes vote down vote up
def convertNodesCoordinates(graph, sourceEPSG, targetEPSG):
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	for node in graph.nodes():
		x, y = graph.node[node]['longitude'], graph.node[node]['latitude']
		point = ogr.Geometry(ogr.wkbPoint)
		point.AddPoint(x,y)
		point.Transform(coordTrans)
		graph.node[node]['longitude'], graph.node[node]['latitude'] = point.GetX(), point.GetY()
	return graph 
Example 32
Project: policosm   Author: ComplexCity   File: convertCoordinates.py    (license) View Source Project 5 votes vote down vote up
def convertCoordinates(coordinate, sourceEPSG, targetEPSG):
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	point = ogr.Geometry(ogr.wkbPoint)
	lon, lat = coordinate
	point.AddPoint(lon,lat)
	point.Transform(coordTrans)
	return point.GetX(), point.GetY() 
Example 33
Project: policosm   Author: ComplexCity   File: convertCoordinates.py    (license) View Source Project 5 votes vote down vote up
def batchConvertCoordinates(sourceCoordinates, sourceEPSG, targetEPSG):
	targetCoordinates = []
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	for lon, lat in sourceCoordinates:
		point = ogr.Geometry(ogr.wkbPoint)
		point.AddPoint(lon,lat)
		point.Transform(coordTrans)
		targetCoordinates.append((point.GetX(), point.GetY()))
	return targetCoordinates 
Example 34
Project: antares   Author: CONABIO   File: split.py    (license) View Source Project 5 votes vote down vote up
def get_convex_hull(shape_name, destination_directory):
    '''
    This method will read all objects from a shape file and create a new one
    with the convex hull of all the geometry points of the first.
    '''
    driver = ogr.GetDriverByName(str('ESRI Shapefile'))
    shape = driver.Open(shape_name, 0)
    layer = shape.GetLayer()
    layer_name = layer.GetName()
    spatial_reference = layer.GetSpatialRef()
    prefix = get_basename(shape_name)
    output_name = create_filename(destination_directory, '%s-hull.shp' % prefix)
    geometries = ogr.Geometry(ogr.wkbGeometryCollection)
    for feature in layer:
        geometries.AddGeometry(feature.GetGeometryRef())
    if os.path.exists(output_name):
        driver.DeleteDataSource(output_name)
    datasource = driver.CreateDataSource(output_name)
    out_layer = datasource.CreateLayer(str('states_convexhull'), spatial_reference, geom_type=ogr.wkbPolygon)
    out_layer.CreateField(ogr.FieldDefn(str('id'), ogr.OFTInteger))
    featureDefn = out_layer.GetLayerDefn()
    feature = ogr.Feature(featureDefn)
    feature.SetGeometry(geometries.ConvexHull())
    feature.SetField(str('id'), 1)
    out_layer.CreateFeature(feature)
    shape.Destroy()
    datasource.Destroy() 
Example 35
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 4 votes vote down vote up
def pixelToGeoCoord(xPix, yPix, inputRaster, sourceSR='', geomTransform='', targetSR=''):
    # If you want to garuntee lon lat output, specify TargetSR  otherwise, geocoords will be in image geo reference
    # targetSR = osr.SpatialReference()
    # targetSR.ImportFromEPSG(4326)
    # Transform can be performed at the polygon level instead of pixel level

    if targetSR =='':
        performReprojection=False
        targetSR = osr.SpatialReference()
        targetSR.ImportFromEPSG(4326)
    else:
        performReprojection=True

    if geomTransform=='':
        srcRaster = gdal.Open(inputRaster)
        geomTransform = srcRaster.GetGeoTransform()

        source_sr = osr.SpatialReference()
        source_sr.ImportFromWkt(srcRaster.GetProjectionRef())




    geom = ogr.Geometry(ogr.wkbPoint)
    xOrigin = geomTransform[0]
    yOrigin = geomTransform[3]
    pixelWidth = geomTransform[1]
    pixelHeight = geomTransform[5]

    xCoord = (xPix * pixelWidth) + xOrigin
    yCoord = (yPix * pixelHeight) + yOrigin
    geom.AddPoint(xCoord, yCoord)


    if performReprojection:
        if sourceSR=='':
            srcRaster = gdal.Open(inputRaster)
            sourceSR = osr.SpatialReference()
            sourceSR.ImportFromWkt(srcRaster.GetProjectionRef())
        coord_trans = osr.CoordinateTransformation(sourceSR, targetSR)
        geom.Transform(coord_trans)

    return (geom.GetX(), geom.GetY()) 
Example 36
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 4 votes vote down vote up
def geoPolygonToPixelPolygonWKT(geom, inputRaster, targetSR, geomTransform, breakMultiPolygonGeo=True,
                                pixPrecision=2):
    # Returns Pixel Coordinate List and GeoCoordinateList

    polygonPixBufferList = []
    polygonPixBufferWKTList = []
    polygonGeoWKTList = []
    if geom.GetGeometryName() == 'POLYGON':
        polygonPix = ogr.Geometry(ogr.wkbPolygon)
        for ring in geom:
            # GetPoint returns a tuple not a Geometry
            ringPix = ogr.Geometry(ogr.wkbLinearRing)

            for pIdx in xrange(ring.GetPointCount()):
                lon, lat, z = ring.GetPoint(pIdx)
                xPix, yPix = latlon2pixel(lat, lon, inputRaster, targetSR, geomTransform)

                xPix = round(xPix, pixPrecision)
                yPix = round(yPix, pixPrecision)
                ringPix.AddPoint(xPix, yPix)

            polygonPix.AddGeometry(ringPix)
        polygonPixBuffer = polygonPix.Buffer(0.0)
        polygonPixBufferList.append([polygonPixBuffer, geom])

    elif geom.GetGeometryName() == 'MULTIPOLYGON':

        for poly in geom:
            polygonPix = ogr.Geometry(ogr.wkbPolygon)
            for ring in poly:
                # GetPoint returns a tuple not a Geometry
                ringPix = ogr.Geometry(ogr.wkbLinearRing)

                for pIdx in xrange(ring.GetPointCount()):
                    lon, lat, z = ring.GetPoint(pIdx)
                    xPix, yPix = latlon2pixel(lat, lon, inputRaster, targetSR, geomTransform)

                    xPix = round(xPix, pixPrecision)
                    yPix = round(yPix, pixPrecision)
                    ringPix.AddPoint(xPix, yPix)

                polygonPix.AddGeometry(ringPix)
            polygonPixBuffer = polygonPix.Buffer(0.0)
            if breakMultiPolygonGeo:
                polygonPixBufferList.append([polygonPixBuffer, poly])
            else:
                polygonPixBufferList.append([polygonPixBuffer, geom])

    for polygonTest in polygonPixBufferList:
        if polygonTest[0].GetGeometryName() == 'POLYGON':
            polygonPixBufferWKTList.append([polygonTest[0].ExportToWkt(), polygonTest[1].ExportToWkt()])
        elif polygonTest[0].GetGeometryName() == 'MULTIPOLYGON':
            for polygonTest2 in polygonTest[0]:
                polygonPixBufferWKTList.append([polygonTest2.ExportToWkt(), polygonTest[1].ExportToWkt()])

    return polygonPixBufferWKTList 
Example 37
Project: rastercube   Author: terrai   File: shputils.py    (license) View Source Project 4 votes vote down vote up
def polygon_to_shapefile(polygons, poly_sr, fname, fields_defs=None,
                         poly_attrs=None):
    """
    Write a set of polygons to a shapefile
    Args:
        polygons: a list of (lat, lng) tuples
        fields_defs: The list of fields for those polygons, as a tuple
                     (name, ogr type) for each field. For example :
                     [('field1', ogr.OFTReal), ('field2', ogr.OFTInteger)]
        poly_attrs: A list of dict containing the attributes for each polygon
                        [{'field1' : 1.0, 'field2': 42},
                         {'field1' : 3.0, 'field2': 60}]
    """
    shp_driver = ogr.GetDriverByName("ESRI Shapefile")
    out_ds = shp_driver.CreateDataSource(fname)
    assert out_ds is not None, "Failed to create temporary %s" % fname
    out_layer = out_ds.CreateLayer(fname, poly_sr, geom_type=ogr.wkbPolygon)

    has_attrs = fields_defs is not None
    if has_attrs:
        attrs_name = []
        for field_name, field_type in fields_defs:
            out_layer.CreateField(ogr.FieldDefn(field_name, field_type))
            attrs_name.append(field_name)

    layer_defn = out_layer.GetLayerDefn()
    for i, poly in enumerate(polygons):
        ring = ogr.Geometry(ogr.wkbLinearRing)
        # gdal uses the (x, y) convention => (lng, lat)
        for point in poly:
            ring.AddPoint(point[1], point[0])
        ring.AddPoint(poly[0][1], poly[0][0])  # re-add the start to close
        p = ogr.Geometry(ogr.wkbPolygon)
        p.AddGeometry(ring)

        out_feature = ogr.Feature(layer_defn)
        out_feature.SetGeometry(p)

        if has_attrs:
            attrs = poly_attrs[i]
            for field_name in attrs_name:
                out_feature.SetField(field_name, attrs[field_name])

        out_layer.CreateFeature(out_feature)

    out_feature.Destroy()
    out_ds.Destroy() 
Example 38
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    (license) View Source Project 4 votes vote down vote up
def filterByCoverage(vectorFile, rasterFile, covPerc):
    
    srcVector = ogr.Open(vectorFile)
    srcLayer = srcVector.GetLayer()
    # merge all features in one geometry (multi polygone)
    multi  = ogr.Geometry(ogr.wkbMultiPolygon)
    for feature in srcLayer:
        geom = feature.GetGeometryRef()
        multi.AddGeometry(geom)
    
    #attributes of raster file
    rasList = raster2array(rasterFile)

    xsize = rasList[4][0]
    ysize = abs(rasList[4][1])

    pixel_area = xsize*ysize

    rows = rasList[0].shape[0]
    cols = rasList[0].shape[1]

    x1 = rasList[2][0]
    y1 = rasList[2][3]
    
    #iterate over raster cells
    for i in range(rows):
        for j in range(cols):
            ring = ogr.Geometry(ogr.wkbLinearRing)

            ring.AddPoint(x1, y1)
            ring.AddPoint(x1 + xsize, y1)
            ring.AddPoint(x1 + xsize, y1 - ysize)
            ring.AddPoint(x1, y1 - ysize)
            ring.AddPoint(x1, y1)

            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            intersect = multi.Intersection(poly)

            if intersect.ExportToWkt() != 'GEOMETRYCOLLECTION EMPTY':
                perc = (intersect.GetArea()/pixel_area)*100
                if perc > covPerc:
                    rasList[0][i][j] = numpy.nan     
            x1 += xsize
        x1 = rasList[2][0]
        y1 -= ysize
    
    return(rasList[0]) #return the filtered array


# numpy array to geo raster 
Example 39
Project: cv4ag   Author: worldbank   File: ogr2ogr.py    (license) View Source Project 4 votes vote down vote up
def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere):
    poGeom = None

    poDS = ogr.Open( pszDS, False )
    if poDS is None:
        return None

    if pszSQL is not None:
        poLyr = poDS.ExecuteSQL( pszSQL, None, None )
    elif pszLyr is not None:
        poLyr = poDS.GetLayerByName(pszLyr)
    else:
        poLyr = poDS.GetLayer(0)

    if poLyr is None:
        print("Failed to identify source layer from datasource.")
        poDS.Destroy()
        return None

    if pszWhere is not None:
        poLyr.SetAttributeFilter(pszWhere)

    poFeat = poLyr.GetNextFeature()
    while poFeat is not None:
        poSrcGeom = poFeat.GetGeometryRef()
        if poSrcGeom is not None:
            eType = wkbFlatten(poSrcGeom.GetGeometryType())

            if poGeom is None:
                poGeom = ogr.Geometry( ogr.wkbMultiPolygon )

            if eType == ogr.wkbPolygon:
                poGeom.AddGeometry( poSrcGeom )
            elif eType == ogr.wkbMultiPolygon:
                for iGeom in range(poSrcGeom.GetGeometryCount()):
                    poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) )

            else:
                print("ERROR: Geometry not of polygon type." )
                if pszSQL is not None:
                    poDS.ReleaseResultSet( poLyr )
                poDS.Destroy()
                return None

        poFeat = poLyr.GetNextFeature()

    if pszSQL is not None:
        poDS.ReleaseResultSet( poLyr )
    poDS.Destroy()

    return poGeom 
Example 40
Project: layman   Author: CCSS-CZ   File: ogr2ogr.py    (license) View Source Project 4 votes vote down vote up
def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere):
    poGeom = None

    poDS = ogr.Open( pszDS, False )
    if poDS is None:
        return None

    if pszSQL is not None:
        poLyr = poDS.ExecuteSQL( pszSQL, None, None )
    elif pszLyr is not None:
        poLyr = poDS.GetLayerByName(pszLyr)
    else:
        poLyr = poDS.GetLayer(0)

    if poLyr is None:
        print("Failed to identify source layer from datasource.")
        poDS.Destroy()
        return None

    if pszWhere is not None:
        poLyr.SetAttributeFilter(pszWhere)

    poFeat = poLyr.GetNextFeature()
    while poFeat is not None:
        poSrcGeom = poFeat.GetGeometryRef()
        if poSrcGeom is not None:
            eType = wkbFlatten(poSrcGeom.GetGeometryType())

            if poGeom is None:
                poGeom = ogr.Geometry( ogr.wkbMultiPolygon )

            if eType == ogr.wkbPolygon:
                poGeom.AddGeometry( poSrcGeom )
            elif eType == ogr.wkbMultiPolygon:
                for iGeom in range(poSrcGeom.GetGeometryCount()):
                    poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) )

            else:
                print("ERROR: Geometry not of polygon type." )
                if pszSQL is not None:
                    poDS.ReleaseResultSet( poLyr )
                poDS.Destroy()
                return None

        poFeat = poLyr.GetNextFeature()

    if pszSQL is not None:
        poDS.ReleaseResultSet( poLyr )
    poDS.Destroy()

    return poGeom 
Example 41
Project: wradlib   Author: wradlib   File: zonalstats.py    (license) View Source Project 4 votes vote down vote up
def _create_dst_features(self, dst, trg, **kwargs):
        """ Create needed OGR.Features in dst OGR.Layer

        Parameters
        ----------
        dst : OGR.Layer
            destination layer
        trg : OGR.Geometry
            target polygon
        """
        # TODO: kwargs necessary?

        # claim and reset source ogr layer
        layer = self.src.ds.GetLayerByName('src')
        layer.ResetReading()

        # if given, we apply a buffer value to the target polygon filter
        trg_index = trg.GetField('index')
        trg = trg.GetGeometryRef()
        trg = trg.Buffer(self._buffer)
        layer.SetSpatialFilter(trg)

        # iterate over layer features
        for ogr_src in layer:
            geom = ogr_src.GetGeometryRef()

            # calculate intersection, if not fully contained
            if not trg.Contains(geom):
                geom = trg.Intersection(geom)

            # checking GeometryCollection, convert to only Polygons,
            #  Multipolygons
            if geom.GetGeometryType() in [7]:
                geocol = georef.ogr_geocol_to_numpy(geom)
                geom = georef.numpy_to_ogr(geocol, 'MultiPolygon')

            # only geometries containing points
            if geom.IsEmpty():
                continue

            if geom.GetGeometryType() in [3, 6, 12]:
                idx = ogr_src.GetField('index')
                georef.ogr_add_geometry(dst, geom, [idx, trg_index]) 
Example 42
Project: pyroSAR   Author: johntruckenbrodt   File: bordernoise.py    (license) View Source Project 4 votes vote down vote up
def crop(seq, maxpoints=20, proximity=100, straighten=False):
    x = map(float, range(0, len(seq)))
    VWpts = simplify(x, seq, maxpoints)
    xn, yn = map(list, zip(*VWpts))
    simple = np.interp(x, xn, yn)
    seq[abs(seq - simple) > proximity] = simple[abs(seq - simple) > proximity]
    points = []
    for xi, yi in enumerate(seq):
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(xi, yi)
        points.append(point)
    points = np.array(points)
    while True:
        poly = createPoly(xn, yn, len(seq), max(seq))
        line = ogr.Geometry(ogr.wkbLineString)
        for xi, yi in zip(xn, yn):
            line.AddPoint(xi, yi)
        dists = np.array([line.Distance(point) for point in points])
        contain = np.array([point.Within(poly) for point in points])
        dists[~contain] = 0
        points = points[(dists > 0)]
        dists = dists[(dists > 0)]
        if len(dists) == 0:
            break
        candidate = points[np.argmax(dists)]
        cp = candidate.GetPoint()
        index = np.argmin(np.array(xn) < cp[0])
        xn.insert(index, cp[0])
        yn.insert(index, cp[1])
    if straighten:
        indices = [i for i in range(0, len(xn)) if (xn[i], yn[i]) in VWpts]
        for i, j in enumerate(indices):
            if i < (len(indices)-1):
                if indices[i+1] > j+1:
                    dx = abs(xn[j] - xn[indices[i + 1]])
                    dy = abs(yn[j] - yn[indices[i + 1]])
                    if dx > dy:
                        seg_y = yn[j:indices[i + 1]+1]
                        print(seg_y)
                        for k in range(j, indices[i + 1]+1):
                            yn[k] = min(seg_y)
    return np.interp(x, xn, yn) 
Example 43
Project: pyroSAR   Author: johntruckenbrodt   File: linesimplify.py    (license) View Source Project 4 votes vote down vote up
def reduce(seq, maxpoints=20, straighten=False):
    if min(seq) == max(seq):
        return np.array(seq)
    x = map(float, range(0, len(seq)))
    # plt.plot(seq)
    VWpts = simplify(x, seq, maxpoints)
    xn, yn = map(list, zip(*VWpts))
    # plt.plot(xn, yn, linewidth=2, color='r')
    simple = np.interp(x, xn, yn)
    seq2 = np.copy(seq)
    points = []
    for xi, yi in enumerate(seq2):
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(xi, yi)
        points.append(point)
    points = np.array(points)
    while True:
        poly = createPoly(xn, yn, len(seq2), max(seq2))
        line = ogr.Geometry(ogr.wkbLineString)
        for xi, yi in zip(xn, yn):
            line.AddPoint(xi, yi)
        dists = np.array([line.Distance(point) for point in points])
        contain = np.array([point.Within(poly) for point in points])
        dists[~contain] = 0
        points = points[(dists > 0)]
        dists = dists[(dists > 0)]
        if len(dists) == 0:
            break
        candidate = points[np.argmax(dists)]
        cp = candidate.GetPoint()
        index = np.argmin(np.array(xn) < cp[0])
        xn.insert(index, cp[0])
        yn.insert(index, cp[1])
    if straighten:
        indices = [i for i in range(0, len(xn)) if (xn[i], yn[i]) in VWpts]
        for i, j in enumerate(indices):
            if i < (len(indices)-1):
                if indices[i+1] > j+1:
                    dx = abs(xn[j] - xn[indices[i + 1]])
                    dy = abs(yn[j] - yn[indices[i + 1]])
                    if dx > dy:
                        seg_y = yn[j:indices[i + 1]+1]
                        for k in range(j, indices[i + 1]+1):
                            yn[k] = min(seg_y)
    # plt.plot(xn, yn, linewidth=2, color='g')
    # plt.show()
    return np.interp(x, xn, yn).astype(int) 
Example 44
Project: unmixing   Author: arthur-e   File: lsma.py    (license) View Source Project 4 votes vote down vote up
def point_to_pixel_geometry(points, source_epsg=None, target_epsg=None, pixel_side_length=30):
    '''
    Where points is a list of X,Y tuples and X and Y are coordinates in
    meters, returns a series of OGR Polygons where each Polygon is the
    pixel extent with a given point at its center. Assumes square pixels.
    Arguments:
        points      Sequence of X,Y numeric pairs or OGR Point geometries
        source_epsg The EPSG code of the source projection (Optional)
        target_epsg The EPSG code of the target projection (Optional)
        pixel_side_length   The length of one side of the intended pixel
    '''
    polys = []
    # Convert points to atomic X,Y pairs if necessary
    if isinstance(points[0], ogr.Geometry):
        points = [(p.GetX(), p.GetY()) for p in points]

    source_ref = target_ref = None
    if all((source_epsg, target_epsg)):
        source_ref = osr.SpatialReference()
        target_ref = osr.SpatialReference()
        source_ref.ImportFromEPSG(source_epsg)
        target_ref.ImportFromEPSG(target_epsg)
        transform = osr.CoordinateTransformation(source_ref, target_ref)

    for p in points:
        r = pixel_side_length / 2 # Half the pixel width
        ring = ogr.Geometry(ogr.wkbLinearRing) # Create a ring
        vertices = [
            (p[0] - r, p[1] + r), # Top-left
            (p[0] + r, p[1] + r), # Top-right, clockwise from here...
            (p[0] + r, p[1] - r),
            (p[0] - r, p[1] - r),
            (p[0] - r, p[1] + r)  # Add top-left again to close ring
        ]

        # Coordinate transformation
        if all((source_ref, target_ref)):
            vertices = [transform.TransformPoint(*v)[0:2] for v in vertices]

        for vert in vertices:
            ring.AddPoint(*vert)
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        polys.append(poly)

    return polys