First time posting here. I'm working on my first remote sensing project to estimate forest structure attributes (e.g. forest height) using satellite images with deep learning. I've struggled for a while constructing the datasets for learning. Specifially, making sure the reference forest maps overlap the satellite images. QGIS actually shows the correct overlap after preprocessing (see picture 1). However, when I load the image in Python using rasterio each index (pixel) does not correspond the same pixel in the other images (see picture 2).
Why does QGIS show nearly perfect alignment while Python show a clear misalignment? Is there something wrong in the way I've preprocessed the image?
QGIS
Reference map (white to red gradient) overlaps background optical image nearly perfectly.
Python
The same images (left: optical, mid: sar, right: reference map) but accessing by the same indices (pixels) show different areas (by up to 50 pixels)
Preprocessing steps
- Reprojecting - EPSG:4326 -> EPSG 32632
gdalwarp -t_srs EPSG:32632 -r lanczos -wo SOURCE_EXTRA=1000 -co COMPRESS=LZW opt.tif opt_reprojected.tif
- Clipping - Clip all images according to a polygon that is an intersection between all images. See
clip_img_to_aoi() - Transforming - Manually fix the last misalignments between imgs
gdalwarp -ct "+proj=affine +xoff=16 +yoff=0" opt.tif transformed_opt.tif
- Resizing - Convert image to correct resolution and size (i.e. numpy shape). See
resize_img()
Post-processing Rasterio Stats
Ground Truth (Reference map)
(1, 8904, 11801)
EPSG:32632
Affine:
| 0.50, 0.00, 290779.00|
| 0.00,-0.50, 6809696.50|
| 0.00, 0.00, 1.00|Opt
(4, 8904, 11801)
EPSG:32632
Affine:
| 0.50, 0.00, 290778.15|
| 0.00,-0.50, 6809677.65|
| 0.00, 0.00, 1.00|
Clip image to area of interest
def clip_img_to_aoi(img_path, aoi_path='data/aoi.geojson'):
with rasterio.open(img_path) as img_ds:
with fiona.open(aoi_path) as geo:
aoi = [feature["geometry"] for feature in geo]
img_clipped, img_transform = rasterio.mask.mask(img_ds, aoi, crop=True)
meta = img_ds.meta.copy()
meta.update({"transform": img_transform,
"height":img_clipped.shape[1],
"width":img_clipped.shape[2]})
with rasterio.open(img_ds.name[:-4] + '_clipped.tif', 'w', **meta) as f:
f.write(img_clipped)
Resize image to appropriate resolution and shape
def resize_img(SAT_img, from_res, to_res):
factor = to_res / from_res
shape = (np.round(SAT_img.shape[1] / factor).astype(int) , np.round(SAT_img.shape[2] / factor).astype(int))
SAT_image_resized = np.zeros((SAT_img.shape[0], shape[0], shape[1]))
for channel in range(0, SAT_img.shape[0]):
print("Resizing channel ", channel)
SAT_image_resized[channel, :, :] = resize(SAT_img[channel, :, :], (shape[0], shape[1]), anti_aliasing=True)
return SAT_image_resized
def change_meta_data(meta_data, factor):
meta_dataNew = meta_data
meta_dataT = meta_dataNew['transform']
meta_dataNew['transform'] = Affine(meta_dataT[0] * factor, meta_dataT[1], meta_dataT[2], \
meta_dataT[3], meta_dataT[4] * factor, meta_dataT[5])
meta_dataNew['width'] = np.round(meta_data['width'] / factor).astype(int)
meta_dataNew['height'] = np.round(meta_data['height'] / factor).astype(int)
return meta_dataNew
