Saya mencoba membuat raster menggunakan pustaka python gdal dan saya sudah sampai pada titik di mana data menjadi output, tetapi data output dibalik pada sumbu x titik asal . Saya tahu saya pasti mengabaikan sesuatu, tetapi saya tidak tahu di mana saya salah. Ada ide?
Saat membuat raster, saya mengatur nilai x / y kiri atas, dan array tampaknya diindeks dari kiri atas dan terus turun ke kanan bawah. Dalam kode di bawah ini saya mengisi array dengan nilai baris.
Saat mencetak array terlihat seperti ini:
[[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1.]
[ 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
2. 2. 2. 2. 2. 2.]
[ 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3.
3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3.
3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3.
3. 3. 3. 3. 3. 3.]
[ 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
4. 4. 4. 4. 4. 4.]
...
Dan data ini berhasil ditulis ke band raster. Namun ketika dilihat di MapWindow GIS , data tampak berlawanan arah dengan titik asal kiri-atas yang awalnya ditetapkan , muncul sebagai nilai kiri-bawah .
Dengan kata lain, data dibalik pada sumbu x titik asal .
import gdal
import osr
import numpy
OUTPUT_FORMAT = "GTiff"
def create_raster(filename="test.tif"):
driver = gdal.GetDriverByName(OUTPUT_FORMAT)
band_type = gdal.GDT_Byte
number_of_bands = 1
x_rotation = 0 # not supported
y_rotation = 0 # not supported
cell_width_meters = 50
cell_height_meters = 50
(min_lon, min_lat, max_lon, max_lat) = _get_point_bounds() # retrieve bounds for point data
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
srs.SetUTM( 54, True) # Set projected coordinate system to handle meters
# create transforms for point conversion
wgs84_coordinate_system = srs.CloneGeogCS() # clone only the geographic coordinate system
wgs84_to_utm_transform = osr.CoordinateTransformation(wgs84_coordinate_system, srs)
# convert to UTM
top_left_x, top_left_y, z = wgs84_to_utm_transform.TransformPoint(min_lon, max_lat, 0)
lower_right_x, lower_right_y, z = wgs84_to_utm_transform.TransformPoint(max_lon, min_lat, 0)
cols, rows = _get_raster_size(top_left_x, lower_right_y, lower_right_x, top_left_y, cell_width_meters, cell_height_meters)
dataset = driver.Create(filename, cols, rows, number_of_bands, band_type) #
# GeoTransform parameters
# --> need to know the area that will be covered to define the geo tranform
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, cell_height_meters ]
dataset.SetGeoTransform(geo_transform)
dataset.SetProjection(srs.ExportToWkt())
dataset_band = dataset.GetRasterBand(1)
data = dataset_band.ReadAsArray(0, 0, cols, rows).astype(numpy.float32) # returns empty array
for row in xrange(rows):
for col in xrange(cols):
data[row][ col] = row + 1
dataset_band.WriteArray(data, 0, 0)
dataset_band.SetNoDataValue(0.0)
dataset_band.FlushCache()
dataset = None # Close file
Saya juga memperhatikan ketika saya menghitung posisi pixel untuk lat / lon yang diberikan nilai-y menghasilkan indeks negatif, yang tampaknya agak benar mengingat array dari kiri atas ke kanan bawah .
inverse_geo_transform = gdal.InvGeoTransform(self.geo_transform)[1] # for mapping lat/lon to pixel
pixel_x, pixel_y = gdal.ApplyGeoTransform(self.inverse_geo_transform, utm_x, utm_y)