From 06b7dcc98e677ec4b22f26ebcd009654be391cbb Mon Sep 17 00:00:00 2001 From: max Date: Fri, 20 Jun 2025 15:41:32 +1000 Subject: [PATCH 1/3] fixing lon lat order --- src/accessvis/earth.py | 107 +++++++++++++++++++++++-------------- tests/test_earth_no_gpu.py | 54 +++++++++---------- 2 files changed, 93 insertions(+), 68 deletions(-) diff --git a/src/accessvis/earth.py b/src/accessvis/earth.py index f1c1f55..ab85de6 100644 --- a/src/accessvis/earth.py +++ b/src/accessvis/earth.py @@ -309,7 +309,7 @@ def split_tex(data, res, flipud=False, fliplr=False): return tiles -def draw_latlon_grid(base_img, out_fn, lat=30, lon=30, linewidth=5, colour=0): +def draw_lonlat_grid(base_img, out_fn, lon=30, lat=30, linewidth=5, colour=0): """ Create lat/lon grid image over provided base image """ @@ -343,7 +343,7 @@ def draw_latlon_grid(base_img, out_fn, lat=30, lon=30, linewidth=5, colour=0): image.save(out_fn) -def latlon_to_uv(lat, lon): +def lonlat_to_uv(lon, lat): """ Convert a decimal longitude, latitude coordinate to a tex coord in an equirectangular image @@ -365,12 +365,12 @@ def uv_to_pixel(u, v, width, height): return int(u * width), int(v * height) -def latlon_to_pixel(lat, lon, width, height): +def lonlat_to_pixel(lon, lat, width, height): """ Convert a decimal latitude/longitude coordinate to a raster image pixel coordinate for given width/height """ - u, v = latlon_to_uv(lat, lon) + u, v = lonlat_to_uv(lon=lon, lat=lat) return uv_to_pixel(u, v, width, height) @@ -472,13 +472,13 @@ def crop_img_uv(img, cropbox): print("Unknown type: ", type(img)) -def crop_img_lat_lon(img, cropbox): +def crop_img_lon_lat(img, top_left, bottom_right): """ Crop an equirectangular image (PIL or numpy array) based on corner coords Provide coords as lat/lon coords in decimal degrees """ - a = latlon_to_uv(*cropbox[0]) - b = latlon_to_uv(*cropbox[1]) + a = lonlat_to_uv(*top_left) + b = lonlat_to_uv(*bottom_right) return crop_img_uv(img, (a, b)) @@ -754,7 +754,9 @@ def load_topography_cubemap( ) -def load_topography(resolution=None, subsample=1, cropbox=None, bathymetry=True): +def load_topography( + resolution=None, subsample=1, top_left=None, bottom_right=None, bathymetry=True +): """ Load topography from pre-saved equirectangular data, can be cropped for regional plots @@ -775,8 +777,10 @@ def load_topography(resolution=None, subsample=1, cropbox=None, bathymetry=True) if subsample > 1: heights = heights[::subsample, ::subsample] - if cropbox: - heights = crop_img_lat_lon(heights, cropbox) + if top_left and bottom_right: + heights = crop_img_lon_lat( + heights, top_left=top_left, bottom_right=bottom_right + ) # Bathymetry? if not bathymetry or not isinstance(bathymetry, bool): @@ -784,7 +788,11 @@ def load_topography(resolution=None, subsample=1, cropbox=None, bathymetry=True) # res_y = resolution//4096 * 10800 # res_y = max(resolution,2048) // 2048 * 10800 mask = load_mask( - res_y=resolution, subsample=subsample, cropbox=cropbox, masktype="oceanmask" + res_y=resolution, + subsample=subsample, + top_left=top_left, + bottom_right=bottom_right, + masktype="oceanmask", ) # print(type(mask), mask.dtype, mask.min(), mask.max()) if bathymetry == "mask": @@ -804,7 +812,8 @@ def load_topography(resolution=None, subsample=1, cropbox=None, bathymetry=True) def plot_region( lv=None, - cropbox=None, + top_left=None, + bottom_right=None, vertical_exaggeration=10, texture="bluemarble", lighting=True, @@ -881,7 +890,9 @@ def plot_region( uniforms[k] = kwargs[k] # Load topo and crop via lat/lon boundaries of data - topo = load_topography(cropbox=cropbox, bathymetry=bathymetry) + topo = load_topography( + top_left=top_left, bottom_right=bottom_right, bathymetry=bathymetry + ) height = np.array(topo) # Scale and apply vertical exaggeration @@ -889,8 +900,8 @@ def plot_region( D = [height.shape[1], height.shape[0]] sverts = np.zeros(shape=(height.shape[0], height.shape[1], 3)) - lat0, lon0 = cropbox[0] if cropbox else [-90, 0] - lat1, lon1 = cropbox[1] if cropbox else [90, 360] + lat0, lon0 = top_left if top_left else [-90, 0] + lat1, lon1 = bottom_right if bottom_right else [90, 360] # TODO: Support crossing zero in longitude # Will probably not support crossing poles xy = lv.grid2d(corners=((lon0, lat1), (lon1, lat0)), dims=D) @@ -922,8 +933,10 @@ def plot_region( ) # , fliptexture=False) img = Image.open(colour_tex) - if cropbox: - cropped_img = crop_img_lat_lon(img, cropbox) + if top_left and bottom_right: + cropped_img = crop_img_lon_lat( + img, top_left=top_left, bottom_right=bottom_right + ) arr = np.array(cropped_img) else: arr = np.array(img) @@ -1455,31 +1468,31 @@ def normalise(vec): return vn -def latlon_normal_vector(lat, lon): +def lonlat_normal_vector(lon, lat): """ If standing at (lat,lon), this unit vector points directly above. """ - return normalise(latlon_to_3D(lat=lat, lon=lon)) + return normalise(lonlat_to_3D(lat=lat, lon=lon)) -def latlon_vector_to_east(lat, lon): +def lonlat_vector_to_east(lon, lat): """ If standing at (lat,lon), this unit vector points east. """ - start = latlon_to_3D(lat=lat, lon=lon) - north = latlon_to_3D(lat=90, lon=0) + start = lonlat_to_3D(lat=lat, lon=lon) + north = lonlat_to_3D(lat=90, lon=0) cross = np.cross(north, start) if np.any(cross): # if start != north return normalise(cross) return np.array([1, 0, 0]) -def latlon_vector_to_north(lat, lon): +def lonlat_vector_to_north(lon, lat): """ If standing at (lat,lon), this unit vector points north. """ - start = latlon_to_3D(lat=lat, lon=lon) - east = latlon_vector_to_east(lat=lat, lon=lon) + start = lonlat_to_3D(lat=lat, lon=lon) + east = lonlat_vector_to_east(lat=lat, lon=lon) return normalise(np.cross(start, east)) @@ -1802,7 +1815,9 @@ def array_to_rgba( bm_tiles = ["A1", "B1", "C1", "D1", "A2", "B2", "C2", "D2"] -def load_mask(res_y=None, masktype="watermask", subsample=1, cropbox=None): +def load_mask( + res_y=None, masktype="watermask", subsample=1, top_left=None, bottom_right=None +): # TODO: Document args """ Loads watermask/oceanmask @@ -1833,8 +1848,8 @@ def load_mask(res_y=None, masktype="watermask", subsample=1, cropbox=None): if subsample > 1: mask = mask[::subsample, ::subsample] - if cropbox: - return crop_img_lat_lon(mask, cropbox) + if top_left and bottom_right: + return crop_img_lon_lat(mask, top_left=top_left, bottom_right=bottom_right) return mask @@ -2420,17 +2435,20 @@ def plot_vectors_xr( alt = max_alt # Basis vector directions, on the 3d model. - normal = latlon_normal_vector(lat=lat, lon=lon) - east = latlon_vector_to_east(lat=lat, lon=lon) - north = latlon_vector_to_north(lat=lat, lon=lon) + normal = lonlat_normal_vector(lat=lat, lon=lon) + east = lonlat_vector_to_east(lat=lat, lon=lon) + north = lonlat_vector_to_north(lat=lat, lon=lon) if rescale_alt: # The bottom-most vector is ground level, top-most is at a height of alt_scale_factor above ground relative_height = (alt - min_alt) / (max_alt - min_alt) - vert = latlon_to_3D(lat, lon) + normal * relative_height * alt_scale_factor + vert = ( + lonlat_to_3D(lat=lat, lon=lon) + + normal * relative_height * alt_scale_factor + ) vertices.append(vert) else: - vert = latlon_to_3D(lat=lat, lon=lon, alt=alt * 1e-6) + vert = lonlat_to_3D(lat=lat, lon=lon, alt=alt * 1e-6) # 1e-6 is to convert from meters to Million meters (3d Scale) vertices.append(vert) @@ -2524,7 +2542,7 @@ def plot_shapefile(lv, fn, features=None, alt=1e-6, label="shape_", **kwargs): # converting to 3d Coordinates. lons, lats = zip(*shape.points) - verts = latlon_to_3D(lat=lats, lon=lons, alt=alt).T + verts = lonlat_to_3D(lat=lats, lon=lons, alt=alt).T # Iterate through each part, e.g. each island. for i, (p1, p2) in enumerate(zip(shape.parts[:-1], shape.parts[1:])): @@ -2597,8 +2615,8 @@ def plot_cross_section( # Calculating the position of all vertices. lats = np.linspace(start[1], end[1], resolution) lons = np.linspace(start[0], end[0], resolution) - lower = latlon_to_3D(lat=lats, lon=lons, alt=min_alt) - upper = latlon_to_3D(lat=lats, lon=lons, alt=max_alt) + lower = lonlat_to_3D(lat=lats, lon=lons, alt=min_alt) + upper = lonlat_to_3D(lat=lats, lon=lons, alt=max_alt) vertices = np.dstack((lower, upper)).T surf.vertices(vertices) @@ -2608,7 +2626,7 @@ def plot_cross_section( return surf -def plot_region_data(lv, data, start, end, alt=0.1, label="data-surface", **kwargs): +def plot_region_data(lv, data, start, end, label="data-surface", **kwargs): """ Plots data on a 2D region of the earth. Use this to plot data after plot_region() is called. @@ -2621,11 +2639,12 @@ def plot_region_data(lv, data, start, end, alt=0.1, label="data-surface", **kwar data: np.ndarray An array of shape [width, height, RGB(A)]. start: tuple[number, number] - (lon, lat) - data[0,0] corrosponds to the start position. + (lon, lat) or (lon, lat, alt) + data[0,0] corresponds to the start position. + alt is optional and defaults to 1e-5. end: tuple[number, number] (lon, lat) - data[-1,-1] corrosponds to the end position. + data[-1,-1] corresponds to the end position. alt: number height above sea level to display the data. label: str @@ -2640,6 +2659,14 @@ def plot_region_data(lv, data, start, end, alt=0.1, label="data-surface", **kwar lavavu.Object The lavavu surface being created. """ + + if len(start) == 3: + alt = start[2] + elif len(end) == 3: + alt = end[2] + else: + alt = 1e-5 + surf = lv.triangles( label, colour="rgba(255,255,255,0)", **kwargs # allows transparent data ) diff --git a/tests/test_earth_no_gpu.py b/tests/test_earth_no_gpu.py index c8a82cd..c5ca566 100644 --- a/tests/test_earth_no_gpu.py +++ b/tests/test_earth_no_gpu.py @@ -67,7 +67,7 @@ def test_normalise_array(): assert np.allclose(arr3, exp3) -def test_latlon_vector_to_north(): +def test_lonlat_vector_to_north(): """ At the equator, all point directly upwards. Should have length 1. @@ -75,54 +75,54 @@ def test_latlon_vector_to_north(): import accessvis import numpy as np - vec1 = accessvis.latlon_vector_to_north(lat=0.0, lon=0.0) - vec2 = accessvis.latlon_vector_to_north(lat=0, lon=90) + vec1 = accessvis.lonlat_vector_to_north(lat=0.0, lon=0.0) + vec2 = accessvis.lonlat_vector_to_north(lat=0, lon=90) vec3 = [0, 1, 0] assert np.allclose(vec1, vec2) assert np.allclose(vec1, vec3) - vec4 = accessvis.latlon_vector_to_north(lat=47.3, lon=88.1) + vec4 = accessvis.lonlat_vector_to_north(lat=47.3, lon=88.1) assert np.allclose(np.linalg.norm(vec4), 1), "doesn't have length 1." -def test_latlon_vector_to_east(): +def test_lonlat_vector_to_east(): import accessvis import numpy as np - vec1 = accessvis.latlon_vector_to_east(lat=0.0, lon=0.0) + vec1 = accessvis.lonlat_vector_to_east(lat=0.0, lon=0.0) vec2 = [1, 0, 0] assert np.allclose(vec1, vec2) - vec3 = accessvis.latlon_vector_to_east(lat=0, lon=-90) + vec3 = accessvis.lonlat_vector_to_east(lat=0, lon=-90) vec4 = [0, 0, 1] assert np.allclose(vec3, vec4) - norm_vec = accessvis.latlon_vector_to_east(lat=87.3, lon=-25.1) + norm_vec = accessvis.lonlat_vector_to_east(lat=87.3, lon=-25.1) assert np.allclose(np.linalg.norm(norm_vec), 1), "doesn't have length 1." -def test_latlon_normal_vector(): +def test_lonlat_normal_vector(): import accessvis import numpy as np - vec1 = accessvis.latlon_normal_vector(lat=0.0, lon=0.0) + vec1 = accessvis.lonlat_normal_vector(lat=0.0, lon=0.0) vec2 = [0, 0, 1] assert np.allclose(vec1, vec2) - vec3 = accessvis.latlon_normal_vector(lat=0, lon=90) + vec3 = accessvis.lonlat_normal_vector(lat=0, lon=90) vec4 = [1, 0, 0] assert np.allclose(vec3, vec4) - vec5 = accessvis.latlon_normal_vector(lat=90, lon=47.2) + vec5 = accessvis.lonlat_normal_vector(lat=90, lon=47.2) vec6 = [0, 1, 0] assert np.allclose(vec5, vec6) - vec7 = accessvis.latlon_normal_vector(lat=-90, lon=47.2) + vec7 = accessvis.lonlat_normal_vector(lat=-90, lon=47.2) vec8 = [0, -1, 0] assert np.allclose(vec7, vec8) - norm_vec = accessvis.latlon_normal_vector(lat=-65.1, lon=-48.1) + norm_vec = accessvis.lonlat_normal_vector(lat=-65.1, lon=-48.1) assert np.allclose(np.linalg.norm(norm_vec), 1), "doesn't have length 1." @@ -252,25 +252,25 @@ def test_crop_img_lat_lon(): arr[:, :10] = 1 arr[:, 10:] = 2 - cropbox1 = (90, -180), (-90, 0) # left side - out1 = accessvis.crop_img_lat_lon(img=arr, cropbox=cropbox1) + tl1, br1 = (-180, 90), (0, -90) # left side + out1 = accessvis.crop_img_lon_lat(img=arr, top_left=tl1, bottom_right=br1) assert out1.shape == (10, 10), f"Left shape {out1.shape}" assert np.all(out1 == 1), "Left sum" - cropbox2 = (45, -90), (-45, 90) # Middle - out2 = accessvis.crop_img_lat_lon(img=arr, cropbox=cropbox2) + tl2, br2 = (-90, 45), (90, -45) # Middle + out2 = accessvis.crop_img_lon_lat(img=arr, top_left=tl2, bottom_right=br2) assert out2.shape == (5, 10), f"Middle shape {out2.shape}" assert np.all(out2[:, :5] == 1), "middle 1" assert np.all(out2[:, 5:] == 2), "middle 2" - cropbox3 = (90, -270), (-90, -90) # overflow left - out3 = accessvis.crop_img_lat_lon(img=arr, cropbox=cropbox3) + tl3, br3 = (-270, 90), (-90, -90) # overflow left + out3 = accessvis.crop_img_lon_lat(img=arr, top_left=tl3, bottom_right=br3) assert out3.shape == (10, 10), "overflow left shape" assert np.all(out3[:, :5] == 2), "overflow left 1" assert np.all(out3[:, 5:] == 1), "overflow left 2" - cropbox4 = (90, 90), (-90, 270) # overflow right - out4 = accessvis.crop_img_lat_lon(img=arr, cropbox=cropbox4) + tl4, br4 = (90, 90), (270, -90) # overflow right + out4 = accessvis.crop_img_lon_lat(img=arr, top_left=tl4, bottom_right=br4) assert out4.shape == (10, 10), "overflow right shape" assert np.all(out4[:, :5] == 2), "overflow right 1" assert np.all(out4[:, 5:] == 1), "overflow right 2" @@ -279,9 +279,9 @@ def test_crop_img_lat_lon(): def test_latlon_to_uv(): import accessvis - assert accessvis.latlon_to_uv(0, 0) == (0.5, 0.5) - assert accessvis.latlon_to_uv(90, -180) == (0, 0) - assert accessvis.latlon_to_uv(-90, 180) == (1, 1) + assert accessvis.lonlat_to_uv(lat=0, lon=0) == (0.5, 0.5) + assert accessvis.lonlat_to_uv(lat=90, lon=-180) == (0, 0) + assert accessvis.lonlat_to_uv(lat=-90, lon=180) == (1, 1) def test_uv_to_pixel(): @@ -295,7 +295,7 @@ def test_latlon_to_pixel(): lon = 0.1 * 360 - 180 lat = (1 - 0.8) * 180 - 90 # same as test_uv_to_pixel - assert accessvis.latlon_to_pixel(lon=lon, lat=lat, width=1017, height=8256) == ( + assert accessvis.lonlat_to_pixel(lon=lon, lat=lat, width=1017, height=8256) == ( 101, 6604, ) @@ -456,8 +456,6 @@ def test_lonlat_grid_3D(): longitudes=np.array([0, 90, 180, 270]), altitude=0, ) - print(out3.shape) - print(out3) R = 6.371 exp3 = np.array( [ From a849a8bba0030826034b326bc17dfbe815896649 Mon Sep 17 00:00:00 2001 From: max Date: Fri, 20 Jun 2025 17:19:26 +1000 Subject: [PATCH 2/3] fixing lonlat for plot_region --- src/accessvis/earth.py | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/accessvis/earth.py b/src/accessvis/earth.py index ab85de6..9253629 100644 --- a/src/accessvis/earth.py +++ b/src/accessvis/earth.py @@ -810,6 +810,21 @@ def load_topography( return heights # * vertical_exaggeration +def lonlat_minmax(a, b): + """ + Takes two lon/lat tuples, + returns two tuples containing the smallest, largest of each lon/lat. + """ + + lon1, lat1 = a + lon2, lat2 = b + if lon1 > lon2: + lon1, lon2 = lon2, lon1 + if lat1 > lat2: + lat1, lat2 = lat2, lat1 + return (lon1, lat1), (lon2, lat2) + + def plot_region( lv=None, top_left=None, @@ -855,6 +870,10 @@ def plot_region( vertical_exaggeration: number Multiplier to topography/bathymetry height """ + + if top_left and bottom_right: + top_left, bottom_right = lonlat_minmax(top_left, bottom_right) + if lv is None: lv = lavavu.Viewer( border=False, @@ -899,9 +918,10 @@ def plot_region( height = height * MtoLL * vertical_exaggeration D = [height.shape[1], height.shape[0]] + sverts = np.zeros(shape=(height.shape[0], height.shape[1], 3)) - lat0, lon0 = top_left if top_left else [-90, 0] - lat1, lon1 = bottom_right if bottom_right else [90, 360] + lon0, lat0 = top_left if top_left else [-90, 0] + lon1, lat1 = bottom_right if bottom_right else [90, 360] # TODO: Support crossing zero in longitude # Will probably not support crossing poles xy = lv.grid2d(corners=((lon0, lat1), (lon1, lat0)), dims=D) From dbb21d272410d09a0def3f3cf57fb2eb1b65e63d Mon Sep 17 00:00:00 2001 From: max Date: Wed, 25 Jun 2025 14:20:13 +1000 Subject: [PATCH 3/3] documentation --- src/accessvis/earth.py | 231 +++++++++++++++++++++++++++++++++++------ 1 file changed, 202 insertions(+), 29 deletions(-) diff --git a/src/accessvis/earth.py b/src/accessvis/earth.py index 9253629..f2ef39d 100644 --- a/src/accessvis/earth.py +++ b/src/accessvis/earth.py @@ -90,7 +90,9 @@ def set_resolution(val): Parameters ---------- - val: Texture and geometry resolution + val: int + Texture and geometry resolution + 1=low ... 4=high """ settings.RES = val settings.TEXRES = pow(2, 10 + val) @@ -287,6 +289,17 @@ def split_tex(data, res, flipud=False, fliplr=False): """ Convert a texture image from equirectangular to a set of 6 cubemap faces (requires py360convert) + + Parameters + ---------- + data: np.ndarray + data to split into cube faces. + res: int + Resolution of each face + flipud: bool|str + if up/down should be flipped (or which faces should be flipped) + fliplr: bool|str + if left/right should be flipped (or which faces should be flipped) """ if isinstance(flipud, bool): # Can provide a list/string of faces to flip, or just a single bool to flip all @@ -312,6 +325,21 @@ def split_tex(data, res, flipud=False, fliplr=False): def draw_lonlat_grid(base_img, out_fn, lon=30, lat=30, linewidth=5, colour=0): """ Create lat/lon grid image over provided base image + + Parameters + ---------- + base_img: + Path to the input image + out_fn: + Path to the output image + lon: int + Distance between lines of longitude (degrees) + lat: int + Distance between lines of latitude (degrees) + linewidth: int + pixel width of the lines + colour: + Colour of the lines """ from PIL import Image @@ -476,6 +504,20 @@ def crop_img_lon_lat(img, top_left, bottom_right): """ Crop an equirectangular image (PIL or numpy array) based on corner coords Provide coords as lat/lon coords in decimal degrees + + Parameters + ---------- + img: np.ndarray or PIL.ImageFile.ImageFile + The image to be cropped (numpy or PIL). + top_left: tuple + (lon, lat) Top left coordinates of the crop box. + If top_left or bottom_right are None, it returns the mask of the full earth. + bottom_right: tuple + (lon, lat) Bottom right coordinates of the crop box. + If top_left or bottom_right are None, it returns the mask of the full earth. + Returns: + ---------- + np.ndarray or PIL.ImageFile.ImageFile """ a = lonlat_to_uv(*top_left) b = lonlat_to_uv(*bottom_right) @@ -760,7 +802,23 @@ def load_topography( """ Load topography from pre-saved equirectangular data, can be cropped for regional plots - TODO: document args + Parameters + ---------- + resolution: + Resolution of the topography in the Y direction. + defaults to settings.FULL_RES_Y (default 10800) + Note: res_x = 2 * res_y + subsample: int + Selects every Nth item from the mask. + e.g. subsample=2 takes every second lat and lon value. + top_left: tuple + (lon, lat) Top left coordinates of the crop box. + If top_left or bottom_right are None, it returns the mask of the full earth. + bottom_right: tuple + (lon, lat) Bottom right coordinates of the crop box. + If top_left or bottom_right are None, it returns the mask of the full earth. + bathymetry: + If the ocean floor should be shown. """ if resolution is None: resolution = settings.FULL_RES_Y @@ -812,8 +870,8 @@ def load_topography( def lonlat_minmax(a, b): """ - Takes two lon/lat tuples, - returns two tuples containing the smallest, largest of each lon/lat. + Takes two lon/lat tuples (lon, lat). + Returns two tuples containing the smallest, largest of each lon/lat. """ lon1, lat1 = a @@ -840,7 +898,6 @@ def plot_region( uniforms={}, shaders=None, background="black", - *args, **kwargs, ): """ @@ -850,16 +907,25 @@ def plot_region( Uses lat/lon as coordinate system, so no use for polar regions, scales heights to equivalent TODO: support using km as unit or other custom units instead with conversions from lat/lon - TODO: FIX LAT/LON ordering - we are using lon/lat in all interfaces from now on - TODO: cropbox: defaults to full earth if None passed, write a test for this - FINISH DOCUMENTING PARAMS + TODO: top_left/bottom_right: defaults to full earth if None passed, write a test for this + + TODO: Implement lighting, waves, shaders arguments (currently unused) + TODO: Implement bathy_exaggeration and topo_exaggeration Note: If you want to plot data, you should use `plot_region_data()` after calling this function. If you want to plot data in a region, but continue to display the entire 3D earth, you should use `earth_2d_plot()` instead. Parameters ---------- + lv: lavavu.Viewer + The viewer object to plot with + top_left: tuple + (lon, lat) Top left coordinates of the crop box. + bottom_right: tuple + (lon, lat) Bottom right coordinates of the crop box. + vertical_exaggeration: number + Multiplier to topography/bathymetry height texture: str Path to textures, face label and texres will be applied with .format(), eg: texture='path/{face}_mytexture_{texres}.png' @@ -869,6 +935,23 @@ def plot_region( Append this label to each face object created vertical_exaggeration: number Multiplier to topography/bathymetry height + when: datetime + Provide a datetime object to set the month for texture sets that vary over the year and time for + position of sun and rotation of earth when calculating sun light position + blendtex: bool + When the texture set has varying seasonal images, enabling this will blend between the current month and next + months texture to smoothly transition between them as the date changes, defaults to True + bathymetry: bool + If the ocean floor should be shown + name: str + The name of the surface + uniforms: dict + Provide a set of uniform variables, these can be used to pass data to a custom shader + background: str + Provide a background colour string, X11 colour name or hex RGB + kwargs: + Additional lavavu arguments. + Covers global props, object props and uniform values """ if top_left and bottom_right: @@ -984,7 +1067,6 @@ def plot_earth( shaders=None, background="black", hemisphere=None, - *args, **kwargs, ): """ @@ -992,15 +1074,10 @@ def plot_earth( and sets up seasonal texture blending and optionally, sun position, based on given or current date and time - TODO: FINISH DOCUMENTING PARAMS - Parameters ---------- - texture: str - Path to textures, face label and texres will be applied with .format(), eg: - texture="path/{face}_mytexture_{texres}.png" - with: texture.format(face="F", texres=settings.TEXRES) - to:"path/F_mytexture_1024.png" + lv: lavavu.Viewer + The viewer object to plot with radius: float Radius of the sphere, defaults to 6.371 Earth's approx radius in Mm vertical_exaggeration: number @@ -1010,6 +1087,11 @@ def plot_earth( bathy_exaggeration: number Multiplier to bathymetry height only texture: str + Path to textures, face label and texres will be applied with .format(), eg: + texture="path/{face}_mytexture_{texres}.png" + with: texture.format(face="F", texres=settings.TEXRES) + to:"path/F_mytexture_1024.png" + Texture set to use, "bluemarble" for the 2004 NASA satellite data, "relief" for a basic relief map or provide a custom set of textures using a filename template with the following variables, only face is required {face} (F/R/B/L/U/D) {month} (name of month, capitialised) {texres} (2048/4096/8192/16384) @@ -1046,6 +1128,9 @@ def plot_earth( "WE" = Prime meridian at centre (Africa/Europe) "E" = Eastern hemisphere - prime meridian to antimeridian (Indian ocean) "W" = Western hemisphere - antimeridian to prime meridian (Americas) + kwargs: + Additional lavavu arguments. + Covers global props, object props and uniform values """ if lv is None: lv = lavavu.Viewer( @@ -1283,6 +1368,17 @@ def earth_2d_plot( Assumes the first coordinate of the array is longitude and second is latitude If order is reversed, need to transpose before passing, eg: data.transpose('latitude', 'longitude') If omitted, will skip texturing + colourmap: str or list or object + If a string, provides the name of a matplotlib colourmap to use + If a list is passed, should contain RGB[A] colours as lists or tuples in range [0,1] or [0,255], + a matplotlib.colors.LinearSegmentedColormap will be created from the list + Otherwise, will assume is a valid matplotlib colormap already + opacitymap: bool or numpy.ndarray + Set to true to use values as an opacity map, top of range will be opaque, bottom transparent + Provide an array to use a different opacity map dataset + opacitypower: float + Power to apply to values when calculating opacity map, + eg: 0.5 equivalent to opacity = sqrt(value), 2 equivalent to opacity = value*value longitudes: ndarray or tuple List of longitude values or min/max longitude as a tuple latitudes: ndarray or tuple @@ -1293,7 +1389,8 @@ def earth_2d_plot( resolution: 2-tuple Resolution of the 2d grid, if provided will use corner points from data only and sample a regular grid between them of provided resolution - + kwargs: + These kwargs are passed through to lv.surface(). Returns ------- lavavu.Object: the drawing object created @@ -1345,7 +1442,25 @@ def earth_2d_plot( def update_earth_datetime( lv, when, name="", texture=None, sunlight=False, blendtex=True ): - # Update date/time for texture blending and lighting + """ + Update date/time for texture blending and lighting + + Parameters + ---------- + lv: lavavu.Viewer + The earth viewer object to plot with. + when: datetime.datetime + timestamp the earth should reflect. + name: str + the name of the surfaces the texture should be applied to. + texture: + "bluemarble" or path to the file containing the texture. + sunlight: bool + True means the light source should move to the correct position to reflect the given time. + blendtex: bool + If the bluemarble texture should be blended between months to approximate the given date. + """ + d = when.day - 1 m = when.month month = when.strftime("%B") @@ -1404,9 +1519,26 @@ def update_earth_datetime( def update_earth_texture( - lv, label, texture, flip=False, shared=True, name="", *args, **kwargs + lv, label, texture, flip=False, shared=True, name="", **kwargs ): - # Update texture values for a specific texture on earth model + """ + Update texture values for a specific texture on earth model + + Parameters + ---------- + lv: lavavu.Viewer + The earth viewer object to plot with. + label: str + flip: bool + Whether or not the textures should be flipped. + shared: bool + If the texture is shared between faces. + name: + The name of the object which the kwargs should be passed through to. + kwargs: + lavavu kwargs to pass through to the uniforms. + """ + if shared: # No need to update each object lv.texture(label, texture, flip) @@ -1425,9 +1557,22 @@ def update_earth_texture( # lv.render() # Required to render a frame which fixes texture glitch -def update_earth_values(lv, name="", flip=False, *args, **kwargs): - # Update uniform values on earth objects via passed kwargs +def update_earth_values(lv, name="", flip=False, **kwargs): + """ + Update uniform values on earth objects via passed kwargs + + Parameters + ---------- + lv: lavavu.Viewer + The earth viewer object to plot with. + name: str + The name of the object which the kwargs should be passed through to. + flip: bool + Whether or not the textures should be flipped. + kwargs: + Lavavu arguments to pass through to the textures. + """ # Replace texture data load shared texture afterwards for k in kwargs: if isinstance(kwargs[k], (np.ndarray, np.generic)) or k == "data": @@ -1838,7 +1983,6 @@ def array_to_rgba( def load_mask( res_y=None, masktype="watermask", subsample=1, top_left=None, bottom_right=None ): - # TODO: Document args """ Loads watermask/oceanmask @@ -1852,7 +1996,24 @@ def load_mask( https://neo.gsfc.nasa.gov/archive/bluemarble/bmng/landmask_new/ https://neo.gsfc.nasa.gov/archive/bluemarble/bmng/landmask/world.watermask.21600x21600.A1.png - masktype = "oceanmask" / "watermask" + Parameters + ---------- + res_y: + defaults to settings.FULL_RES_Y (default 10800) + Resolution of the oceanmask png file. + res_x = 2 * res_y + masktype: str + "relief": lower quality + "oceanmask" or "watermask": better quality + subsample: int + Selects every Nth item from the mask. + e.g. subsample=2 takes every second lat and lon value. + top_left: tuple + (lon, lat) Top left coordinates of the crop box. + If top_left or bottom_right are None, it returns the mask of the full earth. + bottom_right: tuple + (lon, lat) Bottom right coordinates of the crop box. + If top_left or bottom_right are None, it returns the mask of the full earth. """ if res_y is None: res_y = settings.FULL_RES_Y @@ -1942,9 +2103,8 @@ def process_bluemarble(when=None, overwrite=False, redownload=False, blendtex=Tr """ Download and process NASA Blue Marble next gen imagery - month: int - Month [1-12] - If omitted will gather and process all 12 month images + when: datetime.datetime + If None, it defaults to the current time. overwrite: bool Always re-process from source images overwriting any existing redownload: bool @@ -2089,6 +2249,9 @@ def process_landmask(texture, overwrite=False, redownload=False): """ Download and process ocean/water mask imagery + texture: str + "relief": lower quality but smaller. + Otherwise it uses both watermask and oceanmask to generate the mask overwrite: bool Always re-process from source images overwriting any existing redownload: bool @@ -2250,11 +2413,11 @@ def process_landmask(texture, overwrite=False, redownload=False): mimg.save(r_fn) -def process_gebco(cubemap, resolution, overwrite=False, redownload=False): +def process_gebco(cubemap, resolution, overwrite=False): """ # Full res GEBCO .nc grid - This function generates cubemap sections at the desired resolution from the full res GEBCO dataset + This function generates cubemap sections at the desired resolution from the full res GEBCO dataset. https://www.bodc.ac.uk/data/hosted_data_systems/gebco_gridded_bathymetry_data/ @@ -2262,7 +2425,17 @@ def process_gebco(cubemap, resolution, overwrite=False, redownload=False): - See https://download.gebco.net/ - NC version: https://www.bodc.ac.uk/data/open_download/gebco/gebco_2020/zip/ - Sub-ice topo version: https://www.bodc.ac.uk/data/open_download/gebco/gebco_2023_sub_ice_topo/zip/ + + Parameters + ---------- + cubemap: bool + Whether to download cubemap or equirectangular data. + resolution: int + Resolution of each cube face (e.g. 1024). + overwrite: bool + True means the files should be redownloaded. """ + if cubemap: fspec = f"{settings.DATA_PATH}/gebco/gebco_cubemap_{resolution}.npz" if not overwrite and len(glob.glob(fspec)) == 1: