How to use the gdal.InvGeoTransform function in GDAL

To help you get started, we’ve selected a few GDAL 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 nansencenter / nansat / nansat / vrt.py View on Github external
def _update_warped_vrt_xml(self, x_size, y_size, geo_transform, block_size, working_data_type):
        """Update rasterXsize, rasterYsize, geotransform, block_size and working_data_type"""
        node0 = Node.create(str(self.xml))
        node0.replaceAttribute('rasterXSize', str(x_size))
        node0.replaceAttribute('rasterYSize', str(y_size))

        if geo_transform is not None:
            invGeotransform = gdal.InvGeoTransform(geo_transform)
            # convert proper string style and set to the GeoTransform element
            node0.node('GeoTransform').value = str(geo_transform).strip('()')
            node0.node('DstGeoTransform').value = str(geo_transform).strip('()')
            node0.node('DstInvGeoTransform').value = (
                str(invGeotransform[1]).strip('()'))

            if node0.node('SrcGeoLocTransformer'):
                node0.node('BlockXSize').value = str(x_size)
                node0.node('BlockYSize').value = str(y_size)

            if block_size is not None:
                node0.node('BlockXSize').value = str(block_size)
                node0.node('BlockYSize').value = str(block_size)

            if working_data_type is not None:
                node0.node('WorkingDataType').value = working_data_type
github vss-devel / tilers-tools / tilers_tools / tiler_backend.py View on Github external
def create_warped_vrt(self, top_left_coord, res, size):

    #----------------------------

        # generate warp transform
        src_geotr = self.src_ds.GetGeoTransform()
        src_proj = txt2proj4(self.src_ds.GetProjection())
        gcp_proj = None

        if not self.options.tps and src_geotr and src_geotr != (0.0, 1.0, 0.0, 0.0, 0.0, 1.0):
            src_igeotr = gdal.InvGeoTransform(src_geotr)
            src_transform = '%s\n%s' % (warp_src_geotr % src_geotr, warp_src_igeotr % src_igeotr)
        else:
            gcps = self.src_ds.GetGCPs()
            assert gcps, 'Neither geotransform, nor gpcs are in the source file %s' % self.src

            gcp_lst = [(g.Id, g.GCPPixel, g.GCPLine, g.GCPX, g.GCPY, g.GCPZ) for g in gcps]
            ld('src_proj', self.src_ds.GetProjection(), 'gcp_proj', self.src_ds.GetGCPProjection())
            gcp_proj = txt2proj4(self.src_ds.GetGCPProjection())
            if src_proj and gcp_proj != src_proj:
                coords = GdalTransformer(
                    SRC_SRS=gcp_proj,
                    DST_SRS=src_proj
                    ).transform([g[3:6] for g in gcp_lst])

                gcp_lst = [tuple(p[:3] + c) for p, c in zip(gcp_lst, coords)]
github buddebej / ol3-dem / dem-preprocessing / tiler-tools / tiler_backend.py View on Github external
ld('src_proj', self.src_ds.GetProjection(), 'gcp_proj', self.src_ds.GetGCPProjection())
            gcp_proj = txt2proj4(self.src_ds.GetGCPProjection())
            if src_proj and gcp_proj != src_proj:
                coords = GdalTransformer(SRC_SRS=gcp_proj, DST_SRS=src_proj).transform([g[3:6] for g in gcp_lst])
                gcp_lst = [tuple(p[:3]+c) for p, c in zip(gcp_lst, coords)]

            gcp_txt = '\n'.join((gcp_templ % g for g in gcp_lst))
            #src_transform = warp_src_gcp_transformer % (0, gcp_txt)
            src_transform = warp_src_tps_transformer % gcp_txt

        res = self.zoom2res(zoom)
        #ul_ll, lr_ll = self.coords2longlat([ul_c, lr_c])
        ld('max_zoom', zoom, 'size', dst_xsize, dst_ysize, '-tr', res[0], res[1], '-te', ul_c[0], lr_c[1], lr_c[0], ul_c[1], '-t_srs', self.proj_srs)
        dst_geotr = ( ul_c[0], res[0], 0.0,
                    ul_c[1], 0.0, res[1] )
        ok, dst_igeotr = gdal.InvGeoTransform(dst_geotr)
        assert ok
        dst_transform = '%s\n%s' % (warp_dst_geotr % dst_geotr, warp_dst_igeotr % dst_igeotr)

        # generate warp options
        warp_options = []
        def w_option(name, value): # warp options template
            return '    <option name="%s">%s</option>' % (name, value)

        warp_options.append(w_option('INIT_DEST', 'NO_DATA'))

        # generate cut line
        if self.options.cut or self.options.cutline:
            cut_wkt = self.get_cutline()
        else:
            cut_wkt = None
        if cut_wkt:
