Tengo un raster (USGS DEM en realidad) y tengo que dividirlo en trozos más pequeños como la imagen de abajo muestra. Eso se logró en ArcGIS 10.0 utilizando la herramienta Split Raster. Me gustaría un método FOSS para hacer esto. He mirado en GDAL, pensando que seguramente lo haría (de alguna manera con gdal_translate), pero no puedo encontrar nada. En última instancia, me gustaría ser capaz de tomar la trama y decir qué tan grande (4KM por 4KM trozos) me gustaría que se divide en.
Respuesta
¿Demasiados anuncios?
Pranjal Kakati
Puntos
1
import gdal
import os
def create_tiles(tile_size, input_filename, in_path='/media/Data/avinash/input/',
out_path='/home/nesac/PycharmProjects/GIS_Data_Automation/data/tiles/'):
ds = gdal.Open(in_path + input_filename)
band = ds.GetRasterBand(1)
output_filename = 'tile_'
xsize, ysize = (band.XSize, band.YSize)
tile_size_x, tile_size_y = tile_size
tile_list = {}
complete_x = xsize // tile_size_x
complete_y = ysize // tile_size_y
residue_x = xsize % tile_size_x
residue_y = ysize % tile_size_y
# for part A
for i in range(complete_x):
for j in range(complete_y):
Xmin = i * tile_size_x
Xmax = i * tile_size_x + tile_size_x - 1
Ymin = j * tile_size_y
Ymax = j * tile_size_y + tile_size_y - 1
# do patch creation here
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
tile_size_x) + ", " + str(tile_size_y) + " " + \
str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)
x_residue_count = 1
y_residue_count = 1
# for part B
for j in range(complete_y):
Xmin = tile_size_x * complete_x
Xmax = tile_size_x * complete_x + residue_x - 1
Ymin = j * tile_size_y
Ymax = j * tile_size_y + tile_size_y - 1
# do patch creation here
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
residue_x) + ", " + str(tile_size_y) + " " + \
str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)
# for part C
for i in range(complete_x):
Xmin = i * tile_size_x
Xmax = i * tile_size_x + tile_size_x - 1
Ymin = tile_size_y * complete_y
Ymax = tile_size_y * complete_y + residue_y - 1
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
tile_size_x) + ", " + str(residue_y) + " " + \
str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)
# for part D
Xmin = complete_x * tile_size_x
Ymin = complete_y * tile_size_y
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
residue_x) + ", " + str(residue_y) + " " + \
str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)
- Ver respuestas anteriores
- Ver más respuestas
0 votos
Tengo una utilidad que utiliza subprocess.Popen para ejecutar múltiples gdal traduce al mismo tiempo que yo uso para la extracción de una gran trama a los azulejos utilizando una red de pesca, particularmente útil si la entrada y / o salida es altamente comprimido (por ejemplo, LZW o deflate GeoTiff), si ninguno es altamente comprimido el proceso de picos en el acceso al disco duro y no es mucho más rápido que correr uno a la vez. Por desgracia, no es lo suficientemente genérico como para compartirlo debido a la rigidez de las convenciones de nomenclatura, pero de todos modos es algo para reflexionar. La opción -multi para GDALWarp a menudo causa problemas y sólo utiliza 2 hilos (uno de lectura, uno de escritura) no todos disponibles.