Python osgeo.ogr.Geometry() Examples

The following are 30 code examples for showing how to use osgeo.ogr.Geometry(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module osgeo.ogr , or try the search function .

Example 1
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def multiple_union(polygons):
    """
    Takes a list of polygons and returns a polygon of the union of their perimeter

    Parameters
    ----------
    polygons
        A list of ogr.Geometry objects, each containing a single polygon.

    Returns
    -------
        An ogr.Geometry object containing a single polygon

    """
    # Note; I can see this maybe failing(or at least returning a multipolygon)
    # if two consecutive polygons do not overlap at all. Keep eye on.
    running_union = polygons[0]
    for polygon in polygons[1:]:
        running_union = running_union.Union(polygon)
    return running_union.Simplify(0) 
Example 2
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def multiple_intersection(polygons):
    """
    Takes a list of polygons and returns a geometry representing the intersection of all of them

    Parameters
    ----------
    polygons
        A list of ogr.Geometry objects, each containing a single polygon.

    Returns
    -------
        An ogr.Geometry object containing a single polygon
    """
    running_intersection = polygons[0]
    for polygon in polygons[1:]:
        running_intersection = running_intersection.Intersection(polygon)
    return running_intersection.Simplify(0) 
Example 3
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def get_aoi_intersection(raster, aoi):
    """
    Returns a wkbPolygon geometry with the intersection of a raster and a shpefile containing an area of interest

    Parameters
    ----------
    raster
        A raster containing image data
    aoi
        A shapefile with a single layer and feature
    Returns
    -------
    a ogr.Geometry object containing a single polygon with the area of intersection

    """
    raster_shape = get_raster_bounds(raster)
    aoi.GetLayer(0).ResetReading()  # Just in case the aoi has been accessed by something else
    aoi_feature = aoi.GetLayer(0).GetFeature(0)
    aoi_geometry = aoi_feature.GetGeometryRef()
    return aoi_geometry.Intersection(raster_shape) 
Example 4
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def get_raster_intersection(raster1, raster2):
    """
    Returns a wkbPolygon geometry with the intersection of two raster bounding boxes.

    Parameters
    ----------
    raster1, raster2
        A gdal.Image() object
    Returns
    -------
        a ogr.Geometry object containing a single polygon

    """
    bounds_1 = get_raster_bounds(raster1)
    bounds_2 = get_raster_bounds(raster2)
    return bounds_1.Intersection(bounds_2) 
Example 5
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def get_aoi_bounds(aoi):
    """
    Returns a wkbPolygon geometry with the bounding rectangle of a single-polygon shapefile

    Parameters
    ----------
    aoi
        An ogr.Dataset object containing a single layer.

    Returns
    -------

    """
    aoi_bounds = ogr.Geometry(ogr.wkbLinearRing)
    (x_min, x_max, y_min, y_max) = aoi.GetLayer(0).GetExtent()
    aoi_bounds.AddPoint(x_min, y_min)
    aoi_bounds.AddPoint(x_max, y_min)
    aoi_bounds.AddPoint(x_max, y_max)
    aoi_bounds.AddPoint(x_min, y_max)
    aoi_bounds.AddPoint(x_min, y_min)
    bounds_poly = ogr.Geometry(ogr.wkbPolygon)
    bounds_poly.AddGeometry(aoi_bounds)
    return bounds_poly 
Example 6
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def get_poly_size(poly):
    """
    Returns the width and height of a bounding box of a polygon

    Parameters
    ----------
    poly
        A ogr.Geometry object containing the polygon.

    Returns
    -------
    A tuple of (width, height).

    """
    boundary = poly.Boundary()
    x_min, y_min, not_needed = boundary.GetPoint(0)
    x_max, y_max, not_needed = boundary.GetPoint(2)
    out = (x_max - x_min, y_max-y_min)
    return out 
Example 7
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def get_poly_bounding_rect(poly):
    """
    Returns a polygon of the bounding rectangle of an input polygon.
    Parameters
    ----------
    poly
        An ogr.Geometry object containing a polygon

    Returns
    -------
    An ogr.Geometry object with a four-point polygon representing the bounding rectangle.

    """
    aoi_bounds = ogr.Geometry(ogr.wkbLinearRing)
    x_min, x_max, y_min, y_max = poly.GetEnvelope()
    aoi_bounds.AddPoint(x_min, y_min)
    aoi_bounds.AddPoint(x_max, y_min)
    aoi_bounds.AddPoint(x_max, y_max)
    aoi_bounds.AddPoint(x_min, y_max)
    aoi_bounds.AddPoint(x_min, y_min)
    bounds_poly = ogr.Geometry(ogr.wkbPolygon)
    bounds_poly.AddGeometry(aoi_bounds)
    return bounds_poly 
Example 8
Project: osgeopy-code   Author: cgarrard   File: vectorplotter.py    License: MIT License 6 votes vote down vote up
def plot(self, geom_or_lyr, symbol='', name='', **kwargs):
        """Plot a geometry or layer.
        geom_or_lyr - geometry, layer, or filename
        symbol      - optional pyplot symbol to draw the geometries with
        name        - optional name to assign to layer so can access it later
        kwargs      - optional pyplot drawing parameters
        """
        if type(geom_or_lyr) is str:
            lyr, ds = pb._get_layer(geom_or_lyr)
            self.plot_layer(lyr, symbol, name, **kwargs)
        elif type(geom_or_lyr) is ogr.Geometry:
            self.plot_geom(geom_or_lyr, symbol, name, **kwargs)
        elif type(geom_or_lyr) is ogr.Layer:
            self.plot_layer(geom_or_lyr, symbol, name, **kwargs)
        else:
            raise RuntimeError('{} is not supported.'.format(type(geom_or_lyr))) 
Example 9
Project: wradlib   Author: wradlib   File: vector.py    License: MIT License 6 votes vote down vote up
def ogr_add_geometry(layer, geom, attrs):
    """ Copies single OGR.Geometry object to an OGR.Layer object.

    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 10
Project: wradlib   Author: wradlib   File: vector.py    License: MIT License 6 votes vote down vote up
def ogr_to_numpy(ogrobj):
    """Backconvert a gdal/ogr geometry to a numpy vertex array.

    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 11
Project: wradlib   Author: wradlib   File: vector.py    License: MIT License 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 12
Project: SMAC-M   Author: LarsSchy   File: georeference_tif.py    License: MIT License 6 votes vote down vote up
def main():
    args = parse_arguments()
    src_ds = gdal.Open(args.src_file[0])

    driver = gdal.GetDriverByName("GTiff")

    dst_ds = driver.CreateCopy(args.out_file[0], src_ds, 0)

    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(float(args.position[0]), float(args.position[1]))
    point.Transform(transform)

    gt = [point.GetX(), float(args.pixelsize[0]), 0,
          point.GetY(), 0, -float(args.pixelsize[1])]

    dst_ds.SetGeoTransform(gt)

    to_sys = osr.SpatialReference()
    to_sys.ImportFromEPSG(to_srs)

    dest_wkt = to_sys.ExportToWkt()

    dst_ds.SetProjection(dest_wkt) 
Example 13
Project: buzzard   Author: airware   File: _footprint.py    License: Apache License 2.0 5 votes vote down vote up
def raster_to_spatial(self, xy):
        """Convert xy raster coordinates to spatial coordinates

        Parameters
        ----------
        xy: sequence of numbers of shape (..., 2)
           Raster coordinages

        Returns
        -------
        out_xy: np.ndarray
            Spatial coordinates
            with shape = np.asarray(xy).shape
            with dtype = dtype

        """
        # Check xy parameter
        xy = np.asarray(xy)
        if xy.shape[-1] != 2:
            raise ValueError('An array of shape (..., 2) was expected') # pragma: no cover

        workshape = int(xy.size / 2), 2
        xy2 = np.empty(workshape, 'float64')
        xy2[:, :] = xy.reshape(workshape)
        aff = self._aff
        xy2[:, 0], xy2[:, 1] = (
            xy2[:, 0] * aff.a + xy2[:, 1] * aff.b + aff.c,
            xy2[:, 0] * aff.d + xy2[:, 1] * aff.e + aff.f,
        )
        return xy2.reshape(xy.shape)

    # Geometry / Raster conversions ************************************************************* ** 
Example 14
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 5 votes vote down vote up
def get_combined_polygon(rasters, geometry_mode ="intersect"):
    """
    Returns a polygon containing the combined boundary of each raster in rasters.

    Parameters
    ----------
    rasters
        A list of raster objects opened with gdal.Open()
    geometry_mode
        If 'intersect', returns the boundary of the area that all rasters cover.
        If 'union', returns the boundary of the area that any raster covers.

    Returns
    -------
        ogr.Geometry() containing a polygon.

    """
    raster_bounds = []
    for in_raster in rasters:
        raster_bounds.append(get_raster_bounds(in_raster))
    # Calculate overall bounding box based on either union or intersection of rasters
    if geometry_mode == "intersect":
        combined_polygons = multiple_intersection(raster_bounds)
    elif geometry_mode == "union":
        combined_polygons = multiple_union(raster_bounds)
    else:
        raise Exception("Invalid geometry mode")
    return combined_polygons 
Example 15
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 5 votes vote down vote up
def point_to_pixel_coordinates(raster, point, oob_fail=False):
    """
    Returns a tuple (x_pixel, y_pixel) in a georaster raster corresponding to the geographic point in a projection.
     Assumes raster is north-up non rotated.

    Parameters
    ----------
    raster
        A gdal raster object
    point
        One of:
            A well-known text string of a single point
            An iterable of the form (x,y)
            An ogr.Geometry object containing a single point
    Returns
    -------
    A tuple of (x_pixel, y_pixel), containing the indicies of the point in the raster.

    Notes
    -----
    The equation is a rearrangement of the section on affinine geotransform in http://www.gdal.org/gdal_datamodel.html

    """
    if isinstance(point, str):
        point = ogr.CreateGeometryFromWkt(point)
        x_geo = point.GetX()
        y_geo = point.GetY()
    if isinstance(point, list) or isinstance(point, tuple):  # There is a more pythonic way to do this
        x_geo = point[0]
        y_geo = point[1]
    if isinstance(point, ogr.Geometry):
        x_geo = point.GetX()
        y_geo = point.GetY()
    gt = raster.GetGeoTransform()
    x_pixel = int(np.floor((x_geo - floor_to_resolution(gt[0], gt[1]))/gt[1]))
    y_pixel = int(np.floor((y_geo - floor_to_resolution(gt[3], gt[5]*-1))/gt[5]))  # y resolution is -ve
    return x_pixel, y_pixel 
Example 16
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 5 votes vote down vote up
def write_geometry(geometry, out_path, srs_id=4326):
    """
    Saves the geometry in an ogr.Geometry object to a shapefile.

    Parameters
    ----------
    geometry
        An ogr.Geometry object
    out_path
        The location to save the output shapefile
    srs_id
        The projection of the output shapefile. Can be an EPSG number or a WKT string.

    Notes
    -----
    The shapefile consists of one layer named 'geometry'.


    """
    # TODO: Fix this needing an extra filepath on the end
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(out_path)
    srs = osr.SpatialReference()
    if type(srs_id) is int:
        srs.ImportFromEPSG(srs_id)
    if type(srs_id) is str:
        srs.ImportFromWkt(srs_id)
    layer = data_source.CreateLayer(
        "geometry",
        srs,
        geom_type=geometry.GetGeometryType())
    feature_def = layer.GetLayerDefn()
    feature = ogr.Feature(feature_def)
    feature.SetGeometry(geometry)
    layer.CreateFeature(feature)
    data_source.FlushCache()
    data_source = None 
Example 17
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 5 votes vote down vote up
def get_poly_intersection(poly1, poly2):
    """A functional wrapper for ogr.Geometry.Intersection()"""
    return poly1.Intersection(poly2) 
Example 18
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 5 votes vote down vote up
def align_bounds_to_whole_number(bounding_box):
    """
    Creates a new bounding box with it's height and width rounded to the nearest whole number.

    Parameters
    ----------
    bounding_box
        An ogr.Geometry object containing a raster's bounding box as a polygon.

    Returns
    -------
    An ogr.Geometry object containing the aligned bounding box.

    """
    # This should prevent off-by-one errors caused by bad image alignment
    aoi_bounds = ogr.Geometry(ogr.wkbLinearRing)
    (x_min, x_max, y_min, y_max) = bounding_box.GetEnvelope()
    # This will create a box that has a whole number as its height and width
    x_new = x_min + np.floor(x_max-x_min)
    y_new = y_min + np.floor(y_max-y_min)
    aoi_bounds.AddPoint(x_min, y_min)
    aoi_bounds.AddPoint(x_new, y_min)
    aoi_bounds.AddPoint(x_new, y_new)
    aoi_bounds.AddPoint(x_min, y_new)
    aoi_bounds.AddPoint(x_min, y_min)
    bounds_poly = ogr.Geometry(ogr.wkbPolygon)
    bounds_poly.AddGeometry(aoi_bounds)
    return bounds_poly 
Example 19
Project: apls   Author: CosmiQ   File: apls_utils.py    License: Apache License 2.0 5 votes vote down vote up
def _latlon2pixel(lat, lon, input_raster='', targetsr='', geom_transform=''):
    '''
    Convert latitude, longitude coords to pixexl coords.
    From spacenet geotools
    '''

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

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

    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 20
Project: apls   Author: CosmiQ   File: apls_utils.py    License: Apache License 2.0 5 votes vote down vote up
def _wmp2pixel(x, y, input_raster='', targetsr='', geom_transform=''):
    '''
    Convert wmp coords to pixexl coords.
    '''

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

    geom = ogr.Geometry(ogr.wkbPoint)
    geom.AddPoint(x, y)

    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 21
Project: geoio   Author: DigitalGlobe   File: base.py    License: MIT License 5 votes vote down vote up
def get_data_from_coords(self, coords, **kwargs):
        """Method to get pixel data given polygon coordintes in the same projection as
        the image. kwargs can be anything accepted by get_data.

        Parameters
        ----------
        coords : list of lists. polygon coordinates formatted as follows:
            - lat/long (EPSG:4326) projection: [[lon_1, lat_1], [lon_2, lat_2], ...]
            - UTM projection: [[x_1, y_1], [x_2, y_2], ...]

        Returns
        ------
        ndarray
            Three dimensional numpy array of data from region of the image specified
            in the feature geometry.
        """

        # create ogr geometry ring from feature geom
        ring = ogr.Geometry(ogr.wkbLinearRing)

        for point in coords:
            lon, lat = point[0], point[1]
            ring.AddPoint(lon, lat)

        ring.AddPoint(coords[0][0], coords[0][1]) # close geom ring with first point
        geom = ogr.Geometry(ogr.wkbPolygon)
        geom.AddGeometry(ring)

        return self.get_data(geom=geom, **kwargs) 
Example 22
Project: wradlib   Author: wradlib   File: zonalstats.py    License: MIT License 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 23
Project: wradlib   Author: wradlib   File: zonalstats.py    License: MIT License 5 votes vote down vote up
def _get_intersection(self, trg=None, idx=None, buf=0.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.vector.numpy_to_ogr(trg, "Polygon")

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

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

        return intersecs 
Example 24
Project: wradlib   Author: wradlib   File: test_georef.py    License: MIT License 5 votes vote down vote up
def test_ogr_add_geometry(self):
        ds = wradlib.io.gdal_create_dataset("Memory", "test", gdal_type=gdal.OF_VECTOR)
        lyr = georef.ogr_create_layer(
            ds, "test", geom_type=ogr.wkbPoint, fields=[("test", ogr.OFTReal)]
        )
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(1198054.34, 648493.09)
        georef.ogr_add_geometry(lyr, point, [42.42]) 
Example 25
Project: wradlib   Author: wradlib   File: vector.py    License: MIT License 5 votes vote down vote up
def get_vector_points(geom):
    """Extract coordinate points from given ogr geometry as generator object

    If geometries are nested, function recurses.

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

    Returns
    -------
    result : generator object
        expands to Nx2 dimensional nested point arrays
    """
    geomtype = geom.GetGeometryType()
    if geomtype > 1:
        # 1D Geometries, LINESTRINGS
        if geomtype == 2:
            result = np.array(geom.GetPoints())
            yield result
        # RINGS, POLYGONS, MULTIPOLYGONS, MULTILINESTRINGS
        elif geomtype > 2:
            # iterate over geometries and recurse
            for item in geom:
                for result in get_vector_points(item):
                    yield result
    else:
        warnings.warn(
            "unsupported geometry type detected in "
            "wradlib.georef.get_vector_points - skipping"
        ) 
Example 26
Project: wradlib   Author: wradlib   File: vector.py    License: MIT License 5 votes vote down vote up
def numpy_to_ogr(vert, geom_name):
    """Convert a vertex array to gdal/ogr geometry.

    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 27
Project: wradlib   Author: wradlib   File: vector.py    License: MIT License 5 votes vote down vote up
def ogr_geocol_to_numpy(ogrobj):
    """Backconvert a gdal/ogr geometry Collection to a numpy vertex array.

    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"]:
        print(item["type"])
        if item["type"] == "Polygon":
            mpol.append(item["coordinates"])

    return np.squeeze(mpol) 
Example 28
Project: stdm   Author: gltn   File: reader.py    License: GNU General Public License v2.0 5 votes vote down vote up
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 
Example 29
Project: qgisSpaceSyntaxToolkit   Author: SpaceGroupUCL   File: test_shp.py    License: GNU General Public License v3.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 30
Project: RHEAS   Author: nasa   File: modis.py    License: MIT License 5 votes vote down vote up
def findTiles(bbox):
    """Returns the tile IDs that need to be downloaded for
    a given region bounded by *bbox*."""
    log = logging.getLogger(__name__)

    def intersects(bbox, tile):
        if tile[2] != -999.0 and tile[3] != -999.0 and tile[4] != -99.0 and tile[5] != -99.0:
            tiler = ogr.Geometry(ogr.wkbLinearRing)
            tiler.AddPoint(tile[2], tile[4])
            tiler.AddPoint(tile[3], tile[4])
            tiler.AddPoint(tile[3], tile[5])
            tiler.AddPoint(tile[2], tile[5])
            tiler.AddPoint(tile[2], tile[4])
            polyr = ogr.Geometry(ogr.wkbPolygon)
            polyr.AddGeometry(tiler)
            bboxr = ogr.Geometry(ogr.wkbLinearRing)
            bboxr.AddPoint(bbox[0], bbox[1])
            bboxr.AddPoint(bbox[2], bbox[1])
            bboxr.AddPoint(bbox[2], bbox[3])
            bboxr.AddPoint(bbox[0], bbox[3])
            bboxr.AddPoint(bbox[0], bbox[1])
            polyb = ogr.Geometry(ogr.wkbPolygon)
            polyb.AddGeometry(bboxr)
            return polyr.Intersects(polyb)
        else:
            return False
    if bbox is None:
        log.warning("No bounding box provided for MODIS dataset. Skipping download!")
        ids = None
    else:
        ids = [(t[0], t[1]) for t in tiles if intersects(bbox, t)]
    return ids