Skip to content

Commit

Permalink
add proj params, fix python dataset encoding
Browse files Browse the repository at this point in the history
  • Loading branch information
mpiannucci committed Jul 17, 2023
1 parent bb80c04 commit cd31a3d
Show file tree
Hide file tree
Showing 10 changed files with 454 additions and 165 deletions.
2 changes: 1 addition & 1 deletion python/examples/hrrr/hrrr_kerchunk.json

Large diffs are not rendered by default.

403 changes: 290 additions & 113 deletions python/examples/kerchunk_hrrr_subhourly.ipynb

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions python/gribberish/gribberish_backend.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
from xarray.backends.common import BackendEntrypoint, BackendArray
from xarray.core import indexing

from .gribberishpy import parse_grid_dataset, parse_grib_array
from .gribberishpy import parse_grib_dataset, parse_grib_array


DATA_VAR_LOCK = xr.backends.locks.SerializableLock()
Expand All @@ -31,7 +31,7 @@ def open_dataset(
with open(filename_or_obj, 'rb') as f:
raw_data = f.read()

dataset = parse_grid_dataset(
dataset = parse_grib_dataset(
raw_data,
drop_variables=drop_variables,
only_variables=only_variables,
Expand Down
4 changes: 2 additions & 2 deletions python/gribberish/kerchunk/mapper.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

from kerchunk.utils import class_factory, _encode_for_JSON
from .codec import GribberishCodec
from ..gribberishpy import parse_grid_dataset
from ..gribberishpy import parse_grib_dataset


def _split_file(f):
Expand Down Expand Up @@ -113,7 +113,7 @@ def scan_gribberish(
with fsspec.open(url, "rb", **storage_options) as f:
for offset, size, data in _split_file(f):
try:
dataset = parse_grid_dataset(data, perserve_dims=perserve_dims)
dataset = parse_grib_dataset(data, perserve_dims=perserve_dims, encode_coords=True)
except Exception as e:
# Skip messages that gribberish cannot handle yet
continue
Expand Down
53 changes: 42 additions & 11 deletions python/src/dataset.rs
Original file line number Diff line number Diff line change
Expand Up @@ -48,13 +48,14 @@ pub fn build_grib_array<'py>(
}

#[pyfunction]
pub fn parse_grid_dataset<'py>(
pub fn parse_grib_dataset<'py>(
py: Python<'py>,
data: &[u8],
drop_variables: Option<&PyList>,
only_variables: Option<&PyList>,
perserve_dims: Option<&PyList>,
filter_by_attrs: Option<&PyDict>,
encode_coords: Option<bool>,
) -> PyResult<&'py PyDict> {
let drop_variables = if let Some(drop_variables) = drop_variables {
drop_variables
Expand Down Expand Up @@ -91,6 +92,12 @@ pub fn parse_grid_dataset<'py>(
PyDict::new(py)
};

let encode_coords = if let Some(encode_coords) = encode_coords {
encode_coords
} else {
false
};

let mapping = scan_message_metadata(data)
.into_iter()
.filter_map(|(k, v)| {
Expand Down Expand Up @@ -469,17 +476,28 @@ pub fn parse_grid_dataset<'py>(
coords.set_item("x", x).unwrap();

latitude.set_item("dims", vec!["y", "x"]).unwrap();
let lats_array = PyDict::new(py);
lats_array.set_item("shape", [grid_shape.0, grid_shape.1]).unwrap();
lats_array.set_item("offsets", [first.1]).unwrap();
latitude.set_item("values", lats_array).unwrap();
latitude.set_item("shape", [grid_shape.0, grid_shape.1]).unwrap();

longitude.set_item("dims", vec!["y", "x"]).unwrap();
let lngs_array = PyDict::new(py);
lngs_array.set_item("shape", [grid_shape.0, grid_shape.1]).unwrap();
lngs_array.set_item("offsets", [first.1]).unwrap();
longitude.set_item("values", lngs_array).unwrap();

if encode_coords {
let lats_array = PyDict::new(py);
lats_array.set_item("shape", [grid_shape.0, grid_shape.1]).unwrap();
lats_array.set_item("offsets", [first.1]).unwrap();
latitude.set_item("values", lats_array).unwrap();

let lngs_array = PyDict::new(py);
lngs_array.set_item("shape", [grid_shape.0, grid_shape.1]).unwrap();
lngs_array.set_item("offsets", [first.1]).unwrap();
longitude.set_item("values", lngs_array).unwrap();
} else {
let (lat, lng) = first.2.latlng();
let lats = PyArray::from_slice(py, &lat);
let lats = lats.reshape([grid_shape.0, grid_shape.1]).unwrap();
latitude.set_item("values", lats).unwrap();

let lngs = PyArray::from_slice(py, &lng);
let lngs = lngs.reshape([grid_shape.0, grid_shape.1]).unwrap();
longitude.set_item("values", lngs).unwrap();
}

var_dims.iter_mut().for_each(|(_, v)| {
v.push("y".to_string());
Expand Down Expand Up @@ -530,6 +548,19 @@ pub fn parse_grid_dataset<'py>(
.unwrap_or("".to_string()),
)
.unwrap();

let proj_params = PyDict::new(py);
proj_params
.set_item("proj", first.2.projector.proj_name())
.unwrap();
first.2.projector
.proj_params()
.iter()
.for_each(|(k, v)| {
proj_params.set_item(k, v).unwrap();
});
var_metadata.set_item("proj_params", proj_params).unwrap();

var_metadata.set_item("crs", first.2.proj.clone()).unwrap();

let mut filtered = false;
Expand Down
4 changes: 2 additions & 2 deletions python/src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ use pyo3::prelude::*;
use pyo3::wrap_pyfunction;

use crate::dataset::build_grib_array;
use crate::dataset::parse_grid_dataset;
use crate::dataset::parse_grib_dataset;
use crate::message::parse_grib_array;
use crate::message::parse_grib_mapping;
use crate::message::parse_grib_message;
Expand All @@ -20,7 +20,7 @@ fn gribberishpy(_py: Python, m: &PyModule) -> PyResult<()> {
m.add_function(wrap_pyfunction!(parse_grib_message, m)?)?;
m.add_function(wrap_pyfunction!(parse_grib_messages, m)?)?;
m.add_function(wrap_pyfunction!(parse_grib_mapping, m)?)?;
m.add_function(wrap_pyfunction!(parse_grid_dataset, m)?)?;
m.add_function(wrap_pyfunction!(parse_grib_dataset, m)?)?;
m.add_function(wrap_pyfunction!(parse_grib_array, m)?)?;
m.add_function(wrap_pyfunction!(build_grib_array, m)?)?;
Ok(())
Expand Down
4 changes: 4 additions & 0 deletions src/templates/grid_definition/grid_definition_template.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,11 @@
use std::collections::HashMap;

use crate::utils::iter::projection::LatLngProjection;

pub trait GridDefinitionTemplate {
fn proj_string(&self) -> String;
fn proj_name(&self) -> String;
fn proj_params(&self) -> HashMap<String, f64>;
fn crs(&self) -> String;
fn grid_point_count(&self) -> usize;
fn is_regular_grid(&self) -> bool;
Expand Down
26 changes: 24 additions & 2 deletions src/templates/grid_definition/lambert_conformal_template.rs
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ use crate::{
templates::template::{Template, TemplateType},
utils::{
bit_array_from_bytes,
iter::projection::{LatLngProjection, RegularCoordinateIterator},
iter::projection::{LatLngProjection, RegularCoordinateIterator, LambertConformalConicProjection},
read_u32_from_bytes,
},
};
Expand Down Expand Up @@ -222,6 +222,22 @@ impl LambertConformalTemplate {
}

impl GridDefinitionTemplate for LambertConformalTemplate {
fn proj_name(&self) -> String {
"lcc".to_string()
}

fn proj_params(&self) -> std::collections::HashMap<String, f64> {
let mut params = std::collections::HashMap::new();
params.insert(
"lon_0".to_string(),
self.longitude_of_paralell_meridian_to_latitude_increase(),
);
params.insert("lat_0".to_string(), self.latitude_of_dx_dy());
params.insert("lat_1".to_string(), self.latin_1());
params.insert("lat_2".to_string(), self.latin_2());
params
}

fn proj_string(&self) -> String {
format!(
"+proj=lcc lon_0={} lat_0={} lat_1={} lat_2={}",
Expand Down Expand Up @@ -283,6 +299,12 @@ impl GridDefinitionTemplate for LambertConformalTemplate {
self.number_of_points_on_x_axis() as usize,
);

LatLngProjection::LambertConformal(y_iter, x_iter, projection)
LatLngProjection::LambertConformal(LambertConformalConicProjection{
x: x_iter,
y: y_iter,
projection,
projection_name: self.proj_name(),
projection_params: self.proj_params(),
})
}
}
20 changes: 18 additions & 2 deletions src/templates/grid_definition/latlng_template.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
use super::grid_definition_template::GridDefinitionTemplate;
use super::tables::{EarthShape, ScanningModeFlags, ScanningMode};
use crate::templates::template::{Template, TemplateType};
use crate::utils::iter::projection::{RegularCoordinateIterator, LatLngProjection};
use crate::utils::iter::projection::{RegularCoordinateIterator, LatLngProjection, PlateCareeProjection};
use crate::utils::{bit_array_from_bytes, read_u32_from_bytes};

use std::iter::Iterator;
Expand Down Expand Up @@ -153,6 +153,17 @@ impl LatLngTemplate {
}

impl GridDefinitionTemplate for LatLngTemplate {
fn proj_name(&self) -> String {
"latlon".to_string()
}

fn proj_params(&self) -> std::collections::HashMap<String, f64> {
let mut params = std::collections::HashMap::new();
params.insert("a".to_string(), 6367470.0);
params.insert("b".to_string(), 6367470.0);
params
}

fn proj_string(&self) -> String {
format!("+proj=latlon +a=6367470 +b=6367470")
}
Expand Down Expand Up @@ -190,6 +201,11 @@ impl GridDefinitionTemplate for LatLngTemplate {
self.x_count()
);

LatLngProjection::PlateCaree(lat_iter, lon_iter)
LatLngProjection::PlateCaree(PlateCareeProjection {
latitudes: lat_iter,
longitudes: lon_iter,
projection_name: self.proj_name(),
projection_params: self.proj_params(),
})
}
}
Loading

0 comments on commit cd31a3d

Please sign in to comment.