Skip to content

Commit

Permalink
Gpkg, GeoJSON のジオメトリ生成コードを汎化する (#154)
Browse files Browse the repository at this point in the history 
nusamai-gpkg と nusamai-geojson において、(vertices+indices)
な表現のジオメトリと、そうでないふつうのジオメトリを、なるべく共通のコードで扱えるようにする。
  • Loading branch information
@ciscorn
ciscorn authored Jan 7, 2024
1 parent 5e12714 commit eb90961
Show file tree
Hide file tree
Showing 5 changed files with 170 additions and 152 deletions.
142 changes: 105 additions & 37 deletions nusamai-geojson/src/conversion.rs
View file
Original file line number Diff line number Diff line change
@@ -1,71 +1,91 @@
use nusamai_geometry::{MultiLineString, MultiPoint, MultiPolygon, Polygon};
use nusamai_geometry::{CoordNum, MultiLineString, MultiPoint, MultiPolygon};

/// Create a GeoJSON geometry from `nusamai_geometry::MultiPolygon`
/// Create a GeoJSON MultiPolygon from `nusamai_geometry::MultiPolygon`.
pub fn multipolygon_to_geometry(mpoly: &MultiPolygon<3>) -> geojson::Geometry {
let rings_list = mpoly
.iter()
.map(|poly| {
poly.rings()
.map(|c| c.iter_closed().map(|v| v.to_vec()).collect())
.collect::<Vec<_>>()
})
.collect();
geojson::Value::MultiPolygon(rings_list).into()
multipolygon_to_geometry_with_mapping(mpoly, |c| c.to_vec())
}

/// Create a GeoJSON geometry from vertices and MultiPolygon indices.
/// Create a GeoJSON MultiPolygon from vertices and indices.
pub fn indexed_multipolygon_to_geometry(
vertices: &[[f64; 3]],
mpoly: &MultiPolygon<1, u32>,
mpoly_idx: &MultiPolygon<1, u32>,
) -> geojson::Geometry {
let rings_list: Vec<geojson::PolygonType> = mpoly
multipolygon_to_geometry_with_mapping(mpoly_idx, |idx| vertices[idx[0] as usize].to_vec())
}

/// Create a GeoJSON MultiPolygon from `nusamai_geometry::MultiPolygon` with a mapping function.
pub fn multipolygon_to_geometry_with_mapping<const D: usize, T: CoordNum>(
mpoly: &MultiPolygon<D, T>,
mapping: impl Fn(&[T]) -> Vec<f64>,
) -> geojson::Geometry {
let coords: Vec<geojson::PolygonType> = mpoly
.iter()
.map(|poly| indexed_polygon_to_rings(vertices, &poly))
.map(|poly| {
poly.rings()
.map(|ls| {
ls.iter_closed()
.map(&mapping) // Get the actual coord values
.collect()
})
.collect::<Vec<_>>()
})
.collect();
geojson::Value::MultiPolygon(rings_list).into()
geojson::Value::MultiPolygon(coords).into()
}

fn indexed_polygon_to_rings(
vertices: &[[f64; 3]],
poly_idx: &Polygon<1, u32>,
) -> geojson::PolygonType {
poly_idx
.rings()
.map(|ls| {
ls.iter_closed()
.map(|idx| vertices[idx[0] as usize].to_vec()) // Get the actual coord values
.collect()
})
.collect()
/// Create a GeoJSON MultiLineString from `nusamai_geometry::MultiLineString`.
pub fn multilinestring_to_geometry(mls: &MultiLineString<3>) -> geojson::Geometry {
multilinestring_to_geometry_with_mapping(mls, |c| c.to_vec())
}

/// Create a GeoJSON geometry from vertices and MultiLineString indices.
/// Create a GeoJSON MultiLineString from vertices and indices.
pub fn indexed_multilinestring_to_geometry(
vertices: &[[f64; 3]],
mls_idx: &MultiLineString<1, u32>,
) -> geojson::Geometry {
let mls_coords = mls_idx
multilinestring_to_geometry_with_mapping(mls_idx, |idx| vertices[idx[0] as usize].to_vec())
}

/// Create a GeoJSON MultiLineString from `nusamai_geometry::MultiPolygon` with a mapping function.
pub fn multilinestring_to_geometry_with_mapping<const D: usize, T: CoordNum>(
mls: &MultiLineString<D, T>,
mapping: impl Fn(&[T]) -> Vec<f64>,
) -> geojson::Geometry {
let coords = mls
.iter()
.map(|ls_idx| {
ls_idx
.iter()
.map(|idx| vertices[idx[0] as usize].to_vec()) // Get the actual coord values
.map(&mapping) // Get the actual coord values
.collect()
})
.collect();
geojson::Value::MultiLineString(mls_coords).into()
geojson::Value::MultiLineString(coords).into()
}

/// Create a GeoJSON MultiPoint from `nusamai_geometry::MultiPoint`.
pub fn multipoint_to_geometry(mpoint: &MultiPoint<3>) -> geojson::Geometry {
multipoint_to_geometry_with_mapping(mpoint, |c| c.to_vec())
}

/// Create a GeoJSON geometry from vertices and MultiPoint indices.
/// Create a GeoJSON MultiPoint from vertices and indices.
pub fn indexed_multipoint_to_geometry(
vertices: &[[f64; 3]],
mpoint: &MultiPoint<1, u32>,
mpoint_idx: &MultiPoint<1, u32>,
) -> geojson::Geometry {
multipoint_to_geometry_with_mapping(mpoint_idx, |idx| vertices[idx[0] as usize].to_vec())
}

/// Create a GeoJSON MultiPoint from `nusamai_geometry::MultiPoint` with a mapping function.
pub fn multipoint_to_geometry_with_mapping<const D: usize, T: CoordNum>(
mpoint: &MultiPoint<D, T>,
mapping: impl Fn(&[T]) -> Vec<f64>,
) -> geojson::Geometry {
let mpoint_coords = mpoint
let coords = mpoint
.iter()
.map(|p| vertices[p[0] as usize].to_vec()) // Get the actual coord values
.map(&mapping) // Get the actual coord values
.collect();
geojson::Value::MultiPoint(mpoint_coords).into()
geojson::Value::MultiPoint(coords).into()
}

#[cfg(test)]
Expand Down Expand Up @@ -266,6 +286,33 @@ mod tests {
};
}

#[test]
fn test_multilinestring() {
let mut mls = MultiLineString::<3>::new();
mls.add_linestring([[11., 12., 13.], [21., 22., 23.], [31., 32., 33.]]);
mls.add_linestring([[111., 112., 113.], [121., 122., 123.], [131., 132., 133.]]);

let geom = multilinestring_to_geometry(&mls);
let geojson::Value::MultiLineString(mls) = geom.value else {
panic!("The result is not a GeoJSON MultiPolygon");

Check warning on line 297 in nusamai-geojson/src/conversion.rs

View check run for this annotation

Codecov / codecov/patch

nusamai-geojson/src/conversion.rs#L297 

Added line #L297 was not covered by tests
};
assert_eq!(
mls,
vec![
vec![
vec![11., 12., 13.],
vec![21., 22., 23.],
vec![31., 32., 33.],
],
vec![
vec![111., 112., 113.],
vec![121., 122., 123.],
vec![131., 132., 133.],
],
]
);
}

#[test]
fn test_indexed_multilinestring() {
let vertices = vec![
Expand Down Expand Up @@ -318,6 +365,27 @@ mod tests {
}
}

#[test]
fn test_multipoint() {
let mut mpoint = MultiPoint::<3>::new();
mpoint.push(&[11., 12., 13.]);
mpoint.push(&[21., 22., 23.]);
mpoint.push(&[31., 32., 33.]);

let geom = multipoint_to_geometry(&mpoint);
let geojson::Value::MultiPoint(mpoint) = geom.value else {
panic!("The result is not a GeoJSON MultiPolygon");

Check warning on line 377 in nusamai-geojson/src/conversion.rs

View check run for this annotation

Codecov / codecov/patch

nusamai-geojson/src/conversion.rs#L377 

Added line #L377 was not covered by tests
};
assert_eq!(
mpoint,
vec![
vec![11., 12., 13.],
vec![21., 22., 23.],
vec![31., 32., 33.],
]
);
}

#[test]
fn test_indexed_multipoint() {
let vertices = vec![[0., 0., 111.], [1., 2., 222.], [3., 4., 333.]];
Expand Down
164 changes: 55 additions & 109 deletions nusamai-gpkg/src/geometry.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@
//!
//! cf. https://www.geopackage.org/spec130/#gpb_format
use nusamai_geometry::{MultiPolygon, Polygon};
use nusamai_geometry::{CoordNum, MultiPolygon, Polygon};
use std::io::Write;

#[repr(u8)]
pub enum WkbByteOrder {
Expand Down Expand Up @@ -44,136 +45,80 @@ pub enum WkbGeometryType {
GeometryCollectionZM = 3007,
}

fn geometry_header(srs_id: i32) -> Vec<u8> {
let mut header: Vec<u8> = Vec::with_capacity(8);
header.extend_from_slice(&[0x47, 0x50]); // Magic number
header.push(0x00); // Version
header.push(0b00000001); // Flags
header.extend_from_slice(&i32::to_le_bytes(srs_id)); // SRS ID
header
fn write_geometry_header<W: Write>(writer: &mut W, srs_id: i32) -> std::io::Result<()> {
writer.write_all(&[0x47, 0x50])?; // Magic number
writer.write_all(&[
0x00, // Version
0b00000001, // Flags
])?;
writer.write_all(&i32::to_le_bytes(srs_id))?; // SRS ID
Ok(())
}

fn polygon_to_rings(vertices: &[[f64; 3]], poly: &Polygon<1, u32>) -> Vec<Vec<Vec<f64>>> {
let rings: Vec<_> = poly
.rings()
.map(|ls| {
let coords: Vec<_> = ls
.iter_closed()
.map(|idx| vertices[idx[0] as usize].to_vec()) // Get the actual coord values
.collect();
coords
})
.collect();

rings
fn write_polygon_body<W: Write, const D: usize, T: CoordNum>(
writer: &mut W,
poly: &Polygon<D, T>,
mapping: impl Fn(&[T]) -> [f64; 3],
) -> std::io::Result<()> {
// Byte order: Little endian (1)
writer.write_all(&[WkbByteOrder::LittleEndian as u8])?;

// Geometry type: wkbPolygonZ (1003)
writer.write_all(&(WkbGeometryType::PolygonZ as u32).to_le_bytes())?;

// numRings
writer.write_all(&(poly.rings().count() as u32).to_le_bytes())?;

for ring in poly.rings() {
// numPoints
writer.write_all(&(ring.iter_closed().count() as u32).to_le_bytes())?;

for idx in ring.iter_closed() {
let [x, y, z] = mapping(idx);
writer.write_all(&f64::to_le_bytes(x))?;
writer.write_all(&f64::to_le_bytes(y))?;
writer.write_all(&f64::to_le_bytes(z))?;
}
}
Ok(())
}

pub fn multipolygon_to_bytes(
pub fn write_indexed_multipolygon<W: Write>(
writer: &mut W,
vertices: &[[f64; 3]],
mpoly: &MultiPolygon<'_, 1, u32>,
srs_id: i32,
) -> Vec<u8> {
let mut bytes: Vec<u8> = geometry_header(srs_id);
) -> std::io::Result<()> {
write_geometry_header(writer, srs_id)?;
write_multipolygon_body(writer, mpoly, |idx| vertices[idx[0] as usize])?;
Ok(())
}

fn write_multipolygon_body<W: Write, const D: usize, T: CoordNum>(
writer: &mut W,
mpoly: &MultiPolygon<'_, D, T>,
mapping: impl Fn(&[T]) -> [f64; 3],
) -> std::io::Result<()> {
// Byte order: Little endian (1)
bytes.push(WkbByteOrder::LittleEndian as u8);
writer.write_all(&[WkbByteOrder::LittleEndian as u8])?;

// Geometry type: wkbMultiPolygonZ (1006)
bytes.extend_from_slice(&(WkbGeometryType::MultiPolygonZ as u32).to_le_bytes());
writer.write_all(&(WkbGeometryType::MultiPolygonZ as u32).to_le_bytes())?;

// numPolygons
bytes.extend_from_slice(&(mpoly.len() as u32).to_le_bytes());
writer.write_all(&(mpoly.len() as u32).to_le_bytes())?;

for poly in mpoly {
// Byte order: Little endian (1)
bytes.push(WkbByteOrder::LittleEndian as u8);

// Geometry type: wkbPolygonZ (1003)
bytes.extend_from_slice(&(WkbGeometryType::PolygonZ as u32).to_le_bytes());

let rings = polygon_to_rings(vertices, &poly);

// numRings
bytes.extend_from_slice(&(rings.len() as u32).to_le_bytes());

for ring in rings {
// numPoints
bytes.extend_from_slice(&(ring.len() as u32).to_le_bytes());

for coord in ring {
bytes.extend_from_slice(&f64::to_le_bytes(coord[0]));
bytes.extend_from_slice(&f64::to_le_bytes(coord[1]));
bytes.extend_from_slice(&f64::to_le_bytes(coord[2]));
}
}
write_polygon_body(writer, &poly, &mapping)?;
}

bytes
Ok(())
}

#[cfg(test)]
mod tests {
use super::*;

#[test]
fn test_polygon_to_rings() {
let vertices: Vec<[f64; 3]> = vec![
// exterior (vertex 0~3)
[0., 0., 111.],
[5., 0., 111.],
[5., 5., 111.],
[0., 5., 111.],
// interior 1 (vertex 4~7)
[1., 1., 111.],
[2., 1., 111.],
[2., 2., 111.],
[1., 2., 111.],
// interior 2 (vertex 8~11)
[3., 3., 111.],
[4., 3., 111.],
[4., 4., 111.],
[3., 4., 111.],
];

let mut poly = Polygon::<'_, 1, u32>::new();
poly.add_ring([[0], [1], [2], [3]]);
poly.add_ring([[4], [5], [6], [7]]);
poly.add_ring([[8], [9], [10], [11]]);

let rings = polygon_to_rings(&vertices, &poly);

assert_eq!(rings.len(), 3);

for (i, ri) in rings.iter().enumerate() {
match i {
0 => {
assert_eq!(ri.len(), 5);
assert_eq!(ri[0], vec![0., 0., 111.]);
assert_eq!(ri[1], vec![5., 0., 111.]);
assert_eq!(ri[2], vec![5., 5., 111.]);
assert_eq!(ri[3], vec![0., 5., 111.]);
assert_eq!(ri[4], vec![0., 0., 111.]);
}
1 => {
assert_eq!(ri.len(), 5);
assert_eq!(ri[0], vec![1., 1., 111.]);
assert_eq!(ri[1], vec![2., 1., 111.]);
assert_eq!(ri[2], vec![2., 2., 111.]);
assert_eq!(ri[3], vec![1., 2., 111.]);
assert_eq!(ri[4], vec![1., 1., 111.]);
}
2 => {
assert_eq!(ri.len(), 5);
assert_eq!(ri[0], vec![3., 3., 111.]);
assert_eq!(ri[1], vec![4., 3., 111.]);
assert_eq!(ri[2], vec![4., 4., 111.]);
assert_eq!(ri[3], vec![3., 4., 111.]);
}
_ => panic!("Unexpected ring index"),
}
}
}

#[test]
fn test_multipolygon_to_bytes() {
let vertices: Vec<[f64; 3]> = vec![
Expand All @@ -194,7 +139,8 @@ mod tests {
mpoly.add_exterior([[0], [1], [2], [3], [0]]);
mpoly.add_interior([[4], [5], [6], [7], [4]]);

let bytes = multipolygon_to_bytes(&vertices, &mpoly, 1234);
let mut bytes = Vec::new();
write_indexed_multipolygon(&mut bytes, &vertices, &mpoly, 1234).unwrap();

assert_eq!(bytes.len(), 274);

Expand Down
Loading

0 comments on commit eb90961

Please sign in to comment.