Python gdal.ReprojectImage方法代码示例

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import ReprojectImage [as 别名]
def reproject_to_grid_coordinates(self, grid_coordinates, interp=gdalconst.GRA_NearestNeighbour):
        """ Reprojects data in this layer to match that in the GridCoordinates
        object. """
        source_dataset = self.grid_coordinates._as_gdal_dataset()
        dest_dataset = grid_coordinates._as_gdal_dataset()
        rb = source_dataset.GetRasterBand(1)
        rb.SetNoDataValue(NO_DATA_VALUE)
        rb.WriteArray(np.ma.filled(self.raster_data, NO_DATA_VALUE))

        gdal.ReprojectImage(source_dataset, dest_dataset,
                            source_dataset.GetProjection(),
                            dest_dataset.GetProjection(),
                            interp)
        dest_layer = self.clone_traits()
        dest_layer.grid_coordinates = grid_coordinates
        rb = dest_dataset.GetRasterBand(1)
        dest_layer.raster_data = np.ma.masked_values(rb.ReadAsArray(), NO_DATA_VALUE)
        return dest_layer 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import ReprojectImage [as 别名]
def scale_query_to_tile(self, dsquery, dstile):
        """Scales down query dataset to the tile dataset"""

        querysize = dsquery.RasterXSize
        tilesize = dstile.RasterXSize
        tilebands = dstile.RasterCount

        if self.options.resampling == 'average':
            for i in range(1,tilebands+1):
                res = gdal.RegenerateOverview( dsquery.GetRasterBand(i),
                    dstile.GetRasterBand(i), 'average' )
                if res != 0:
                    self.error("RegenerateOverview() failed")
        else:
            # Other algorithms are implemented by gdal.ReprojectImage().
            dsquery.SetGeoTransform( (0.0, tilesize / float(querysize), 0.0, 0.0, 0.0, tilesize / float(querysize)) )
            dstile.SetGeoTransform( (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) )
            
            res = gdal.ReprojectImage(dsquery, dstile, None, None, self.resampling)    
            if res != 0:
                self.error("ReprojectImage() failed on %s, error %d" % (tilefilename, res)) 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import ReprojectImage [as 别名]
def reproject_image(in_raster, out_raster_path, new_projection,  driver = "GTiff",  memory = 2e3, do_post_resample=True):
    """
    Creates a new, reprojected image from in_raster using the gdal.ReprojectImage function.

    Parameters
    ----------
    in_raster
        Either a gdal.Raster object or a path to a raster
    out_raster_path
        The path to the new output raster.
    new_projection
        The new projection in .wkt
    driver
        The format of the output raster.
    memory
        The amount of memory to give to the reprojection.
    do_post_resample
        If set to false, do not resample the image back to the original projection.

    Notes
    -----
    The GDAL reprojection routine changes the size of the pixels by a very small amount; for example, a 10m pixel image
    can become a 10.002m pixel resolution image. To stop alignment issues,
    by deafult this function resamples the images back to their original resolution

    """
    if type(new_projection) is int:
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(new_projection)
        new_projection = proj.ExportToWkt()
    log = logging.getLogger(__name__)
    log.info("Reprojecting {} to {}".format(in_raster, new_projection))
    if type(in_raster) is str:
        in_raster = gdal.Open(in_raster)
    res = in_raster.GetGeoTransform()[1]
    gdal.Warp(out_raster_path, in_raster, dstSRS=new_projection, warpMemoryLimit=memory, format=driver)
    # After warping, image has irregular gt; resample back to previous pixel size
    # TODO: Make this an option
    if do_post_resample:
        resample_image_in_place(out_raster_path, res)
    return out_raster_path 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import ReprojectImage [as 别名]
def create_display_layer(class_path, out_path, class_color_key):
    srs = """PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]"""
    class_raster = gdal.Open(class_path)
    display_raster = pyeo.raster_manipulation.create_matching_dataset(class_raster, out_path, bands=3, datatype=gdal.GDT_Byte)
    display_array = display_raster.GetVirtualMemArray(eAccess=gdal.GF_Write)

    class_array = class_raster.GetVirtualMemArray()
    for index, class_pixel in np.ndenumerate(class_array):
        display_array[:, index[0], index[1]] =\
            [class_row[1:4] for class_row in class_color_key if class_row[0] == str(class_pixel)][0]
    display_array = None
    class_array = None
    gdal.ReprojectImage(display_raster, dst_wkt=srs)
    display_raster = None
    class_raster = None 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import ReprojectImage [as 别名]
