How to use the ogr.Geometry function in ogr

To help you get started, we’ve selected a few ogr examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github vss-devel / tilers-tools / tilers_tools / tiler_functions.py View on Github external
def sasplanet_hlg2ogr(fname):
    with open(fname) as f:
        lines = f.readlines(4096)
        if not lines[0].startswith('[HIGHLIGHTING]'):
            return None
        coords = [[], []]
        for l in lines[2:]:
            val = float(l.split('=')[1].replace(',','.'))
            coords[1 if 'Lat' in l else 0].append(val)
        points = zip(*coords)
        ld('sasplanet_hlg2ogr', 'points', points)

    ring = ogr.Geometry(ogr.wkbLinearRing)
    for p in points:
        ring.AddPoint(*p)
    polygon = ogr.Geometry(ogr.wkbPolygon)
    polygon.AddGeometry(ring)

    ds = ogr.GetDriverByName('Memory').CreateDataSource( 'wrk' )
    assert ds is not None, 'Unable to create datasource'

    #~ src_srs = txt2srs('EPSG:4326')
    src_srs = txt2srs('+proj=latlong +a=6378137 +b=6378137 +datum=WGS84  +nadgrids=@null +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +no_defs')

    layer = ds.CreateLayer('sasplanet_hlg', srs=src_srs)

    feature = ogr.Feature(layer.GetLayerDefn())
    feature.SetGeometry(polygon)
    layer.CreateFeature(feature)

    del feature
    del polygon
github Kitware / Danesfield / tools / align_buildings.py View on Github external
print("Shifting vector -> {}".format(os.path.basename(outputVectorFile)))
    outVector = outDriver.CreateDataSource(outputVectorFile)
    outSrs = osr.SpatialReference(outProjection)
    # create layer
    outLayer = outVector.CreateLayer(os.path.basename(outputLayerName),
                                     srs=outSrs, geom_type=ogr.wkbPolygon)
    outFeatureDef = outLayer.GetLayerDefn()
    # create rings from input rings by shifting points
    for feature in inputFeatures:
        # create the poly
        outPoly = ogr.Geometry(ogr.wkbPolygon)
        poly = feature.GetGeometryRef()
        for ring_idx in range(poly.GetGeometryCount()):
            ring = poly.GetGeometryRef(ring_idx)
            # create the ring
            outRing = ogr.Geometry(ogr.wkbLinearRing)
            for i in range(0, ring.GetPointCount()):
                pt = ring.GetPoint(i)
                outRing.AddPoint(pt[0] + offsetGeo[0], pt[1] + offsetGeo[1])
            outPoly.AddGeometry(outRing)
        # create feature
        outFeature = ogr.Feature(outFeatureDef)
        outFeature.SetGeometry(outPoly)
        outLayer.CreateFeature(outFeature)
github OpenDataAnalytics / gaia / gaia / geo / gdal_functions.py View on Github external
def world_to_pixel_poly(rpc_dict, geometry):
        """
        Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
        the pixel location of a geospatial coordinate
        """
        pixelRing = ogr.Geometry(ogr.wkbLinearRing)
        geoRing = geometry.GetGeometryRef(0)
        numPoints = geoRing.GetPointCount()
        rpc = rpc_from_gdal_dict(rpc_dict)

        for p in range(numPoints):
            point = numpy.array(map(float, geoRing.GetPoint(p)))
            pixel, line = rpc.project(point)
            pixelRing.AddPoint(pixel, line)

        pixelPoly = ogr.Geometry(ogr.wkbPolygon)
        pixelPoly.AddGeometry(pixelRing)

        return pixelPoly
