orthographic/xp_scenery.py

424 lines
17 KiB
Python

# -------------------------------------------------------------------
# ORTHOGRAPHIC
# Your personal aerial satellite. Always on. At any altitude.*
# Developed by MarStrMind
# License: Open Software License 3.0
# Up to date version always on marstr.online
# -------------------------------------------------------------------
# xp_scenery.py
# This class builds an elevation mesh from provided DEM data, and
# generates all required files that are needed to provide a usable
# scenery package for X-Plane.
# -------------------------------------------------------------------
import os
import math
import urllib.request
import numpy
from defines import *
from log import *
from PIL import Image, ImageFilter, ImageEnhance
class mstr_xp_scenery:
# Set required variables
def __init__(self, lat, lng, mlat, mlng, vstep, latlngfld):
mstr_msg("xp_scenery", "[X-Plane] Scenery generator instantiated")
self._lat = lat
self._lng = lng
self._mlat = mlat
self._mlng = mlng
self._vstep = vstep
self._latlngfld = latlngfld
self._demfn = self.build_dem_filename()
self._dsfstring = ""
self._demdata = None # To be populated when the mesh is built
self._demcoord = None # Also to be populated when mesh is built
self._waterdata = [] # So that we know where to implement water
#self.load_water_data()
# Build the correct file name for the elevation model
def build_dem_filename(self, xes=False):
fn = ""
if self._lat > 0:
fn = fn + "N"
else:
fn = fn + "S"
if abs(self._lat) < 10: fn = fn + "0" + str(self._lat)
if abs(self._lat) >= 10 and abs(self._lat) <= 90: fn = fn + str(self._lat)
if self._lng > 0:
fn = fn + "E"
else:
fn = fn + "W"
if abs(self._lng) < 10: fn = fn + "00" + str(self._lng)
if abs(self._lng) >= 10 and abs(self._lng) <= 99: fn = fn + "0" + str(self._lng)
if abs(self._lng) >= 100 : fn = fn + str(self._lng)
if xes == False:
fn = fn + ".hgt"
if xes == True:
fn = fn + ".xes"
mstr_msg("xp_scenery", "[X-Plane] DEM file name constructed: " + fn)
return fn
# Load the water data before we generate the mesh
def load_water_data(self):
fn = mstr_datafolder + "z_orthographic/data/" + self._latlngfld + "/wtrfile"
with open(fn) as file:
for line in file:
ln = line.replace(" ", "_")
ln = ln.replace("\n", "")
ln = ln.replace("\r", "")
self._waterdata.append(ln)
# Check if ortho has water
def does_ortho_have_water(self, ortho):
wtr = False
if ortho in self._waterdata: wtr = True
return wtr
# Build the DSF for the ortho photo overlay
def build_and_convert_dsf(self):
end = self.find_earthnavdata_number()
llf = self.xplane_latlng_folder(end)
meshtxt = mstr_datafolder + "_cache/mesh_"+self._latlngfld+".txt"
cmd = mstr_xp_dsftool + " --text2dsf " + meshtxt + " '" + mstr_datafolder + "z_orthographic/Earth nav data/" + llf + "/" + self._latlngfld + ".dsf'"
os.system(cmd)
# Find exact with of longitude
def find_width_of_longitude(self, lat):
dm = math.cos(math.radians(lat)) * 111.321 # <- 1 deg width at equator in km
return round(dm * 1000, 3)
# Find the next "by-ten" numbers for the current latitude and longitude
def find_earthnavdata_number(self):
earthnavdata = []
lat = abs(int(self._lat / 10) * 10)
lng = abs(int(self._lng / 10) * 10)
earthnavdata.append(lat)
earthnavdata.append(lng)
return earthnavdata
# Construct an X-Plane compatible folder name for latitude and longitude
def xplane_latlng_folder(self, numbers):
fstr = ""
if numbers[0] >= 0: fstr = "+"
if numbers[0] < 0: fstr = "-"
if abs(numbers[0]) < 10: fstr = fstr + "0" + str(numbers[0])
if abs(numbers[0]) >= 10 and numbers[0] <= 90: fstr = fstr + str(numbers[0])
if numbers[1] >= 0: fstr = fstr + "+"
if numbers[1] < 0: fstr = fstr + "-"
if abs(numbers[1]) < 10: fstr = fstr + "00" + str(numbers[1])
if abs(numbers[1]) >= 10 and numbers[0] <= 99: fstr = fstr + "0" + str(numbers[1])
if abs(numbers[1]) >= 100 : fstr = fstr + str(numbers[1])
return fstr
# Acquires the elevation data as needed -
# you can either acquire the older STRM set with lower resolution,
# or a modern LIDAR scan with higher resolution. However, the
# LIDAR files are much larger and generates a more taxing mesh.
# If you are testing, acquire the low resolution file.
# Both versions are coming from my repository at marstr.online .
def acquire_elevation_data(self):
mstr_msg("xp_scenery", "[X-Plane] Acquiring DEM model data")
url = "https://marstr.online/dem/"
url = url + self._demfn
urllib.request.urlretrieve(url, mstr_datafolder + "_cache/" + self._demfn)
mstr_msg("xp_scenery", "[X-Plane] DEM data acquired")
# Download the X-Plane definition file from my server
def acquire_xes_data(self):
mstr_msg("xp_scenery", "[X-Plane] Acquiring XES file")
url = "https://marstr.online/xes/"
xesfn = self.build_dem_filename(True)
xp_folder = self.xplane_latlng_folder([self._lat, self._lng])
url = url + xp_folder + ".xes"
urllib.request.urlretrieve(url, mstr_datafolder + "_cache/" + xesfn)
mstr_msg("xp_scenery", "[X-Plane] XES data acquired")
# This generates all .ter files
def build_ter_files(self):
mstr_msg("xp_scenery", "[X-Plane] Generating and writing terrain (.ter) files")
cur_lat = self._lat
cur_lng = self._lng
xp_folder = self.xplane_latlng_folder([self._lat, self._lng])
for lat in range(1, self._mlat+1):
for lng in range(1, self._mlng+1):
wdt = self.find_width_of_longitude(cur_lat)
dmt = wdt * mstr_zl_18
cnt_x = cur_lat + (self._vstep/2)
cnt_y = cur_lng + (mstr_zl_18/2)
terstr = ""
terstr = terstr + "A\r\n"
terstr = terstr + "800\r\n"
terstr = terstr + "TERRAIN\r\n"
terstr = terstr + "\r\n"
terstr = terstr + "LOAD_CENTER " + str(cnt_x) + " " + str(cnt_y) + " " + str(dmt) + " 2048\r\n"
terstr = terstr + "BASE_TEX_NOWRAP ../../orthos/" + self._latlngfld + "/" + str(lat)+"_"+str(lng)+".dds\r\n"
if mstr_xp_scn_normalmaps:
terstr = terstr + "NORMAL_TEX 1.0 ../../normals/" + self._latlngfld + "/" + str(lat)+"_"+str(lng)+".png\r\n"
terfln = mstr_datafolder + "z_orthographic/terrain/" + self._latlngfld + "/" + str(lat)+"_"+str(lng)+".ter"
with open(terfln, 'w') as textfile:
textfile.write(terstr)
cur_lng = round(cur_lng + mstr_zl_18, 6)
cur_lng = self._lng
cur_lat = round(cur_lat + self._vstep, 6)
mstr_msg("xp_scenery", "[X-Plane] Terrain files written")
# This generates the entire terrain mesh
def generate_terrain_mesh(self):
# Get the DEM model file name, and acquire important info about the data
meshfn = mstr_datafolder + "_cache/" + self.build_dem_filename()
siz = os.path.getsize(meshfn)
dim = int(math.sqrt(siz/2))
assert dim*dim*2 == siz, 'Invalid file size'
self._demdata = numpy.fromfile(meshfn, numpy.dtype('>i2'), dim*dim).reshape((dim, dim))
self._demdata = self._demdata[::-1] # Invert order so that we can start from bottom left
# We want to achieve perfect stepping for each data point in the DEM.
demstep = round( 1 / len(self._demdata), 6)
# Generate an array which contains only the coordinates
self._demcoord = []
for r in range(0, len(self._demdata)):
row = []
for c in range(0, len(self._demdata)):
lat = round(self._lat + r * demstep, 6)
lng = round(self._lng + c * demstep, 6)
crd = [ lat, lng, self._demdata[r][c]]
#crd = [ lat, lng ]
row.append(crd)
self._demcoord.append(row)
mstr_msg("xp_scenery", "[X-Plane] Populating DSF information file")
# The complete string to write into the DSF txt file
dsf_str = ""
dsf_str = dsf_str + "PROPERTY sim/west " + str(int(self._lng)) + "\r\n"
dsf_str = dsf_str + "PROPERTY sim/east " + str((int(self._lng) + 1)) + "\r\n"
dsf_str = dsf_str + "PROPERTY sim/south " + str(int(self._lat)) + "\r\n"
dsf_str = dsf_str + "PROPERTY sim/north " + str((int(self._lat) + 1)) + "\r\n"
dsf_str = dsf_str + "PROPERTY sim/require_object 0/6\r\n"
dsf_str = dsf_str + "PROPERTY planet earth\r\n"
dsf_str = dsf_str + "PROPERTY sim/creation_agent Orthographic\r\n"
#dsf_str = dsf_str + "TERRAIN_DEF terrain_Water\r\n"
# The file to be converted into DSF later
meshtxt = mstr_datafolder + "_cache/mesh_"+self._latlngfld+".txt"
with open(meshtxt, 'w') as textfile:
textfile.write(dsf_str)
dsf_str = ""
# Orthos
for lat in range(1, self._mlat+1):
for lng in range(1, self._mlng+1):
# Write the line only if an ortho exists of course.
ddsf = mstr_datafolder + "z_orthographic/orthos/" + self._latlngfld + "/" + str(lat) + "_" + str(lng) + ".dds"
if os.path.isfile(ddsf):
dsf_str = dsf_str + "TERRAIN_DEF terrain/" + self._latlngfld + "/" + str(lat) + "_" + str(lng) + ".ter\r\n"
with open(meshtxt, 'a') as textfile:
textfile.write(dsf_str)
# OK. So. Let's build the mesh.
# Current patch
curpatch = 0
for lat in range(1, self._mlat+1):
for lng in range(1, self._mlng+1):
# Create the patch only if the matching ortho exists.
# This way we make sure that we hit the same order as the .ter files.
# We can also detect which lat and lng coord we are on.
ddsf = mstr_datafolder + "z_orthographic/orthos/" + self._latlngfld + "/" + str(lat) + "_" + str(lng) + ".dds"
if os.path.isfile(ddsf):
# Base coords for this ortho
base_lat = self._lat + ((lat-1) * self._vstep)
base_lng = self._lng + ((lng-1) * mstr_zl_18)
# Begin a new patch
mstr_msg("xp_scenery", "[X-Plane] Processing ortho patch " + str(curpatch))
with open(meshtxt, 'a') as textfile:
textfile.write("BEGIN_PATCH " + str(curpatch) + " 0.000000 -1.000000 1 7\r\n")
# Step for each ortho vertex
odiv = 4
latstep = self._vstep/odiv
lngstep = mstr_zl_18 /odiv
uv_step = 1 / odiv
# Generate the ortho tile
for y in range(0,odiv):
for x in range(0,odiv):
# Coordinates
lat_b = round(base_lat + (y*latstep), 6)
lat_t = round(base_lat + ((y+1)*latstep), 6)
lng_l = round(base_lng + (x*lngstep), 6)
lng_r = round(base_lng + ((x+1)*lngstep), 6)
# Minimal adjustment
if x == 0:
lng_l = base_lng
if y == 0:
lat_b = base_lat
if y == 3:
lat_t = base_lat + self._vstep
if x == 3:
lng_r = base_lng + mstr_zl_18
# Corrections, just in case
if lat_b > self._lat + 1: lat_b = self._lat+1
if lat_t > self._lat + 1: lat_t = self._lat+1
if lng_l > self._lng + 1: lng_l = self._lng+1
if lng_r > self._lng + 1: lng_r = self._lng+1
# Height indexes
hgt_bl_idx = self.find_height_for_coord([lat_b, lng_l])
hgt_br_idx = self.find_height_for_coord([lat_b, lng_r])
hgt_tr_idx = self.find_height_for_coord([lat_t, lng_r])
hgt_tl_idx = self.find_height_for_coord([lat_t, lng_l])
hgt_bl = round(self._demcoord[ hgt_bl_idx[0] ][ hgt_bl_idx[1] ][2], 6)
hgt_br = round(self._demcoord[ hgt_br_idx[0] ][ hgt_br_idx[1] ][2], 6)
hgt_tr = round(self._demcoord[ hgt_tr_idx[0] ][ hgt_tr_idx[1] ][2], 6)
hgt_tl = round(self._demcoord[ hgt_tl_idx[0] ][ hgt_tl_idx[1] ][2], 6)
# Coords of triangle vertices
# 0 - Longitude
# 1 - Latitude
# 2 - Height in m
t1_v1 = [ lng_r, lat_b, hgt_br ]
t1_v2 = [ lng_l, lat_t, hgt_tl ]
t1_v3 = [ lng_r, lat_t, hgt_tr ]
t2_v1 = [ lng_l, lat_t, hgt_tl ]
t2_v2 = [ lng_r, lat_b, hgt_br ]
t2_v3 = [ lng_l, lat_b, hgt_bl ]
# Write down the two triangles
t_str = ""
t_str = t_str + "BEGIN_PRIMITIVE 0\r\n"
t_str = t_str + "PATCH_VERTEX " + str(t1_v1[0]) + " " + str(t1_v1[1]) + " " + str(t1_v1[2]) + " 0.000015 0.000015 " + str((x+1) * uv_step) + " " + str(y*uv_step) + "\r\n"
t_str = t_str + "PATCH_VERTEX " + str(t1_v2[0]) + " " + str(t1_v2[1]) + " " + str(t1_v2[2]) + " 0.000015 0.000015 " + str(x * uv_step) + " " + str((y+1)*uv_step) + "\r\n"
t_str = t_str + "PATCH_VERTEX " + str(t1_v3[0]) + " " + str(t1_v3[1]) + " " + str(t1_v3[2]) + " 0.000015 0.000015 " + str((x+1) * uv_step) + " " + str((y+1)*uv_step) + "\r\n"
t_str = t_str + "END_PRIMITIVE 0\r\n"
t_str = t_str + "BEGIN_PRIMITIVE 0\r\n"
t_str = t_str + "PATCH_VERTEX " + str(t2_v1[0]) + " " + str(t2_v1[1]) + " " + str(t2_v1[2]) + " 0.000015 0.000015 " + str(x * uv_step) + " " + str((y+1)*uv_step) + "\r\n"
t_str = t_str + "PATCH_VERTEX " + str(t2_v2[0]) + " " + str(t2_v2[1]) + " " + str(t2_v2[2]) + " 0.000015 0.000015 " + str((x+1) * uv_step) + " " + str(y*uv_step) + "\r\n"
t_str = t_str + "PATCH_VERTEX " + str(t2_v3[0]) + " " + str(t2_v3[1]) + " " + str(t2_v3[2]) + " 0.000015 0.000015 " + str(x * uv_step) + " " + str(y*uv_step) + "\r\n"
t_str = t_str + "END_PRIMITIVE 0\r\n"
# Send to the file
with open(meshtxt, 'a') as textfile:
textfile.write(t_str)
t_str = ""
# End this patch
with open(meshtxt, 'a') as textfile:
textfile.write("END PATCH\r\n")
# Increase patch number
curpatch = curpatch + 1
# Find the next best matching height for a point
def find_height_for_coord(self, coord):
idx = [0,0]
dst = 99999
ste = self.find_height_scan_start_end_points(coord)
for r in range(ste[0], ste[1]+1):
for d in range(ste[2], ste[3]+1):
dist = math.dist(coord, [self._demcoord[r][d][0], self._demcoord[r][d][1]])
if dist < dst:
dst = dist
idx = [r,d]
return idx
# Find the starting and end points to scan for heights in the DEM grid
def find_height_scan_start_end_points(self, stc):
startend = [0,0,0,0]
stp = 1 / len(self._demdata)
# Bottom
lt = self._lat
while lt < stc[0]:
lt = lt + stp
startend[0] = startend[0] + 1
# Top
lt = self._lat
while lt < stc[0]+self._vstep:
lt = lt+stp
startend[1] = startend[1] + 1
# Left
ln = self._lng
while ln < stc[1]:
ln = ln + stp
startend[2] = startend[2] + 1
# Right
ln = self._lng
while ln < stc[1]+mstr_zl_18:
ln = ln + stp
startend[3] = startend[3] + 1
# Make sure we have everything
startend[0] = startend[0]-1
startend[1] = startend[1]+1
startend[2] = startend[2]-1
startend[3] = startend[3]+1
# Some corrections
if startend[0] < 0: startend[0] = 0
if startend[1] > len(self._demdata)-1: startend[1] = startend[1] = len(self._demdata)-1
if startend[2] < 0: startend[2] = 0
if startend[3] > len(self._demdata)-1: startend[3] = startend[3] = len(self._demdata)-1
return startend