import os
import numpy as np
from osgeo import gdal,osr
import matplotlib.pyplot as plt
from polsartools.utils.utils import conv2d,time_it
# from pyproj import CRS
import os
def read_bin(file):
ds = gdal.Open(file)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
return arr
def norm_data(data):
data = 10*np.log10(data)
data[data==-np.inf] = np.nan
data[data==np.inf] = np.nan
data = (data - np.nanmin(data)) / (np.nanmax(data) - np.nanmin(data))
p5 = np.nanpercentile(data, 5)
p95 = np.nanpercentile(data, 95)
data = np.clip(data, p5, p95)
data = (data - p5) / (p95 - p5)
return data
def read_and_normalize(file_path):
"""Reads binary data and normalizes it."""
return norm_data(read_bin(file_path)) if os.path.isfile(file_path) else None
# def create_prj(ref_raster_path, output_png_path):
# """
# Creates a .prj file with WKT from the given EPSG code.
# """
# ds = gdal.Open(ref_raster_path)
# wkt = ds.GetProjection()
# # if not wkt:
# # raise ValueError("Reference raster does not contain spatial projection information.")
# crs = CRS.from_wkt(wkt)
# wkt_text = crs.to_wkt()
# prj_path = os.path.splitext(output_png_path)[0] + ".prj"
# with open(prj_path, "w") as f:
# f.write(wkt_text)
def generate_rgb_png(red, green, blue,georef_file, output_path):
# Get geotransform and projection from reference file
ref = gdal.Open(georef_file)
ref_band = ref.GetRasterBand(1).ReadAsArray()
alpha_channel = np.where(ref_band == 0, 0, 255).astype(np.uint8)
del ref_band
"""Creates an RGBA image and saves it."""
rgb_uint8 = (np.dstack((red, green, blue)) * 255).astype(np.uint8)
# alpha_channel = np.where(np.all(rgb_uint8 == 0, axis=2), 0, 255).astype(np.uint8)
rgba_uint8 = np.dstack((rgb_uint8, alpha_channel))
plt.imsave(output_path, rgba_uint8)
# print(f"Pauli RGB image saved as {output_path}")
fig, ax = plt.subplots()
plt.imshow(rgba_uint8, vmin=0, vmax=255)
ax.axis('off')
plt.savefig(output_path.replace(".png", "_thumb.png"), format='png',
bbox_inches='tight', pad_inches=0, transparent=True)
def generate_rgb_tif(red, green, blue, georef_file, output_path):
driver = gdal.GetDriverByName('GTiff')
height, width = red.shape
dataset = driver.Create(output_path, width, height, 4, gdal.GDT_Byte,options=['COMPRESS=DEFLATE','PREDICTOR=2','ZLEVEL=9', 'BIGTIFF=YES','TILED=YES'])
# Prepare the image bands
red_uint8 = (red * 255).astype(np.uint8)
green_uint8 = (green * 255).astype(np.uint8)
blue_uint8 = (blue * 255).astype(np.uint8)
# alpha = np.where((red_uint8 == 0) & (green_uint8 == 0) & (blue_uint8 == 0), 0, 255).astype(np.uint8)
# Get geotransform and projection from reference file
ref = gdal.Open(georef_file)
dataset.SetGeoTransform(ref.GetGeoTransform())
dataset.SetProjection(ref.GetProjection())
ref_band = ref.GetRasterBand(1).ReadAsArray()
alpha = np.where(ref_band == 0, 0, 255).astype(np.uint8)
del ref_band
dataset.GetRasterBand(1).WriteArray(red_uint8)
dataset.GetRasterBand(2).WriteArray(green_uint8)
dataset.GetRasterBand(3).WriteArray(blue_uint8)
dataset.GetRasterBand(4).WriteArray(alpha)
dataset.FlushCache()
ref = None
# print(f"RGB GeoTIFF saved to {output_path}")
def create_pgw(reference_file, output_png_path):
from osgeo import gdal
import os
ds = gdal.Open(reference_file)
gt = ds.GetGeoTransform()
pgw_values = [
gt[1], # pixel width
gt[2], # rotation (typically 0)
gt[4], # rotation (typically 0)
gt[5], # pixel height (typically negative)
gt[0] + gt[1] / 2, # x of center of top-left pixel
gt[3] + gt[5] / 2 # y of center of top-left pixel
]
pgw_path = os.path.splitext(output_png_path)[0] + ".pgw"
with open(pgw_path, "w") as f:
for val in pgw_values:
f.write(f"{val:.10f}\n")
# print(f"PGW file created at {pgw_path}")
[docs]
@time_it
def pauliRGB(infolder,tif_flag=False):
"""
Generate Pauli RGB image from polarimetric SAR data (C4, C3, T4, T3, or S2 matrices).
This function creates a Pauli decomposition RGB composite image from full polarimetric SAR data.
Examples
--------
>>> # Generate Pauli RGB from a polarimetric folder
>>> pauliRGB("/path/to/polSAR_data")
>>> # Output both PNG and GeoTIFF versions
>>> pauliRGB("/path/to/polSAR_data", tif_flag=True)
Parameters
----------
infolder : str
Path to the folder containing polarimetric matrix files (.bin or .tif).
Supports full- and dual-pol data in formats: C4, C3, T4, T3, or S2.
tif_flag : bool, default=False
If True, generates both PNG and GeoTIFF (.tif) versions of the Pauli RGB image.
Returns
-------
None
Writes output to:
- PauliRGB.png: RGB composite with world file for georeferencing
- PauliRGB.tif: Optional GeoTIFF output if `tif_flag=True`
"""
def find_file(infolder, base_name):
for ext in [".bin", ".tif"]:
path = os.path.join(infolder, f"{base_name}{ext}")
if os.path.isfile(path):
return path
return None
c4_keys = ["C11", "C22", "C33", "C44", "C14_real"]
c4_files = [find_file(infolder, key) for key in c4_keys]
if all(c4_files):
C11, C22, C33, C44 = [read_bin(f) for f in c4_files[:4]]
C14_real = read_bin(c4_files[4])
red = norm_data(C11 + C44 - 2 * C14_real)
green = norm_data(C22 + C33)
blue = norm_data(C11 + C44 + 2 * C14_real)
georef_file = c4_files[0]
# C3 fallback
elif all(find_file(infolder, f"C{i}") for i in [11, 22, 33]):
C11 = read_bin(find_file(infolder, "C11"))
C22 = read_bin(find_file(infolder, "C22"))
C33 = read_bin(find_file(infolder, "C33"))
C13_real = read_bin(find_file(infolder, "C13_real"))
blue = norm_data(0.5 * (C11 + C33 + 2 * C13_real))
red = norm_data(0.5 * (C11 + C33 - 2 * C13_real))
green = norm_data(C22)
georef_file = find_file(infolder, "C11")
# T3 fallback
elif all(find_file(infolder, f"T{i}") for i in [11, 22, 33]):
red = read_and_normalize(find_file(infolder, "T22"))
green = read_and_normalize(find_file(infolder, "T33"))
blue = read_and_normalize(find_file(infolder, "T11"))
georef_file = find_file(infolder, "T11")
# S2 fallback
elif all(find_file(infolder, f"s{i}") for i in [11, 12, 22]):
S11 = read_bin(find_file(infolder, "s11"))
S12 = read_bin(find_file(infolder, "s12"))
S22 = read_bin(find_file(infolder, "s22"))
blue = norm_data(np.abs((1 / np.sqrt(2)) * (S11 + S22))**2)
red = norm_data(np.abs((1 / np.sqrt(2)) * (S11 - S22))**2)
green = norm_data(np.abs((2 / np.sqrt(2)) * S12)**2)
georef_file = find_file(infolder, "s11")
else:
raise ValueError("No matching dataset found for C4, C3, T3, or S2!")
output_path = os.path.join(infolder, "PauliRGB.png")
generate_rgb_png(red, green, blue,georef_file, output_path)
create_pgw(georef_file, output_path)
# create_prj(georef_file, output_path)
print(f"Pauli RGB image saved as {output_path}")
if tif_flag:
output_path = os.path.join(infolder, "PauliRGB.tif")
generate_rgb_tif(red, green, blue, georef_file, output_path)
print(f"Pauli RGB image saved as {output_path}")