github GeoscienceAustralia / agdc / api / source / main / python / datacube / api / utils.py View on Github external
def latlon_to_xy(lat, lon, transform):
    """
    Convert lat/lon to x/y for raster
    NOTE: No projection done - assumes raster has native lat/lon projection

    :param lat: latitude
    :param lon: longitude
    :param transform: GDAL GeoTransform
    :return: x, y pair
    """
    # Get the reverse direction GeoTransform
    _, transform = gdal.InvGeoTransform(transform)
    
    '''
    ulx, uly = transform[0], transform[3]
    psx, psy = transform[1], transform[5]
    x = int(math.floor(ulx + psx * lon))
    y = int(math.floor(uly + psy * lat))
    '''
    # This is to fix to make nearest rounding before casting to integer and apply gdal utility
    x, y= gdal.ApplyGeoTransform(transform, lon, lat)
    return int(round(x)), int(round(y))
github Wireless-Innovation-Forum / Spectrum-Access-System / src / data / nlcd_origin.py View on Github external
"""Initializes the Tile Information.

    Inputs:
      tile_path: path the tile file.
    """
    self.tile_path = tile_path

    ds = gdal.Open(tile_path)
    self.width = ds.RasterXSize
    self.height = ds.RasterYSize

    # Gets the gdal geo transformation tuples
    gdal_version = osgeo.gdal.__version__
    self._txf = ds.GetGeoTransform()
    if gdal_version[0] == '1':
      self._inv_txf = gdal.InvGeoTransform(self._txf)[1]
    else:
      self._inv_txf = gdal.InvGeoTransform(self._txf)
    # Gets the transformation from lat/lon to coordinates
    wgs84_ref = osr.SpatialReference()
    wgs84_ref.ImportFromEPSG(4326)   # WGS84
    sref = osr.SpatialReference()
    sref.ImportFromWkt(ds.GetProjection())
    self._transform = osr.CoordinateTransformation(wgs84_ref, sref)
    inv_transform = osr.CoordinateTransformation(sref, wgs84_ref)
    # Find a loose lat/lon bounding box  for quick check without
    # having to do full coordinates transformation
    corners = []
    for x in [0, self.width]:
      for y in [0, self.height]:
        corners.append([self._txf[0] + self._txf[1] * x + self._txf[2] * y,
                        self._txf[3] + self._txf[4] * x + self._txf[5] * y])
github Wireless-Innovation-Forum / Spectrum-Access-System / src / geo / ned_indexer.py View on Github external
def LoadTileForLatLng(self, lat, lng):
    latf = int(math.ceil(lat))
    lngf = int(math.floor(lng))
    k = '%s.%s' % (latf, lngf)
    if k in self.tile_cache:
      return

    filename = self.latlng_file[k]
    print 'Loading tile %s from %s' % (k, filename)
    dataset = gdal.Open(filename)
    # Store the inverse geo transform which will map (lat, lng) to array indices
    # in this tile.
    tx = dataset.GetGeoTransform()
    self.txf[k] = gdal.InvGeoTransform(tx)

    self.tile_cache[k] = dataset.ReadAsArray().astype(numpy.float)
    self.tile_lru[k] = time.clock()

    # evict a tile if there are more than tile_lru_size elements in the cache
    if len(self.tile_lru) > self.tile_lru_size:
      mint = 0
      mink = ''
      for k in self.tile_lru.keys():
        if self.tile_lru[k] > mint:
          mint = self.tile_lru[k]
          mink = k
      print 'Evicting tile %s' % k
      del self.tile_cache[k]
      del self.tile_lru[k]
