Skip to content

Commit

Permalink
added OSM source and adapted the tests
Browse files Browse the repository at this point in the history
set correct semver version

updated bitmask and added the creating/conversion scripts

updated bitmask and added the creating/conversion scripts

.

fixed tests and added images for both masks

fix errors related to paths
  • Loading branch information
TheSylex committed Sep 10, 2024
1 parent e156961 commit da0b180
Show file tree
Hide file tree
Showing 16 changed files with 257 additions and 144 deletions.
File renamed without changes.
File renamed without changes.
Binary file added assets/osm.wkb.xz
Binary file not shown.
Binary file added assets/osm_mask.tbmap.xz
Binary file not shown.
38 changes: 18 additions & 20 deletions build.rs
Original file line number Diff line number Diff line change
Expand Up @@ -4,35 +4,31 @@ use std::fs;
use std::io::prelude::*;
use std::path::Path;

pub static GSHHS_F: &str = "gshhs_f_-180.000000E-90.000000N180.000000E90.000000N.wkb.xz";
pub static GSHHS_F_CS: &str = "05bdf3089407b9829a7a5be7ee43f1e4205f2bbc641e4778af77e4814be216da";
pub static GSHHG: &str = "gshhg.wkb.xz";
pub static GSHHG_SHA256: &str = "05bdf3089407b9829a7a5be7ee43f1e4205f2bbc641e4778af77e4814be216da";

pub static GSHHG_MASK: &str = "gshhg_mask.tbmap.xz";
pub static GSHHG_MASK_SHA256: &str =
"5ea0e772ffc6ca8ad10c5de02be50670cbaedcff20b3541df6b78d3e1fdf48a1";
pub static GSHHG_MASK_SHA256: &str = "5ea0e772ffc6ca8ad10c5de02be50670cbaedcff20b3541df6b78d3e1fdf48a1";

pub static OSM: &str = "osm.wkb.xz";
pub static OSM_SHA256: &str = "7cbbbb56dc8f6a339d837e57aac4c50c9f54e7ac1118803274725cf61226b727";

pub static OSM_MASK: &str = "osm_mask.tbmap.xz";
pub static OSM_MASK_SHA256: &str =
"e60dd30737ad8480619d727bb246a1107d30a66563b73628337dc3f92255b684";
pub static OSM_MASK_SHA256: &str = "e60dd30737ad8480619d727bb246a1107d30a66563b73628337dc3f92255b684";