github brycefrank / pyfor / sampler.py View on Github external
if os.path.exists(path):
            os.remove(path)

        outDriver = ogr.GetDriverByName('ESRI Shapefile')
        outDataSource = outDriver.CreateDataSource(path)
        outLayer = outDataSource.CreateLayer(path, geom_type=ogr.wkbPolygon)
        featureDefn = outLayer.GetLayerDefn()

        self.vert_list = vert_list(xs,ys)
        for vert_set in self.vert_list:
            ring = ogr.Geometry(ogr.wkbLinearRing)
            first_point = vert_set[0]
            for point in vert_set:
                ring.AddPoint(point[0], point[1])
            ring.AddPoint(first_point[0], first_point[1])
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)
            outFeature = ogr.Feature(featureDefn)
            outFeature.SetGeometry(poly)

            outLayer.CreateFeature(outFeature)
            outFeature.Destroy()
        if verts ==True:
            import gisexport
            gisexport.export_coords_to_shp(self.vert_list, "gridverts.shp")
github GEOF-OSGL / CartoLineGen / generalize.py View on Github external
angles[i-1] = segments_angle(points[i-2], points[i-1], points[i])
                if i == p_len-2:    
                    angles[0] = segments_angle(points[p_len-2], points[0], points[1])
                else:    
                    angles[i+2] = segments_angle(points[i+1], points[i+2], points[i+3])
                p_len = p_len + 1
            else:
                angles[i] =  2*MIN_ANGLE  #angle was sharp but segments were too long, skip it in next iteration     
            #if segments are smaller then threshold
        # if smallest angle is on starting point        
        i = np.argmin(angles) #index of sharpest angle
        min_a = np.amin(angles) #sharpest angle
    #while angle is smaller then threshold
        
    if geom_type == 0:
        ring = ogr.Geometry(ogr.wkbLinearRing)    
    else:
        ring = ogr.Geometry(ogr.wkbLineString)    
    for i in range (0,p_len):
        ring.AddPoint(points[i,0], points[i,1])
    if points[0,0]!=points[p_len-1,0] or points[0,1]!=points[p_len-1,1]:
        ring.AddPoint(points[0,0], points[0,1])
    return 1,ring
github GEOF-OSGL / CartoLineGen / generalize.py View on Github external
#output: gen_ring=1 indicates that some geometry is preserved after generalisation
                #output: out_geom will receive all generalised rings in it
                #input: True means that rings are allways closed
                #input: first 0 means that ogrLinearRing is geometry type for creation
                #input: out_geom is reference to geometry in which new generalised geometry of right type will be stored
                #input: alg_type - generalisation (simplification+smoothing), simplification or smoothing
                #input: second 0 means that this is multigeometry (polygon with multiple rings)
                gen_ring,out_geom = Decide(ring,True,0,out_geom,alg_type,0) #perform generalisation

                #there were some geometry preserved, store this in gen so output feature will be created
                #some parts could be deleted due to too small area for given scale
                if gen_ring == 1:
                    gen = 1
                    
        elif geom.GetGeometryName() == 'MULTILINESTRING':
            out_geom = ogr.Geometry(ogr.wkbMultiLineString) #create output geometry of given type
            for i in range(0, geom.GetGeometryCount()): #iterate over lines
                line = geom.GetGeometryRef(i)
                #check if it closed polyline, if so, generalize it as ring, which means 
                #that neither vertex is considered as fixed
                #if line is open, starting and ending points are preserved
                ps = line.GetPoint(0)
                pe = line.GetPoint(line.GetPointCount()-1)
                closed = False
                if (ps[0] == pe[0]) and (ps[1] == pe[1]):
                    closed = True
                    
                #output: gen_line=1 indicates that some geometry is preserved after generalisation
                #output: out_geom will receive all generalised lines in it
                #input: closed indicates whether line is closed or not
                #input: 1 means that ogrLineString is geometry type for creation
                #input: out_geom is reference to geometry in which new generalised geometry of right type will be stored
