edit: even more specific, the problem lies with .astype(uint16). Without this conversion, it loads fine. Also show(limg.read()) also works
original:
The problem lies with the show(rio.open(georeferenced.tif)). The data is saved properly, and when loaded all operations you would want to perform should work fine. if you're only interested in showing the image, you can show(limg.read([1,2,3])
Further, an (ugly) workaround would be to re-load and save the data. Then you can also resample so your rasters overlap (for calculations). When performed directly on the PNG, this crashed my kernel (sounds like a bug).
I also changed the arguments for the rasterio.open('w') to use the shape of your img in stead of the reference image in case they do not exactly match. optionally, you could maybe resample your georeferenced image, so
import rasterio as rio
from rasterio.enums import Resampling
from rasterio.plot import show
#Input png image, to convert as geotiff
img = rio.open('/path/to.png')
img = img.read([1,2,3])
img = img.astype('uint16')
show(img) #shows true color
# Input image for coordinate reference
with rio.open('/path/to/reference.tif') as naip:
#open georeferenced.tif for writing
with rio.open(
'georeferenced.tif',
'w',
driver='GTiff',
count=img.shape[0],
height=img.shape[1],
width=img.shape[2],
dtype=img.dtype,
crs=naip.crs,
transform=naip.transform,
) as dst:
dst.write(img)
with rio.open('georeferenced.tif') as limg:
show(limg) #1 band only shows
show(limg.read([1,2,3])) #shows true color
#resample so pixels overlap with reference
limg = limg.read(out_shape=(3,naip.shape[0],naip.shape[1]),
resampling=Resampling.nearest)
with rio.open('resampled.tif','w',
driver='GTiff',
count=limg.shape[0],
height=limg.shape[1],
width=limg.shape[2],
dtype=limg.dtype,
crs=naip.crs,
transform=naip.transform,
) as dst:
dst.write(limg)
The writing is a bit verbose, but a passing of **naip.profile does not work, since the count is different, unless your reference layer has the same number of bands. an alternative would be storing naip.profile, and then changing the count manually.
Best Answer
Perhaps the easiest way is to use the GDAL tools, which underpin many GIS and mapping systems out there. If you're dealing with lots of files, it is a fairly easy task to write a batch file or script to automate the process.
I'll make some assumptions here about your level of comfort with the command line, but hopefully it should be easy to follow.
First, obtain GDAL. For most Linux distros, there should be a package in your repository. For Windows, it is easiest to download either FWTools or OSGeo4W.
Next, you need to run the program
gdal_translatewhich will do the actual hard work for you. For your needs, it is a fairly simple set of options:This just tells GDAL to output a JPEG image given your source file and a destination filename. It will create two files, the JPEG image itself, and a .jpw (or possibly .wld) file which contains the georeferencing information, without which it would just be a picture.