github vss-devel / tilers-tools / tilers_tools / tiler_backend.py View on Github external
).transform([g[3:6] for g in gcp_lst])

                gcp_lst = [tuple(p[:3] + c) for p, c in zip(gcp_lst, coords)]

            gcp_txt = '\n'.join((gcp_templ % g for g in gcp_lst))
            #src_transform = warp_src_gcp_transformer % (0, gcp_txt)
            src_transform = warp_src_tps_transformer % gcp_txt


        #tl_ll, br_ll = self.coords2longlat([tl_c, br_c])
        #~ ld('create_target_dataset', 'max_zoom', self.max_zoom, 'size', size[0], size[1], '-tr', res[0], res[1], '-te', tl_c[0], br_c[1], br_c[0], tl_c[1], '-t_srs', self.proj_srs)

        left, top  = top_left_coord
        dst_geotr = ( left, res[0], 0.0,
                      top, 0.0, -res[1] )
        dst_igeotr = gdal.InvGeoTransform(dst_geotr)
        dst_transform = '%s\n%s' % (warp_dst_geotr % dst_geotr, warp_dst_igeotr % dst_igeotr)

        # generate warp options
        warp_options = []
        def w_option(name, value): # warp options template
            return '    <option name="%s">%s</option>' % (name, value)

        warp_options.append(w_option('INIT_DEST', 'NO_DATA'))

        # generate cut line
        if self.options.cut or self.options.cutline:
            cut_wkt = self.get_cutline()
        else:
            cut_wkt = None
        if cut_wkt:
            warp_options.append(w_option('CUTLINE', cut_wkt))
github buddebej / dem2tiles / tiler-tools / tiler_backend.py View on Github external
# base zoom level raster size
        dst_xsize = lr_pix[0]-ul_pix[0]
        dst_ysize = lr_pix[1]-ul_pix[1]

        ld('target Upper Left', self.bounds[0], ul_c, self.proj2geog.transform([self.bounds[0], ul_c]))
        ld('target Lower Right', self.bounds[1], lr_c, self.proj2geog.transform([self.bounds[1], lr_c]))

        # create VRT for base image warp

        # generate warp transform
        src_geotr = self.src_ds.GetGeoTransform()
        src_proj = txt2proj4(self.src_ds.GetProjection())
        gcp_proj = None

        if not self.options.tps and src_geotr and src_geotr != (0.0, 1.0, 0.0, 0.0, 0.0, 1.0):
            ok, src_igeotr = gdal.InvGeoTransform(src_geotr)
            assert ok
            src_transform = '%s\n%s' % (warp_src_geotr % src_geotr, warp_src_igeotr % src_igeotr)
        else:
            gcps = self.src_ds.GetGCPs()
            assert gcps, 'Neither geotransform, nor gpcs are in the source file %s' % self.src

            gcp_lst = [(g.Id, g.GCPPixel, g.GCPLine, g.GCPX, g.GCPY, g.GCPZ) for g in gcps]
            ld('src_proj', self.src_ds.GetProjection(), 'gcp_proj', self.src_ds.GetGCPProjection())
            gcp_proj = txt2proj4(self.src_ds.GetGCPProjection())
            if src_proj and gcp_proj != src_proj:
                coords = GdalTransformer(SRC_SRS=gcp_proj, DST_SRS=src_proj).transform([g[3:6] for g in gcp_lst])
                gcp_lst = [tuple(p[:3]+c) for p, c in zip(gcp_lst, coords)]

            gcp_txt = '\n'.join((gcp_templ % g for g in gcp_lst))
            #src_transform = warp_src_gcp_transformer % (0, gcp_txt)
            src_transform = warp_src_tps_transformer % gcp_txt