Python gdal.ReprojectImage() Examples
The following are 8
code examples of gdal.ReprojectImage().
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 also want to check out all available functions/classes of the module
gdal
, or try the search function
.
.
Example #1
| Source File: my_types.py From pydem with Apache License 2.0 | 6 votes |
|
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
Example #2
| Source File: gdal2cesium.py From gdal2cesium with GNU General Public License v2.0 | 6 votes |
|
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))
Example #3
| Source File: raster_manipulation.py From pyeo with GNU General Public License v3.0 | 5 votes |
|
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
Example #4
| Source File: create_eolabs_layers.py From pyeo with GNU General Public License v3.0 | 5 votes |
|
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
Example #5
| Source File: modis_tg_halli.py From hydrology with GNU General Public License v3.0 | 5 votes |
|
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"
Example #6
| Source File: functions.py From hants with Apache License 2.0 | 4 votes |
|
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
Example #7
| Source File: functions.py From wa with Apache License 2.0 | 4 votes |
|
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
Example #8
| Source File: gdal2tiles_parallel.py From geopackage-python with GNU General Public License v3.0 | 4 votes |
|
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))
# -------------------------------------------------------------------------