def reproject_dataset(dataset, output_file_name, wkt_from="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs", epsg_to=32643, pixel_spacing=1000, file_format='GTiff'):
    '''
    :param dataset: Modis subset name use gdal.GetSubdatasets()
    :param output_file_name: file location along with proper extension
    :param wkt_from: Modis wkt (default)
    :param epsg_to: integer default(32643)
    :param pixel_spacing: in metres
    :param file_format: default GTiff
    :return:
    '''
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(epsg_to)
    modis = osr.SpatialReference()
    modis.ImportFromProj4(wkt_from)
    tx = osr.CoordinateTransformation(modis, wgs84)
    g = gdal.Open(dataset)
    geo_t = g.GetGeoTransform()
    print geo_t
    x_size = g.RasterXSize
    y_size = g.RasterYSize
    (ulx, uly, ulz) = tx.TransformPoint(geo_t[0], geo_t[3])
    (lrx, lry, lrz) = tx.TransformPoint(geo_t[0] + (geo_t[1]*x_size), geo_t[3]+ (geo_t[5]*y_size))
    mem_drv = gdal.GetDriverByName(file_format)
    dest = mem_drv.Create(output_file_name, int((lrx-ulx)/pixel_spacing), int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)
    new_geo = ([ulx, pixel_spacing, geo_t[2], uly, geo_t[4], -pixel_spacing])
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(wgs84.ExportToWkt())
    gdal.ReprojectImage(g, dest, modis.ExportToWkt(), wgs84.ExportToWkt(), gdal.GRA_Bilinear)
    print "reprojected" 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import ReprojectImage [as 别名]
def Raster_to_Array(input_tiff, ll_corner, x_ncells, y_ncells,
                    values_type='float32'):
    """
    Loads a raster into a numpy array
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_data_type = inp_band.DataType

    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    ll_x = ll_corner[0]
    ll_y = ll_corner[1]

    top_left_x = ll_x
    top_left_y = ll_y - cellsize_y*y_ncells

    # Change start point
    temp_path = tempfile.mkdtemp()
    temp_driver = gdal.GetDriverByName('GTiff')
    temp_tiff = os.path.join(temp_path, os.path.basename(input_tiff))
    temp_source = temp_driver.Create(temp_tiff, x_ncells, y_ncells,
                                     1, inp_data_type)
    temp_source.GetRasterBand(1).SetNoDataValue(NoData_value)
    temp_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                 top_left_y, rot_2, cellsize_y))
    temp_source.SetProjection(inp_srs)

    # Snap
    gdal.ReprojectImage(inp_lyr, temp_source, inp_srs, inp_srs,
                        gdal.GRA_Bilinear)
    temp_source = None

    # Read array
    d_type = pd.np.dtype(values_type)
    out_lyr = gdal.Open(temp_tiff)
    array = out_lyr.ReadAsArray(0, 0, out_lyr.RasterXSize,
                                out_lyr.RasterYSize).astype(d_type)
    array[pd.np.isclose(array, NoData_value)] = pd.np.nan
    out_lyr = None

    return array 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import ReprojectImage [as 别名]
def scale_query_to_tile(self, dsquery, dstile, tilefilename=''):
        """Scales down query dataset to the tile dataset"""

        querysize = dsquery.RasterXSize
        tilesize = dstile.RasterXSize
        tilebands = dstile.RasterCount

        if self.options.resampling == 'average':

            # Function: gdal.RegenerateOverview()
            for i in range(1, tilebands + 1):
                # Black border around NODATA
                #if i != 4:
                #   dsquery.GetRasterBand(i).SetNoDataValue(0)
                res = gdal.RegenerateOverview(
                    dsquery.GetRasterBand(i), dstile.GetRasterBand(i),
                    'average')
                if res != 0:
                    self.error("RegenerateOverview() failed on %s, error %d" %
                               (tilefilename, res))

        elif self.options.resampling == 'antialias':

            # Scaling by PIL (Python Imaging Library) - improved Lanczos
            array = numpy.zeros((querysize, querysize, tilebands), numpy.uint8)
            for i in range(tilebands):
                array[:, :, i] = gdalarray.BandReadAsArray(
                    dsquery.GetRasterBand(i + 1), 0, 0, querysize, querysize)
            im = Image.fromarray(array, 'RGBA')  # Always four bands
            im1 = im.resize((tilesize, tilesize), Image.ANTIALIAS)
            if path.exists(tilefilename):
                im0 = Image.open(tilefilename)
                im1 = Image.composite(im1, im0, im1)
            im1.save(tilefilename, self.tiledriver)

        else:

            # Other algorithms are implemented by gdal.ReprojectImage().
            dsquery.SetGeoTransform((0.0, tilesize / float(querysize), 0.0,
                                     0.0, 0.0, tilesize / float(querysize)))
            dstile.SetGeoTransform((0.0, 1.0, 0.0, 0.0, 0.0, 1.0))
            
            
            res = gdal.ReprojectImage(dsquery, dstile, None, None,
                                      self.resampling)
            if res != 0:
                self.error("ReprojectImage() failed on %s, error %d" %
                           (tilefilename, res))

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