fn main() {
println!("hello");

let out_dir = env::var_os("OUT_DIR").unwrap();
println!("outdir: {:?}", out_dir);

let assets_dir = Path::new(&out_dir).join("gshhs");
let assets_dir = Path::new(&out_dir).join("assets");

// write assets script
let assets = Path::new(&out_dir).join("gshhs.rs");
{
let mut fd = fs::File::create(assets).unwrap();
write!(
fd,
"
let mut fd = fs::File::create(Path::new(&out_dir).join("source_data.rs")).unwrap();
write!(
fd,
"
use rust_embed::RustEmbed;
#[derive(RustEmbed)]
#[folder = \"{}\"]
Expand All @@ -44,17 +40,19 @@ pub struct OsmData;
",
assets_dir.to_slash().unwrap(),
assets_dir.to_slash().unwrap()
)
.unwrap();
).unwrap();

if !assets_dir.exists() {
fs::create_dir(assets_dir).unwrap();
}


// copy or download files
if env::var("DOCS_RS").is_err() {
copy_or_download(GSHHS_F, GSHHS_F_CS);
copy_or_download(MASK, MASK_CS);
copy_or_download(GSHHG, GSHHG_SHA256);
copy_or_download(GSHHG_MASK, GSHHG_MASK_SHA256);
copy_or_download(OSM, OSM_SHA256);
copy_or_download(OSM_MASK, OSM_MASK_SHA256);
} else {
println!("not downloading anything when on docs.rs.");
}
Expand All @@ -64,16 +62,16 @@ fn copy_or_download(from: impl AsRef<Path>, csum: &str) {
let from = from.as_ref();

let out_dir = env::var_os("OUT_DIR").unwrap();
let full_from = Path::new("gshhs").join(&from);
let full_to = Path::new(&out_dir).join("gshhs").join(&from);
let full_from = Path::new("assets").join(&from);
let full_to = Path::new(&out_dir).join("assets").join(&from);

if !full_to.exists() {
if full_from.exists() {
println!("copying {:?}..", &from);
fs::copy(&full_from, &full_to).unwrap();
} else {
let url = format!(
"https://github.com/gauteh/roaring-landmask/raw/main/gshhs/{}",
"https://github.com/gauteh/roaring-landmask/raw/main/assets/{}",
from.to_str().unwrap()
);
println!("downloading {:?} from {:?}..", &from, &url);
Expand Down
36 changes: 22 additions & 14 deletions src/devel/make_demo_plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,27 +3,35 @@
import numpy as np
import cartopy.crs as ccrs
from roaring_landmask import RoaringLandmask
from roaring_landmask import LandmaskProvider


fig = plt.figure(dpi = 150, figsize = (20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

ax.coastlines()
ax.gridlines(draw_labels = True)

x = np.arange(-180, 180, .5)
y = np.arange(-90, 90, .5)
x = np.arange(-180, 180, .1)
y = np.arange(-90, 90, .1)

xx, yy = np.meshgrid(x,y)

l = RoaringLandmask.new()
land = l.contains_many(xx.ravel(), yy.ravel())
gshhg_landmask = RoaringLandmask.new_with_provider(LandmaskProvider.Gshhg)
land = gshhg_landmask.contains_many(xx.ravel(), yy.ravel())
land = land.reshape(xx.shape)

cmap = matplotlib.colors.ListedColormap(['b', 'g'])
plt.pcolormesh(xx, yy, land, cmap = cmap, alpha = .7, transform = ccrs.PlateCarree())

plt.title('The Earth')
fig = plt.figure(dpi = 600, figsize = (20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())


plt.savefig('the_earth.png')
# plt.show()
ax.pcolormesh(xx, yy, land, cmap = cmap, alpha = .7, transform = ccrs.PlateCarree())
ax.gridlines(draw_labels = True)
ax.set_title('The Earth GSHHG', pad=20)
fig.savefig('the_earth_gshhg.png')

osm_landmask = RoaringLandmask.new_with_provider(LandmaskProvider.Osm)
land = osm_landmask.contains_many(xx.ravel(), yy.ravel())
land = land.reshape(xx.shape)

ax.clear()
ax.pcolormesh(xx, yy, land, cmap = cmap, alpha = .7, transform = ccrs.PlateCarree())
ax.gridlines(draw_labels = True)
ax.set_title('The Earth OSM', pad=20)
fig.savefig('the_earth_osm.png')
19 changes: 0 additions & 19 deletions src/devel/make_mask_old.py

This file was deleted.

113 changes: 113 additions & 0 deletions src/devel/roaring_scripts.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
import geopandas as gpd
import shapely
from shapely.geometry import MultiPolygon
from shapely.ops import unary_union
import numpy as np
import time
import regionmask
import xarray as xr


def wkb_to_wkt(wkb_file_path, wkt_file_path):
# Read the WKB file
with open(wkb_file_path, "rb") as wkb_file:
wkb_data = wkb_file.read()

# Parse WKB and create a shapely geometry object
geometry = shapely.wkb.loads(wkb_data)

# Convert geometry to WKT
wkt_data = shapely.wkt.dumps(geometry)

# Write WKT to file
with open(wkt_file_path, "w") as wkt_file:
wkt_file.write(wkt_data)

print(f"Conversion complete. WKT file saved as {wkt_file_path}")


# wkb_to_wkt('output.wkb', 'output.wkt')


def shapefile_to_wkb(shp_path, wkb_path, tolerance=0.0005, chunk_size=1000):
print("Reading shapefile...")
start_time = time.time()
gdf = gpd.read_file(shp_path)
print(f"Shapefile read in {time.time() - start_time:.2f} seconds.")
print(f"Number of geometries: {len(gdf)}")

simplified_geometries = []
total_geometries = len(gdf)

for i in range(0, total_geometries, chunk_size):
chunk = gdf.iloc[i : i + chunk_size]
print(
f"Processing chunk {i//chunk_size + 1} of {(total_geometries-1)//chunk_size + 1}..."
)

# Simplify each geometry individually
simplified_chunk = chunk.geometry.simplify(tolerance, preserve_topology=True)
simplified_geometries.extend(list(simplified_chunk))

print(
f"Processed {min(i+chunk_size, total_geometries)} out of {total_geometries} geometries."
)

print("Combining all simplified geometries...")
# Use unary_union to combine all geometries efficiently
combined_geometry = unary_union(simplified_geometries)

# Ensure the result is a MultiPolygon
if combined_geometry.geom_type == "Polygon":
final_multi_polygon = MultiPolygon([combined_geometry])
elif combined_geometry.geom_type == "MultiPolygon":
final_multi_polygon = combined_geometry
else:
raise ValueError(f"Unexpected geometry type: {combined_geometry.geom_type}")

print("Getting WKB representation...")
wkb = final_multi_polygon.wkb

print("Writing WKB to file...")
with open(wkb_path, "wb") as wkb_file:
wkb_file.write(wkb)

print(f"Process completed. Output written to: {wkb_path}")


# shapefile_to_wkb('land-polygons-complete-4326/land_polygons.shp', 'output.wkb')


def wkb_to_mask(wkb_file_path, chunk_size=1000):
# Read the WKB file as binary
with open(wkb_file_path, "rb") as file:
wkb_data = file.read()

# Convert WKB to Shapely geometry
geometry = shapely.wkb.loads(wkb_data)

# Create a GeoDataFrame
gdf = gpd.GeoDataFrame({"geometry": [geometry]})

lon = np.linspace(-180, 180, 86400, endpoint=False)
lat = np.linspace(-90, 90, 43200, endpoint=False)

total_rows = len(lat)

with open("mask.bin", "wb") as f:
for start_row in range(0, total_rows, chunk_size):
end_row = min(start_row + chunk_size, total_rows)
lat_chunk = lat[start_row:end_row]

mask_chunk = regionmask.mask_geopandas(gdf, lon, lat_chunk)
mask_chunk = xr.where(mask_chunk == 0.0, 1.0, mask_chunk)
mask_chunk = mask_chunk.fillna(0.0)
mask_chunk = mask_chunk.astype(np.uint8)

mask_chunk.values.tofile(f)
print(
f"{end_row / total_rows * 100.0}% - ({start_row} -> {end_row}) / {total_rows}"
)


# wkb_to_mask("output.wkb")
Loading

0 comments on commit da0b180

Please sign in to comment.