Skip to content

Commit b5fe814

Browse files
OKaluzaOKaluza
committed
Update earth model
1 parent 2b5e07c commit b5fe814

File tree

4 files changed

+742
-39
lines changed

4 files changed

+742
-39
lines changed

‎Earth/Earth-Model.ipynb

+525-38
Large diffs are not rendered by default.

‎Earth/SphereMesh.ipynb

+143
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "b959fb73-7ccf-4fa0-a5b7-d23b26a3e634",
6+
"metadata": {},
7+
"source": [
8+
"# Sphere Mesh\n",
9+
"\n",
10+
"Generate a spherical mesh for use in plotting 3D globe data\n",
11+
"\n",
12+
"Uses specified radius, add an offset to plot data just about topography etc"
13+
]
14+
},
15+
{
16+
"cell_type": "code",
17+
"execution_count": null,
18+
"id": "c086c0ed-e736-4fb5-917f-312507f6ef60",
19+
"metadata": {},
20+
"outputs": [],
21+
"source": [
22+
"#Earth's radius ~6371km - we'll use 1000's of km units (Mm)\n",
23+
"radius = 6371 * 1e-3 #Radius in Mm\n",
24+
"\n",
25+
"#Radius + small offset for stratosphere\n",
26+
"s_radius = radius * 1.016\n",
27+
"\n",
28+
"#Sphere mesh quality (higher = more triangles)\n",
29+
"QUALITY=256"
30+
]
31+
},
32+
{
33+
"cell_type": "code",
34+
"execution_count": null,
35+
"id": "e63255a8-e12b-4316-b18d-85bdf46940b3",
36+
"metadata": {},
37+
"outputs": [],
38+
"source": [
39+
"import lavavu\n",
40+
"lv = lavavu.Viewer(border=False, axis=False, resolution=[1280,720], background='black')"
41+
]
42+
},
43+
{
44+
"cell_type": "code",
45+
"execution_count": null,
46+
"id": "e695a61a-4197-4abd-bd91-f78e7d281559",
47+
"metadata": {},
48+
"outputs": [],
49+
"source": [
50+
"tris0 = lv.spheres(\"sphere\", scaling=s_radius, segments=QUALITY, colour=\"grey\", vertices=[0,0,0], fliptexture=False)\n",
51+
"tris0['rotate'] = [0,-90,0] #This rotates the sphere coords to align with [0,360] longitude texture\n",
52+
"tris0['texture'] = 'blank.png' #Need an initial texture or texcoords will not be generated\n",
53+
"tris0['renderer'] = 'sortedtriangles'\n",
54+
"lv.render()\n",
55+
"\n",
56+
"#Generate and extract sphere vertices, texcoords etc\n",
57+
"lv.bake(7)\n",
58+
"#lv['cullface'] = False #Must disable this for the ozone plot\n",
59+
"#tris0[\"rotate\"] = [0,0,0]\n",
60+
"#tris0[\"alpha\"] = 0.6\n",
61+
"lv.display((400,400))"
62+
]
63+
},
64+
{
65+
"cell_type": "code",
66+
"execution_count": null,
67+
"id": "f39d7030-cc23-4783-8e9b-47cc5a77c794",
68+
"metadata": {},
69+
"outputs": [],
70+
"source": [
71+
"tris0.data[0]\n",
72+
"#DrawData(\"triangles\") ==> {'vertices': (65792, 3), 'normals': (65792, 3), 'indices': (196608,), 'texcoords': (65792, 2)}"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": null,
78+
"id": "e52c2419-e3e6-4349-be13-0f84cba06339",
79+
"metadata": {},
80+
"outputs": [],
81+
"source": [
82+
"sdata = {}\n",
83+
"element = tris0.data[0]\n",
84+
"keys = element.available.keys()\n",
85+
"for k in keys:\n",
86+
" print(k)\n",
87+
" sdata[k] = tris0.data[k][0]\n",
88+
"print(sdata)"
89+
]
90+
},
91+
{
92+
"cell_type": "code",
93+
"execution_count": null,
94+
"id": "85d8ff63-7e99-4dad-947c-cad318efa25f",
95+
"metadata": {},
96+
"outputs": [],
97+
"source": [
98+
"#Save compressed vertex data\n",
99+
"import numpy as np\n",
100+
"np.savez_compressed(f'Sphere_{s_radius:.4f}', **sdata)"
101+
]
102+
},
103+
{
104+
"cell_type": "code",
105+
"execution_count": null,
106+
"id": "718b2386-1daf-4865-bde9-1b57fa25570e",
107+
"metadata": {},
108+
"outputs": [],
109+
"source": [
110+
"!ls -ltrh *.npz"
111+
]
112+
},
113+
{
114+
"cell_type": "code",
115+
"execution_count": null,
116+
"id": "64718aa0-da17-424f-b579-a98d1588d378",
117+
"metadata": {},
118+
"outputs": [],
119+
"source": []
120+
}
121+
],
122+
"metadata": {
123+
"kernelspec": {
124+
"display_name": "Python 3 (ipykernel)",
125+
"language": "python",
126+
"name": "python3"
127+
},
128+
"language_info": {
129+
"codemirror_mode": {
130+
"name": "ipython",
131+
"version": 3
132+
},
133+
"file_extension": ".py",
134+
"mimetype": "text/x-python",
135+
"name": "python",
136+
"nbconvert_exporter": "python",
137+
"pygments_lexer": "ipython3",
138+
"version": "3.10.12"
139+
}
140+
},
141+
"nbformat": 4,
142+
"nbformat_minor": 5
143+
}

‎Earth/blank.png

1.02 KB
Loading

‎accessvis.py

+74-1
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,9 @@
22
# #!pip install py360convert
33
import numpy as np
44
import py360convert
5+
from PIL import Image
6+
import os
7+
Image.MAX_IMAGE_PIXELS = 1061683200
58

69
def latlon_to_3D(lat, lon, alt=0):
710
"""
@@ -40,6 +43,10 @@ def latlon_to_3D_true(lat, lon, alt=0, flattening=1.0/298.257223563):
4043
#return (x, y, z)
4144

4245
def split_tex(data, res, flip=[]):
46+
"""
47+
Convert a texture image from equirectangular to a set of 6 cubemap faces
48+
(requires py360convert)
49+
"""
4350
if len(data.shape) == 2:
4451
data = data.reshape(data.shape[0],data.shape[1],1)
4552
channels = data.shape[2]
@@ -55,7 +62,9 @@ def split_tex(data, res, flip=[]):
5562
return tiles
5663

5764
def draw_latlon_grid(base_img, out_fn, lat=30, lon=30, linewidth=5, colour=0):
58-
#Create lat/lon grid image over provided base image
65+
"""
66+
Create lat/lon grid image over provided base image
67+
"""
5968
from PIL import Image
6069

6170
#Open base image
@@ -85,3 +94,67 @@ def draw_latlon_grid(base_img, out_fn, lat=30, lon=30, linewidth=5, colour=0):
8594
display(image)
8695
image.save(out_fn)
8796

97+
def latlon_to_uv(lat, lon):
98+
"""
99+
Convert a decimal longitude, latitude coordinate
100+
to a tex coord in an equirectangular image
101+
"""
102+
#X/u E-W Longitude - [-180,180]
103+
u = 0.5 + lon/360.0
104+
105+
#Y/v N-S Latitude - [-90,90]
106+
v = 0.5 - lat/180.0
107+
108+
return u,v
109+
110+
def uv_to_pixel(u, v, width, height):
111+
"""
112+
Convert tex coord [0,1]
113+
to a raster image pixel coordinate for given width/height
114+
"""
115+
return int(u * width), int(v * height)
116+
117+
def latlon_to_pixel(lat, lon, width, height):
118+
"""
119+
Convert a decimal latitude/longitude coordinate
120+
to a raster image pixel coordinate for given width/height
121+
"""
122+
u, v = latlon_to_uv(lat, lon)
123+
return uv_to_pixel(u, v, width, height)
124+
125+
def crop_img_uv(img, top_left, bottom_right):
126+
"""
127+
Crop an image (PIL or numpy array) based on corner coords
128+
Provide coords as texture coords in [0,1]
129+
"""
130+
u0 = top_left[0]
131+
u1 = bottom_right[0]
132+
v0 = top_left[1]
133+
v1 = bottom_right[1]
134+
#Swap coords if order incorrect
135+
if u0 > u1:
136+
u0, u1 = u1, u0
137+
if v0 > v1:
138+
v0, v1 = v1, v0
139+
#Supports numpy array or PIL image
140+
if isinstance(img, np.ndarray):
141+
#Assumes [lat][lon]
142+
lat = int(v0*img.shape[0]), int(v1*img.shape[0])
143+
lon = int(u0*img.shape[1]), int(u1*img.shape[1])
144+
print(lat, lon)
145+
return img[lat[0] : lat[1], lon[0] : lon[1]]
146+
elif hasattr(img, 'crop'):
147+
area = (int(u0*img.size[0]), int(v0*img.size[1]), int(u1*img.size[0]), int(v1*img.size[1]))
148+
return img.crop(area)
149+
else:
150+
print("Unknown type: ", type(img))
151+
152+
def crop_img_lat_lon(img, top_left, bottom_right):
153+
"""
154+
Crop an equirectangular image (PIL or numpy array) based on corner coords
155+
Provide coords as lat/lon coords in decimal degrees
156+
"""
157+
a = latlon_to_uv(*top_left)
158+
b = latlon_to_uv(*bottom_right)
159+
return crop_img_uv(img, a, b)
160+

0 commit comments

Comments
 (0)