Give a geographic location, download and read LIDAR from GeoNB

.gitignore

@@ -0,0 +1,3 @@
+/target
+*~
+\#*#

Cargo.toml

@@ -0,0 +1,19 @@
+[package]
+name = "pteropus"
+version = "0.1.0"
+edition = "2021"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
+serde_json = "~1.0"
+tokio = {version="~1.17", features=["rt-multi-thread", "macros"]}
+tokio-util = {version="~0.7", features=["io"]}
+reqwest = {version="~0.11", features=["json", "stream"]}
+proj = "~0.25"
+geo-types = "~0.7"
+clap = {version="~3.1", features=["derive"]}
+anyhow = "~1.0"
+futures = "~0.3"
+bytes = "~1.1"
+las = {version="~0.7", features=["laz"]}

src/geonb.rs

@@ -0,0 +1,111 @@
+use {futures::stream::TryStreamExt, geo_types::Point, tokio_util::io::StreamReader};
+
+#[derive(Debug)]
+pub struct Error {
+    pub message: String,
+}
+
+impl Error {
+    pub fn new(message: &str) -> Error {
+        let message = message.to_string();
+        Error { message }
+    }
+}
+
+impl std::fmt::Display for Error {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> Result<(), std::fmt::Error> {
+        write!(f, "GeoNB Error: \"{}\"", self.message)
+    }
+}
+
+impl std::error::Error for Error {}
+fn convert_err(err: reqwest::Error) -> std::io::Error {
+    todo!()
+}
+pub async fn get_lidar_tile_around_point(
+    location: Point<f64>,
+) -> Result<impl tokio::io::AsyncRead, anyhow::Error> {
+    let client = reqwest::Client::new();
+    let query_response_json = &client
+        .get("https://geonb.snb.ca/arcgis/rest/services/GeoNB_SNB_LidarIndex/MapServer/1/query")
+        .query(&[
+            ("f", "json"),
+            ("geometryType", "esriGeometryPoint"),
+            ("geometry", &format!("{},{}", location.x(), location.y())),
+            ("returnIdsOnly", "true"),
+        ])
+        .send()
+        .await?
+        .json::<serde_json::Value>()
+        .await?;
+
+    println!(
+        "{}",
+        serde_json::to_string_pretty(&query_response_json).expect("JSON")
+    );
+
+    let object_id = query_response_json
+        .get("objectIds")
+        .and_then(|object_ids| object_ids.get(0))
+        .and_then(&serde_json::Value::as_i64)
+        .ok_or_else(|| Error::new("Could not find \"objectId\" in response."))?;
+
+    let object_response_json = client
+        .get(format!(
+            "https://geonb.snb.ca/arcgis/rest/services/GeoNB_SNB_LidarIndex/MapServer/1/{}",
+            object_id
+        ))
+        .query(&[("f", "json")])
+        .send()
+        .await?
+        .json::<serde_json::Value>()
+        .await?;
+
+    println!(
+        "{}",
+        serde_json::to_string_pretty(&object_response_json).expect("JSON")
+    );
+
+    let laz_file_url = object_response_json
+        .get("feature")
+        .ok_or_else(|| Error::new("Could not find \"feature\" in response."))?
+        .get("attributes")
+        .ok_or_else(|| Error::new("Could not find \"attributes\" in response."))?
+        .get("FILE_URL")
+        .ok_or_else(|| Error::new("Could not find \"FILE_URL\" in response."))?
+        .as_str()
+        .ok_or_else(|| Error::new("Expected \"FILE_URL\" to be a string but it was not."))?;
+
+    println!("LAZ URL: {}", laz_file_url);
+    Ok(StreamReader::new(
+        client
+            .get(laz_file_url)
+            .send()
+            .await?
+            .bytes_stream()
+            .map_err(convert_err),
+    ))
+}
+
+pub async fn test() -> Result<(), reqwest::Error> {
+    let client = reqwest::Client::new();
+    print!(
+        "{}",
+        serde_json::to_string_pretty(&client
+            .get("https://geonb.snb.ca/arcgis/rest/services/GeoNB_SNB_LidarIndex/MapServer/1/query")
+            .query(&[("f", "json"),("geometryType","esriGeometryPoint"),("geometry","2470000,7443000"),("returnIdsOnly","true")])
+            .send()
+            .await?
+            .json::<serde_json::Value>().await?).expect("JSON")
+    );
+    print!(
+        "{}",
+        serde_json::to_string_pretty(&client
+            .get("https://geonb.snb.ca/arcgis/rest/services/GeoNB_SNB_LidarIndex/MapServer/1/14601")
+            .query(&[("f", "json")])
+            .send()
+            .await?
+            .json::<serde_json::Value>().await?).expect("JSON")
+    );
+    Ok(())
+}

src/main.rs

@@ -0,0 +1,48 @@
+use {clap::Parser, geo_types::Point, las::Read as LasRead, proj::Proj, tokio::io::AsyncReadExt};
+
+mod geonb;
+
+#[derive(Parser, Debug)]
+#[clap(author, version, about, long_about=None)]
+struct Args {
+    // Latitude to fetch LIDAR tile at
+    #[clap(long)]
+    latitude: f64,
+
+    // Longitude to fetch LIDAR tile at
+    #[clap(long)]
+    longitude: f64,
+}
+
+#[tokio::main]
+async fn main() -> Result<(), anyhow::Error> {
+    let args = Args::parse();
+
+    let location = Proj::new_known_crs("+proj=longlat +datum=WGS84", "EPSG:2953", None)
+        .unwrap()
+        .convert(Point::new(args.longitude, args.latitude))
+        .unwrap();
+    println!("{:?}", location);
+    let mut las_reader =
+        tokio::io::BufReader::new(geonb::get_lidar_tile_around_point(location).await?);
+    let mut las_bytes = Vec::new();
+    let mut buffer = [0_u8; 4096];
+    let mut byte_count = 0;
+    loop {
+        let num_bytes = las_reader.read(&mut buffer).await?;
+        if num_bytes == 0 {
+            break;
+        }
+        byte_count += num_bytes;
+        print!("{} bytes read\r", byte_count);
+        las_bytes.extend_from_slice(&buffer[0..num_bytes]);
+    }
+    println!();
+    let mut las_reader = las::Reader::new(std::io::Cursor::new(las_bytes))?;
+    for wrapped_point in las_reader.points().take(10) {
+        let point = wrapped_point.unwrap();
+        println!("Point coordinates: ({}, {}, {})", point.x, point.y, point.z);
+    }
+
+    Ok(())
+}