github gltn / stdm / data / importexport / reader.py View on Github external
def to_ogr_multi_type(self, geom, ogr_type):
        """
        Convert ogr geometry to multi-type when supplied the ogr.type.
        :param geom: The ogr geometry
        :type geom: OGRGeometry
        :param ogr_type: The ogr geometry type
        :type ogr_type: String
        :return: WkB geometry and the output layer geometry type.
        :rtype: Tuple
        """
        multi_geom = ogr.Geometry(ogr_type)
        multi_geom.AddGeometry(geom)
        geom_wkb = multi_geom.ExportToWkt()
        geom_type = multi_geom.GetGeometryName()
        return geom_wkb, geom_type
github sourcepole / ogrtools / ogrtools / pyogr / ogr2ogr.py View on Github external
elif EQUAL(args[iArg], "-clipsrcwhere") and iArg < nArgc - 1:
            pszClipSrcWhere = args[iArg + 1]
            iArg = iArg + 1

        elif EQUAL(args[iArg], "-clipdst") and iArg < nArgc - 1:

            if IsNumber(args[iArg + 1]) and iArg < nArgc - 4:
                oRing = ogr.Geometry(ogr.wkbLinearRing)

                oRing.AddPoint_2D(float(args[iArg + 1]), float(args[iArg + 2]))
                oRing.AddPoint_2D(float(args[iArg + 1]), float(args[iArg + 4]))
                oRing.AddPoint_2D(float(args[iArg + 3]), float(args[iArg + 4]))
                oRing.AddPoint_2D(float(args[iArg + 3]), float(args[iArg + 2]))
                oRing.AddPoint_2D(float(args[iArg + 1]), float(args[iArg + 2]))

                poClipDst = ogr.Geometry(ogr.wkbPolygon)
                poClipDst.AddGeometry(oRing)
                iArg = iArg + 4

            elif (len(args[iArg + 1]) >= 7 and EQUAL(args[iArg + 1][0:7], "POLYGON") ) or \
                    (len(args[iArg + 1]) >= 12 and EQUAL(args[iArg + 1][0:12], "MULTIPOLYGON")):
                poClipDst = ogr.CreateGeometryFromWkt(args[iArg + 1])
                if poClipDst is None:
                    print(
                        "FAILURE: Invalid geometry. Must be a valid POLYGON or MULTIPOLYGON WKT\n")
                    return Usage()

                iArg = iArg + 1

            elif EQUAL(args[iArg + 1], "spat_extent"):
                iArg = iArg + 1
github meetar / dotmap / bin / makedots.py View on Github external
def make_ogr_point(x,y):
	return ogr.Geometry(wkt="POINT(%f %f)"%(x,y))
github GEOF-OSGL / CartoLineGen / generalize.py View on Github external
for i in range(1, p_len):
        p = g.GetPoint(i)
        #remove duplicate neighbour points
        if points[j,0] != p[0] or points[j,1] != p[1]:
            j = j+1
            points[j,0] = p[0]
            points[j,1] = p[1]
    p_len = j+1 #final number of points in point list
    
    #malformed closed geometry
    if p_len < 4:
        return 0,g
    #triangle, nothing to generalize, return cleaned    
    if p_len == 4:    
        if geom_type == 0:
            ring = ogr.Geometry(ogr.wkbLinearRing)    
        else:
            ring = ogr.Geometry(ogr.wkbLineString)    
        for i in range (0,p_len):
            ring.AddPoint(points[i,0], points[i,1])
        if points[0,0]!=points[p_len-1,0] or points[0,1]!=points[p_len-1,1]:
            ring.AddPoint(points[0,0], points[0,1])
        return -1,ring

    #calculate squared segments lengths, algorithm always change line where shorthest segment is found
    segments = np.zeros(p_len-1)
    for i in range(0, p_len-1):
        segments[i] = squared_length(points[i],points[i+1])
           
    i = np.argmin(segments) #index of shorthest segment
    min_s = np.amin(segments) #squared length of shorthest segment
    while min_s < SQR_TOL_LENGTH and p_len > 4: