maxlatlng = [ mlat, mlng ]\r
\r
while grid_lat <= maxlatlng[0]:\r
- # Reset these two\r
- bb_lat = self._lat + ((grid_lat-1)*self._vstep)\r
- bb_lng = self._long + ((grid_lng-1)*mstr_zl_18)\r
- bb_lat_edge = self._lat + ((grid_lat-1)*self._vstep) + self._vstep\r
- bb_lng_edge = self._long + ((grid_lng-1)*mstr_zl_18) + mstr_zl_18\r
+ ddsf = mstr_datafolder + "z_orthographic/orthos/" + self._latlngfld + "/" + str(grid_lat) + "_" + str(grid_lng) + ".dds"\r
+ if os.path.isfile(ddsf) == False:\r
+ # Reset these two\r
+ bb_lat = self._lat + ((grid_lat-1)*self._vstep)\r
+ bb_lng = self._long + ((grid_lng-1)*mstr_zl_18)\r
+ bb_lat_edge = self._lat + ((grid_lat-1)*self._vstep) + self._vstep\r
+ bb_lng_edge = self._long + ((grid_lng-1)*mstr_zl_18) + mstr_zl_18\r
\r
- osmxml = mstr_osmxml()\r
- osmxml.adjust_bbox(bb_lat, bb_lng, bb_lat_edge, bb_lng_edge)\r
- osmxml.acquire_osm(grid_lat, grid_lng)\r
+ osmxml = mstr_osmxml()\r
+ osmxml.adjust_bbox(bb_lat, bb_lng, bb_lat_edge, bb_lng_edge)\r
+ osmxml.acquire_osm(grid_lat, grid_lng)\r
\r
- # Let the user know\r
- mstr_msg("orthographic", "Generating orthophoto " + str(grid_lat) + "-" + str(grid_lng))\r
+ # Let the user know\r
+ mstr_msg("orthographic", "Generating missing orthophoto " + str(grid_lat) + "-" + str(grid_lng))\r
\r
- # Check for work to be done\r
- layers = self.determineLayerWork(osmxml)\r
+ # Check for work to be done\r
+ layers = self.determineLayerWork(osmxml)\r
\r
- # We need to walk through the array of layers,\r
- # in their z-order.\r
- # For each layer, we will generate the mask, the layer image\r
- # itself, and finally, compose the ortho photo.\r
- mstr_msg("orthographic", "Beginning generation of layers")\r
+ # We need to walk through the array of layers,\r
+ # in their z-order.\r
+ # For each layer, we will generate the mask, the layer image\r
+ # itself, and finally, compose the ortho photo.\r
+ mstr_msg("orthographic", "Beginning generation of layers")\r
\r
- # In here we store the layers\r
- photolayers = []\r
+ # In here we store the layers\r
+ photolayers = []\r
\r
- # The masks are handed to layergen in sequence. The layers are then\r
- # in turn handed to photogen.\r
+ # The masks are handed to layergen in sequence. The layers are then\r
+ # in turn handed to photogen.\r
\r
- curlyr = 1\r
- for layer in layers:\r
- # Let the user know\r
- mstr_msg("orthographic", "Processing layer " + str(curlyr) + " of " + str(len(layers)))\r
-\r
- # Generate the mask\r
- mg = mstr_maskgen( [self._lat, grid_lat, self._long, grid_lng], self._vstep, layer[0], layer[1], layer[2])\r
- if layer[0] == "building":\r
- mg.set_tile_width(self._findWidthOfLongitude(bb_lat))\r
- mg.set_latlng_numbers(self._lat, grid_lat, self._long, grid_lng)\r
- mask = mg._build_mask(osmxml)\r
- \r
- # Generate the layer\r
- lg = mstr_layergen(layer[0], layer[1], self._lat, grid_lat, self._long, grid_lng, layer[2])\r
- lg.set_max_latlng_tile(maxlatlng)\r
- lg.set_latlng_folder(self._latlngfld)\r
- #lg.open_db()\r
- lg.open_tile_info()\r
- photolayers.append(lg.genlayer(mask, osmxml))\r
- curlyr = curlyr+1\r
- mstr_msg("orthographic", "All layers created")\r
-\r
- # We should have all layers now.\r
- # Snap a photo with our satellite :)\r
- mstr_msg("orthographic", "Generating ortho photo")\r
- pg = mstr_photogen(self._lat, self._long, grid_lat, grid_lng, maxlatlng[0], maxlatlng[1])\r
- pg.genphoto(photolayers)\r
- mstr_msg("orthographic", " -- Ortho photo generated -- ")\r
- print("")\r
- print("")\r
+ curlyr = 1\r
+ for layer in layers:\r
+ # Let the user know\r
+ mstr_msg("orthographic", "Processing layer " + str(curlyr) + " of " + str(len(layers)))\r
+\r
+ # Generate the mask\r
+ mg = mstr_maskgen( [self._lat, grid_lat, self._long, grid_lng], self._vstep, layer[0], layer[1], layer[2])\r
+ if layer[0] == "building":\r
+ mg.set_tile_width(self._findWidthOfLongitude(bb_lat))\r
+ mg.set_latlng_numbers(self._lat, grid_lat, self._long, grid_lng)\r
+ mask = mg._build_mask(osmxml)\r
+ \r
+ # Generate the layer\r
+ lg = mstr_layergen(layer[0], layer[1], self._lat, grid_lat, self._long, grid_lng, layer[2])\r
+ lg.set_max_latlng_tile(maxlatlng)\r
+ lg.set_latlng_folder(self._latlngfld)\r
+ #lg.open_db()\r
+ lg.open_tile_info()\r
+ photolayers.append(lg.genlayer(mask, osmxml))\r
+ curlyr = curlyr+1\r
+ mstr_msg("orthographic", "All layers created")\r
+\r
+ # We should have all layers now.\r
+ # Snap a photo with our satellite :)\r
+ mstr_msg("orthographic", "Generating ortho photo")\r
+ pg = mstr_photogen(self._lat, self._long, grid_lat, grid_lng, maxlatlng[0], maxlatlng[1])\r
+ pg.genphoto(photolayers)\r
+ mstr_msg("orthographic", " -- Ortho photo generated -- ")\r
+ print("")\r
+ print("")\r
\r
# Perform adjustment of grid position\r
n_lng = grid_lng + step\r
\r
\r
\r
+ # Generates X-Plane 11/12 scenery with\r
+ # - the finished orthos\r
+ # - a current LIDAR scan of the terrain\r
+ def generate_xp_scenery(self):\r
+ mstr_msg("orthographic", "[X-Plane] Generation of scenery started")\r
\r
+ # This call appears quite often... surely this can be done better\r
+ mlat = 1\r
+ mlng = 1\r
+ bb_lat = self._lat\r
+ bb_lng = self._long\r
+ bb_lat_edge = self._lat+self._vstep\r
+ bb_lng_edge = self._long+mstr_zl_18\r
+ while bb_lat < self._lat + 1:\r
+ bb_lat = bb_lat + self._vstep\r
+ mlat = mlat+1\r
+ while bb_lng < self._long + 1:\r
+ bb_lng = bb_lng + mstr_zl_18\r
+ mlng = mlng+1\r
+ mstr_msg("orthographic", "Max lat tile: " + str(mlat) + " - max lng tile: " + str(mlng))\r
+ maxlatlng = [ mlat, mlng ]\r
+ \r
+ # The object that handles it all\r
+ xpscn = mstr_xp_scenery(self._lat, self._long, maxlatlng[0], maxlatlng[1], self._vstep, self._latlngfld)\r
+ mstr_msg("orthographic", "[X-Plane] Scenery object instantiated")\r
+\r
+ # Generate the script\r
+ #xpscn.build_mesh_script()\r
+ #mstr_msg("orthographic", "[X-Plane] Mesh script written")\r
+\r
+ # Download LIDAR scan from our endpoint\r
+ #xpscn.acquire_elevation_data()\r
+ #mstr_msg("orthographic", "[X-Plane] LIDAR scan acquired")\r
+\r
+ # Download required XES data\r
+ #xpscn.acquire_xes_data()\r
+ #mstr_msg("orthographic", "[X-Plane] MeshTool XES data acquired")\r
+\r
+ # Generate the .ter files\r
+ xpscn.build_ter_files()\r
+ mstr_msg("orthographic", "[X-Plane] Terrain files (.ter) generated and written")\r
+\r
+ # Build mesh\r
+ #xpscn.build_mesh()\r
+ #xpscn._dsf_test()\r
+ #xpscn.build_and_convert_dsf()\r
+\r
+ # And lastly, generate the mesh\r
+ xpscn.generate_terrain_mesh()\r
+ mstr_msg("orthographic", "[X-Plane] Scenery mesh constructed")\r
+\r
+ # Convert the DSF\r
+ xpscn.build_and_convert_dsf()\r
+ mstr_msg("orthographic", "[X-Plane] DSF generated")\r
+ \r
+ \r
# Checks which layers need to be generated, and what kind of layer it is\r
def determineLayerWork(self, xmlobj):\r
\r
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
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
# Build the correct file name for the elevation model
return fn
- # Generate the mesh script for the ortho photos
- def build_mesh_script(self):
- scr = mstr_datafolder + "z_orthographic/data/meshscript.txt"
- # Before we blast all these lines into the file, we need to make sure they do not exist already
- write_lines = True
-
- if os.path.isfile(scr) == True:
- fnlines = []
- with open(scr) as textfile:
- fnlines = textfile.readlines()
-
- for line in fnlines:
- l = line.split(" ")
- if l[2] == str(self._lng) and l[3] == str(self._lat):
- write_lines = False
- break
- else:
- open(scr, 'a').close()
-
- # If we did not find the initial corner coordinate in the script, we can go ahead
- if write_lines == True:
- mstr_msg("xp_scenery", "[X-Plane] Writing mesh script file")
- # We basically run through all tiles and note down the position of the orthos
- # as needed by X-Plane.
- cur_lat = self._lat
- cur_lng = self._lng
- 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.
- if os.path.isfile(mstr_datafolder + "z_orthographic/" + self._latlngfld + "/orthos/" + str(lat) + "_" + str(lng) + ".dds" ) == True:
- # The '1' after 'ORTHOPHOTO' defines we want water underneath transparent parts of the DDS texture/ortho.
- # This ensures that even if the mesh does not include information for there being a water body,
- # we will get 100% correct representation of the water bodies.
- scrtxt = "ORTHOPHOTO 1 " + str(cur_lng) + " " + str(cur_lat) + " " + str(round(cur_lng+mstr_zl_18, 6)) + " " + str(cur_lat) + " " + str(round(cur_lng+mstr_zl_18, 6)) + " " + str(round(cur_lat+self._vstep, 6)) + " " + str(cur_lng) + " " + str(round(cur_lat+self._vstep, 6)) + " terrain/" + self._latlngfld + "/" + str(lat) + "_" + str(lng) + ".ter\n"
- with open(scr, 'a') as textfile:
- textfile.write(scrtxt)
+ # 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"
+ scr = mstr_datafolder + "z_orthographic/data/" + self._latlngfld + "/meshscript.txt"
+ cmd = mstr_xp_dsftool + " --text2dsf " + meshtxt + " '" + mstr_datafolder + "z_orthographic/Earth nav data/" + llf + "/" + self._latlngfld + ".dsf'"
+ os.system(cmd)
- 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] Mesh script completed")
+ # 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
mstr_msg("xp_scenery", "[X-Plane] XES data acquired")
- # This builds the entire mesh in one go
- def build_mesh(self):
- mstr_msg("xp_scenery", "[X-Plane] Building DSF mesh")
- end_bt = self.find_earthnavdata_number()
- btlfn = str(self.xplane_latlng_folder(end_bt))
- xp_folder = self.xplane_latlng_folder([self._lat, self._lng])
- scr = mstr_datafolder + "z_orthographic/data/meshscript.txt"
- wd = mstr_datafolder + "z_orthographic/data"
- dsf = mstr_datafolder + "z_orthographic/Earth nav data/" + btlfn + "/" + xp_folder
- xesfn = self.build_dem_filename(True)
-
- # The main command to build the mesh
- cmd = mstr_xp_meshtool + " \"" + scr + "\" \"" + mstr_datafolder + "_cache/" + xesfn + "\"" + " \"" + mstr_datafolder + "_cache/" + self._demfn + "\" \"" + wd + "\" \"" + dsf + ".dsf\""
-
- os.system(cmd)
-
- mstr_msg("xp_scenery", "[X-Plane] Mesh construction complete")
-
# This generates all .ter files
def build_ter_files(self):
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\n"
- terstr = terstr + "800\n"
- terstr = terstr + "TERRAIN\n"
- terstr = terstr + "\n"
- terstr = terstr + "BASE_TEX_NOWRAP ../orthos/"+xp_folder+"/"+str(lat)+"_"+str(lng)+".dds\n"
+ 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 + "TEXTURE_NOWRAP ../../orthos/" + self._latlngfld + "/" + str(lat)+"_"+str(lng)+".dds\r\n"
if mstr_xp_scn_normalmaps == True:
- terstr = terstr + "TEXTURE_NORMAL ../normals/"+xp_folder+"/"+str(lat)+"_"+str(lng)+".dds\n"
+ terstr = terstr + "NORMAL_TEX 1.0 ../../normals/" + self._latlngfld + "/" + str(lat)+"_"+str(lng)+".png\r\n"
- terfln = mstr_datafolder + "z_orthographic/terrain/"+xp_folder+"/"+str(lat)+"_"+str(lng)+".ter"
+ 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 6/0\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)
+
+ 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) == True:
+ dsfstr = "TERRAIN_DEF terrain/" + self._latlngfld + "/" + str(lat) + "_" + str(lng) + ".ter\r\n"
+
+ # Let's check if this tile needs water beneath
+ needs_water = False
+ thistile = Image.open(ddsf)
+ tile_pix = thistile.load()
+ for y in range(thistile.height):
+ for x in range(thistile.width):
+ clr = tile_pix[x,y]
+ if clr[3] == 0:
+ needs_water = True
+ break
+
+ if needs_water == True:
+ dsfstr = dsfstr + "TERRAIN_DEF terrain_Water\r\n"
+
+ with open(meshtxt, 'a') as textfile:
+ textfile.write(dsfstr)
+
+
+ # OK. So. Let's build the mesh.
+
+ """
+ # First, the ground water mesh
+ with open(meshtxt, 'a') as textfile:
+ textfile.write("BEGIN_PATCH 0 0.000000 -1.000000 1 5\r\n")
+
+ # Vertical row (Latitude Row)
+ for lat_r in range(0, len(self._demcoord)-2):
+
+ # Horizontal row (Longitude Column)
+ for lng_c in range(0, len(self._demcoord)-2):
+
+ # Lat/lng coordinate
+ lat_crd = self._demcoord[lat_r][lng_c][0]
+ lat_crd_t = self._demcoord[lat_r+1][lng_c][0]
+ lng_crd = self._demcoord[lat_r][lng_c][1]
+ lng_crd_r = self._demcoord[lat_r][lng_c+1][1]
+
+ # Coords of triangle vertices
+ # 0 - Longitude
+ # 1 - Latitude
+ # 2 - Height in m
+ t1_v1 = [ lng_crd_r, lat_crd, 0 ]
+ t1_v2 = [ lng_crd, lat_crd_t, 0 ]
+ t1_v3 = [ lng_crd_r, lat_crd_t, 0 ]
+ t2_v1 = [ lng_crd, lat_crd_t, 0 ]
+ t2_v2 = [ lng_crd_r, lat_crd, 0 ]
+ t2_v3 = [ lng_crd, lat_crd, 0 ]
+
+
+ t1_v1 = [ lng_crd_r, lat_crd, self._demcoord[lat_r][lng_c+1][2] ]
+ t1_v2 = [ lng_crd, lat_crd_t, self._demcoord[lat_r+1][lng_c][2] ]
+ t1_v3 = [ lng_crd_r, lat_crd_t, self._demcoord[lat_r+1][lng_c+1][2] ]
+ t2_v1 = [ lng_crd, lat_crd_t, self._demcoord[lat_r+1][lng_c][2] ]
+ t2_v2 = [ lng_crd_r, lat_crd, self._demcoord[lat_r][lng_c+1][2] ]
+ t2_v3 = [ lng_crd, lat_crd, self._demcoord[lat_r][lng_c][2] ]
+
+
+ # 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\r\n"
+ t_str = t_str + "PATCH_VERTEX " + str(t1_v2[0]) + " " + str(t1_v2[1]) + " " + str(t1_v2[2]) + " 0.000015 0.000015\r\n"
+ t_str = t_str + "PATCH_VERTEX " + str(t1_v3[0]) + " " + str(t1_v3[1]) + " " + str(t1_v3[2]) + " 0.000015 0.000015\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\r\n"
+ t_str = t_str + "PATCH_VERTEX " + str(t2_v2[0]) + " " + str(t2_v2[1]) + " " + str(t2_v2[2]) + " 0.000015 0.000015\r\n"
+ t_str = t_str + "PATCH_VERTEX " + str(t2_v3[0]) + " " + str(t2_v3[1]) + " " + str(t2_v3[2]) + " 0.000015 0.000015\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 = ""
+
+ # Water mesh ends
+ with open(meshtxt, 'a') as textfile:
+ textfile.write("END PATCH\r\n")
+ """
+
+ # 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 + "/1_1.dds"
+ ddsf = mstr_datafolder + "z_orthographic/orthos/" + self._latlngfld + "/" + str(lat) + "_" + str(lng) + ".dds"
+ if os.path.isfile(ddsf) == True:
+
+ scangrid = self.find_height_scan_start_end_points([ self._lat+((lat-1)*self._vstep), self._lng+((lng-1)*mstr_zl_18) ])
+ #sloped = self.build_sloped_scangrid(scangrid)
+
+ # 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
+
+ # Height values
+ hgt_bl = 0
+ hgt_br = 0
+ hgt_tr = 0
+ hgt_tl = 0
+
+ # 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 = ""
+
+
+ # Height value:
+ # hgtindex = self.find_height_for_coord([lat, lng])
+ # height = self._demcoord[hgtindex[0]][hgtindex[1]][2]
+
+
+ # End this patch
+ with open(meshtxt, 'a') as textfile:
+ textfile.write("END PATCH\r\n")
+
+ # Increase patch number
+ curpatch = curpatch + 1
+
+
+ # Let's check if this tile needs water beneath
+ needs_water = False
+ thistile = Image.open(ddsf)
+ tile_pix = thistile.load()
+ for y in range(thistile.height):
+ for x in range(thistile.width):
+ clr = tile_pix[x,y]
+ if clr[3] == 0:
+ needs_water = True
+ break
+
+
+ if needs_water == True:
+
+ # Begin a new patch
+ with open(meshtxt, 'a') as textfile:
+ textfile.write("BEGIN_PATCH " + str(curpatch) + " 0.000000 -1.000000 1 5\r\n")
+
+ # 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
+
+ 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] - .1, 6)
+ hgt_br = round(self._demcoord[ hgt_br_idx[0] ][ hgt_br_idx[1] ][2] - .1, 6)
+ hgt_tr = round(self._demcoord[ hgt_tr_idx[0] ][ hgt_tr_idx[1] ][2] - .1, 6)
+ hgt_tl = round(self._demcoord[ hgt_tl_idx[0] ][ hgt_tl_idx[1] ][2] - .1, 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\r\n"
+ t_str = t_str + "PATCH_VERTEX " + str(t1_v2[0]) + " " + str(t1_v2[1]) + " " + str(t1_v2[2]) + " 0.000015 0.000015\r\n"
+ t_str = t_str + "PATCH_VERTEX " + str(t1_v3[0]) + " " + str(t1_v3[1]) + " " + str(t1_v3[2]) + " 0.000015 0.000015\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\r\n"
+ t_str = t_str + "PATCH_VERTEX " + str(t2_v2[0]) + " " + str(t2_v2[1]) + " " + str(t2_v2[2]) + " 0.000015 0.000015\r\n"
+ t_str = t_str + "PATCH_VERTEX " + str(t2_v3[0]) + " " + str(t2_v3[1]) + " " + str(t2_v3[2]) + " 0.000015 0.000015\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)
+
+ # 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
+
+
+ """
+ t = self._lat
+ while t < startcoord[0]:
+ t = t + self._vstep
+ startend[0] = startend[0]+1
+
+ t = self._lat
+ while t < startcoord[0]+self._vstep:
+ t = t + self._vstep
+ startend[1] = startend[1]+1
+
+ t = self._lng
+ while t < startcoord[1]:
+ t = t + mstr_zl_18
+ startend[2] = startend[2]+1
+
+ t = self._lng
+ while t < startcoord[1]+mstr_zl_18:
+ t = t + mstr_zl_18
+ startend[3] = startend[3]+1
+
+ # Some corrections
+ startend[0] = startend[0]-1
+ if startend[0] < 0: startend[0] = 0
+
+ startend[1] = startend[1]+1
+ if startend[1] > len(self._demdata)-1: startend[1] = startend[1] = len(self._demdata)-1
+
+ startend[2] = startend[2]-1
+ if startend[2] < 0: startend[2] = 0
+
+ startend[3] = startend[3]+1
+ if startend[3] > len(self._demdata)-1: startend[3] = startend[3] = len(self._demdata)-1
+ """
+
+ return startend
+
+
+ # Function to subdivide between two vectors
+ def subdivide_vectors(self, v1, v2, subdivisions):
+ #return np.linspace(v1, v2, subdivisions + 2, axis=0) # +2 to include endpoints
+ return numpy.linspace(v1, v2, subdivisions + 2, axis=0) # +2 to include endpoints
+
+
+ # This build a scangrid with increased resolution, extrapolated from existing points.
+ # With this we can accurately depict the height of a point a long the slope of the ground mesh.
+ def build_sloped_scangrid(self, grid):
+
+ # Contains the data as defined by the grid passed in
+ tmp_dem = []
+
+ # Acquire original grid data
+ for l in range(grid[2], grid[3]+1):
+ row = []
+ for c in range(grid[0], grid[1]+1):
+ row.append(self._demcoord[l][c])
+ tmp_dem.append(row)
+
+ # Subdivide the array
+ subdivisions = 10
+ result = []
+ for i in range(len(tmp_dem) - 1):
+ for j in range(len(tmp_dem[i]) - 1):
+ subdivided = self.subdivide_vectors(tmp_dem[i][j], tmp_dem[i][j + 1], subdivisions)
+ if i > 0: # Avoid duplicating start point
+ subdivided = subdivided[1:]
+ result.append(subdivided)
+
+ # Combine all subdivisions into one array
+ result = numpy.vstack(result)
+
+ return result
\ No newline at end of file
marstr / orthographic / commitdiff