2

I'm clipping a raster to a shapefile, which returned a masked array, and my goal is to write that masked array to a new .tif file with the correct metadata. I can view the masked array and it looks great. Then, when I write it to a new file using rasterio's write function, something goes wrong, because when I load it back in, it's the wrong shape and it won't plot. See below for details and images:

  1. The original raster and polygon shapefile watershed boundary that I want to clip it to.

enter image description here

  1. Clipped output array and image:
out_img, out_transform = mask(rastDat, shapes=coords, crop=True,indexes=1)

out_mask = np.ma.masked_where(out_img<0,out_img)

fig, ax = plt.subplots(figsize=(8,12))
outPlot = ax.imshow(out_mask,cmap='terrain')

enter image description here

  1. Looks great at this point, and I write it to a new file
#copying metadata from original raster
out_meta = rastDat.meta.copy()

#amending original metadata
out_meta.update({'nodata': -99,
                 'height' : out_img.shape[0],
                 'width' : out_img.shape[1],
                 'transform' : out_transform}) 

#writing and then re-reading the output data to see if it looks good

with rio.open(outLoc,'w',**out_meta) as dst:
    dst.write(out_fill,1)

newClip = rio.open(outLoc)
clipIM = newClip.read(1,masked=True)

  1. however, the new data is not plotting and has a weird shape. Am I doing something wrong on the writing, or reading, or both? How can I write the raster so that the data will come in and plot like the first figure in this post?

enter image description here

  1. To field any questions, I have tried writing as dst.write(out_fill) as well as reading the data in with and without the mask.

closed as off-topic by user2856, Vince, Jochen Schwarze, nmtoken, MrXsquared Oct 4 at 13:08

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This problem cannot or can no longer be reproduced. Changes to the system or to the asker's circumstances may have rendered the question obsolete, or the question does not include a procedure to enable potential answerers to reproduce the same symptoms. Such questions are off-topic as they are unlikely to help future readers, but editing them to include more details can lead to re-opening." – user2856, Vince, Jochen Schwarze, nmtoken, MrXsquared
If this question can be reworded to fit the rules in the help center, please edit the question.

  • Your code works for me (on my data), though I don't plot it, I just checked outLoc displays correctly in QGIS and that clipIM.shape was correct. Does outLoc display correctly if you open it in QGIS (if you have it installed)? – user2856 Oct 3 at 21:25
  • Well, I don't have QGIS, but when I load it in an entirely new .py (or ipynb) file, it works! Wow, silly. Thanks so much!!! =D – swimdonovan Oct 3 at 22:01
0

You may be prescribing the new height and width incorrectly. What happens when you call out_image.shape before writing the raster? It may be the shape is actually 1 x 6094 x 4000 (or whatever the width is) and so when you prescribe shape[0] as the height you are actually giving it a height of 1. The value of 1 in the zeroth element of the shape tuple is indicating this is a single-band raster (rather than a three-band RGB image, for example).

Try changing the elements of the shape tuple you are assigning to the height and width keys of the metadata:

# amending original metadata out_meta.update({'nodata': -99, 'height' : out_img.shape[1], 'width' : out_img.shape[2], 'transform' : out_transform})

  • Hi Charlie, thanks for that suggestion. I just looked into it to make sure, but it seems that it's correct when I use shape[0] and shape[1]. if I call: out_fill.shape I get: (6094, 4172) out_meta` And all looks well {'driver': 'GTiff', 'dtype': 'int16', 'nodata': -99, 'width': 4172, 'height': 6094, 'count': 1, 'crs': CRS.from_dict(init='epsg:27200'), 'transform': Affine(26.683493831244363, 0.0, 2029176.4496368268, 0.0, -26.683493831244363, 5591892.947473064)} – swimdonovan Oct 3 at 16:55
  • how are you creating out_fill ? I don't see it in the code in your example. – Charlie Parr Oct 3 at 17:04
  • Ohh, sorry. That's just a temporary copy i made of out_mask, which is plotted in (2) of the original post. I had made a copy when I was playing around and called it 'fill'. – swimdonovan Oct 3 at 17:22
  • hmm well i am stumped for now - clearly the raster you are writing as the wrong shape (1 x 6094) so something is off. share the data you are working with if possible and maybe myself or someone else can dig deeper – Charlie Parr Oct 3 at 17:31
  • It's uploaded to a dropbox! see link here: dropbox.com/sh/zq9s6fx4ol9quza/AAAqmcs69-BFkMCbJSvrlXn_a?dl=0 Just change the directory on the input files and you're good to go! Thanks in advance for helping me out!!! – swimdonovan Oct 3 at 17:45

Not the answer you're looking for? Browse other questions tagged or ask your own question.