Python osgeo.ogr.Geometry() Examples

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) 

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))) 

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 

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 

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 

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) 

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) 

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) 

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) 

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"]) 

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] 

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) 

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) 

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 

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"
        ) 

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) 

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 

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 

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]) 

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 

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() 

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 ************************************************************* ** 

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) 

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)


############################################################################### 

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)


############################################################################### 

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) 

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 

def get_poly_intersection(poly1, poly2):
    """A functional wrapper for ogr.Geometry.Intersection()"""
    return poly1.Intersection(poly2) 

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 

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