-1

I am trying to write some code that takes a collection of png files + some metadata information (not from a geotiff) and creates a geotiff file. That one geotiff file would have a number of bands/layers equal to the number of png files passed as input. When I try to open the file in QGIS, it is unviewable. By comparing my geotiff file to other, established geotiff files, I noticed that the extent field was empty, which seems wrong (looking at the image information from qgis shown in image below).

enter image description here

I am utilizing the rasterio python library to create the metadata object. By printing out the metadata of my established raster file, I was able to figure out what the profile fields were and what I should set them to (lines 1-20)

I have cannibalized some code from here to add the png files as layers. Even though I am not setting the layers to layers from a different geotiff file, the way the png is formatted looks identical to what a layer should be formatted as. (lines 22-26).

I am confused because when I open my resulting tiff file with a generic png viewer on Windows 10, I can see the general shape of my image the way I would expect, and since I am setting the top-left corner coordinates with the transform field, and the size of the resulting image with the width and height fields, that makes sense to me. However, it would seem that qgis needs the extent information to be specified somehow.

I have scoured rasterio's github and documentation, and am struggling to find any helpful information. My code is given below; any help would be much appreciated!

  1  png = cv2.imread(source_files[0], 0)
  2   
  3  new_width = png.shape[0]
  4  new_height = png.shape[1]
  5  new_count = len(source_files)
  6  new_crs = rio.crs.CRS.from_string("EPSG:32719")
  7  new_transform = Affine(0.0, 0.0, -74.262057, 0, 0, -11.920552)
  8   
  9  new_metadata = {'driver': 'GTiff',
 10                  'dtype': 'float32',
 11                  'nodata': None,
 12                  'width': new_width,
 13                  'height': new_height,
 14                  'count': new_count,
 15                  'crs': new_crs,
 16                  'transform': new_transform,
 17                  'tiled': False,
 18                  'interleave': 'band'}
 19  
 20  metadata = rio.profiles.Profile(data=new_metadata)   
 21
 22  with rio.open('../output/pngs/output.tif', 'w', **metadata) as dst :
 23      for i, layer in enumerate(source_files, start=1):
 24          png = cv2.imread(layer, 0)
 25          png = np.array(png, np.float32) / 255.0
 26          dst.write_band(i, png)

Edit: By checking the bounds field of the established geotiff file and my created geotiff file, I can confirm that this is the issue; the bounding box in my created file has no width or height, strangely. However, it would appear that this field is not writable, so my question still stands.

1 Answers1

0

My issue was in my transformation matrix. I did not want there to be any rotation, so I foolishly set the rotation matrix to be 0. Instead I should have set it to the identity matrix. All is well.

Before: new_transform = Affine(0.0, 0.0, -74.262057, 0, 0, -11.920552)

After: new_transform = Affine(1.0, 0.0, -74.262057, 0.0, 1.0, -11.920552)

  • You have not changed the rotation by making that adjustment; rather you have set the width and height of a pixel to 1.0, rather than 0.0. The rotation parameters (b and d) are still 0.0. https://www.perrygeo.com/python-affine-transforms.html – alphabetasoup Mar 17 '20 at 00:54