Skip to content

Commit

Permalink
add get_extent to cog mod
Browse files Browse the repository at this point in the history
  • Loading branch information
sharkAndshark committed Feb 2, 2025
1 parent 31b781c commit 1be025c
Show file tree
Hide file tree
Showing 2 changed files with 128 additions and 1 deletion.
12 changes: 12 additions & 0 deletions martin/src/cog/errors.rs
Original file line number Diff line number Diff line change
Expand Up @@ -42,4 +42,16 @@ pub enum CogError {

#[error("Striped tiff file is not supported, the tiff file is {0}")]
NotSupportedChunkType(PathBuf),

#[error("No geospatial information found in {0}. Either ModelTransformationTag or both ModelTiepointTag and ModelPixelScaleTag are required")]
MissingGeospatialInfo(PathBuf),

#[error("Invalid ModelTiepointTag: expected at least 6 values, got {0}")]
InvalidModelTiepoint(usize),

#[error("Invalid ModelPixelScaleTag: expected at least 3 values, got {0}")]
InvalidModelPixelScale(usize),

#[error("Invalid ModelTransformationTag: expected at least 12 values for 3x4 matrix, got {0}")]
InvalidModelTransformation(usize),
}
117 changes: 116 additions & 1 deletion martin/src/cog/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,20 @@ impl CogSource {
pub fn new(id: String, path: PathBuf) -> FileResult<Self> {
let tileinfo = TileInfo::new(Format::Png, martin_tile_utils::Encoding::Uncompressed);
let meta = get_meta(&path)?;
let tilejson = tilejson! {
let mut tilejson = tilejson! {
tiles: vec![],
minzoom: meta.min_zoom,
maxzoom: meta.max_zoom
};

// Add COG-specific metadata to tilejson.other
let mut cog_info = serde_json::Map::new();
// Add extent bounds
cog_info.insert("extent".to_string(), serde_json::json!(meta.extent));
// Add all COG info to tilejson.other
tilejson
.other
.insert("cogMetadata".to_string(), serde_json::json!(cog_info));
Ok(CogSource {
id,
path,
Expand All @@ -60,6 +69,7 @@ struct Meta {
max_zoom: u8,
zoom_and_ifd: HashMap<u8, usize>,
zoom_and_tile_across_down: HashMap<u8, (u32, u32)>,
extent: [f64; 4], // [minx, miny, maxx, maxy] bounds
nodata: Option<f64>,
}

Expand Down Expand Up @@ -308,6 +318,28 @@ fn get_meta(path: &PathBuf) -> Result<Meta, FileError> {
None
};

let model_transformation = decoder.get_tag_f64_vec(Tag::ModelTransformationTag).ok();
let model_tiepoint = decoder.get_tag_f64_vec(Tag::ModelTiepointTag).ok();
let pixel_scale = decoder.get_tag_f64_vec(Tag::ModelPixelScaleTag).ok();

let (width, height) = decoder.dimensions().map_err(|e| {
CogError::TagsNotFound(
e,
vec![Tag::ImageWidth.to_u16(), Tag::ImageLength.to_u16()],
0,
path.to_path_buf(),
)
})?;

let extent = get_extent(
model_transformation,
model_tiepoint,
pixel_scale,
width,
height,
path.clone(),
)?;

let images_ifd = get_images_ifd(&mut decoder, path);

for (idx, image_ifd) in images_ifd.iter().enumerate() {
Expand All @@ -333,6 +365,7 @@ fn get_meta(path: &PathBuf) -> Result<Meta, FileError> {
max_zoom: images_ifd.len() as u8 - 1,
zoom_and_ifd,
zoom_and_tile_across_down,
extent,
nodata,
})
}
Expand All @@ -350,6 +383,88 @@ fn get_grid_dims(
Ok((tiles_across, tiles_down))
}

/// Calculate the extent [minx, miny, maxx, maxy] of a GeoTIFF image
fn get_extent(
model_transformation: Option<Vec<f64>>,
model_tiepoint: Option<Vec<f64>>,
pixel_scale: Option<Vec<f64>>,
width: u32,
height: u32,
path: PathBuf,
) -> Result<[f64; 4], CogError> {
match (model_transformation, model_tiepoint, pixel_scale) {
(Some(transform), _, _) => get_extent_from_transform(&transform, width, height),
(None, Some(tiepoint), Some(pixel_scale)) => {
get_extent_from_tiepoint(&tiepoint, &pixel_scale, width, height)
}
_ => Err(CogError::MissingGeospatialInfo(path)),
}
}

/// Calculate the extent [minx, miny, maxx, maxy] using model transformation matrix
fn get_extent_from_transform(
transform: &[f64],
width: u32,
height: u32,
) -> Result<[f64; 4], CogError> {
// ModelTransformationTag should have at least 12 values for a 3x4 matrix
if transform.len() < 12 {
return Err(CogError::InvalidModelTransformation(transform.len()));
}

let corners = [
(0.0, 0.0),
(0.0, height as f64),
(width as f64, 0.0),
(width as f64, height as f64),
];

let mut xs = Vec::with_capacity(4);
let mut ys = Vec::with_capacity(4);

// Apply transformation matrix to each corner
for (i, j) in corners {
let x = transform[3] + (transform[0] * i) + (transform[1] * j);
let y = transform[7] + (transform[4] * i) + (transform[5] * j);
xs.push(x);
ys.push(y);
}

Ok([
*xs.iter().min_by(|a, b| a.total_cmp(b)).unwrap(),
*ys.iter().min_by(|a, b| a.total_cmp(b)).unwrap(),
*xs.iter().max_by(|a, b| a.total_cmp(b)).unwrap(),
*ys.iter().max_by(|a, b| a.total_cmp(b)).unwrap(),
])
}

/// Calculate the extent [minx, miny, maxx, maxy] using model tiepoint and pixel scale
fn get_extent_from_tiepoint(
tiepoint: &[f64],
pixel_scale: &[f64],
width: u32,
height: u32,
) -> Result<[f64; 4], CogError> {
// ModelTiepointTag should have at least 6 values (I,J,K,X,Y,Z)
if tiepoint.len() < 6 {
return Err(CogError::InvalidModelTiepoint(tiepoint.len()));
}

// ModelPixelScaleTag should have 3 values (ScaleX, ScaleY, ScaleZ)
if pixel_scale.len() < 3 {
return Err(CogError::InvalidModelPixelScale(pixel_scale.len()));
}

let x1 = tiepoint[3]; // Origin X
let y1 = tiepoint[4]; // Origin Y

// Calculate max extent using resolution/pixel scale
let x2 = x1 + (pixel_scale[0] * width as f64);
let y2 = y1 + (pixel_scale[1] * height as f64);

Ok([x1.min(x2), y1.min(y2), x1.max(x2), y1.max(y2)])
}

fn get_image_dims(
decoder: &mut Decoder<File>,
path: &Path,
Expand Down

0 comments on commit 1be025c

Please sign in to comment.