# /********************************************************************** # * * # * Copyright (c) 2021 Javier Braña * # * * # * This program is free software; you can redistribute it and/or modify* # * it under the terms of the GNU Lesser General Public License (LGPL) * # * as published by the Free Software Foundation; either version 2 of * # * the License, or (at your option) any later version. * # * for detail see the LICENCE text file. * # * * # * This program is distributed in the hope that it will be useful, * # * but WITHOUT ANY WARRANTY; without even the implied warranty of * # * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * # * GNU Library General Public License for more details. * # * * # * You should have received a copy of the GNU Library General Public * # * License along with this program; if not, write to the Free Software * # * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307* # * USA * # * * # *********************************************************************** import json import urllib.request import Draft import FreeCAD import FreeCADGui from PySide import QtCore, QtGui from PySide.QtCore import QT_TRANSLATE_NOOP try: _fromUtf8 = QtCore.QString.fromUtf8 except AttributeError: def _fromUtf8(s): return s import os from PVPlantResources import DirIcons as DirIcons import PVPlantSite def get_elevation_from_oe(coordinates): # v1 deepseek """Obtiene elevaciones de Open-Elevation API y devuelve vectores FreeCAD en coordenadas UTM. Args: coordinates (list): Lista de tuplas con coordenadas (latitud, longitud) Returns: list: Lista de vectores FreeCAD con coordenadas UTM y elevación (en milímetros) o lista vacía en caso de error. """ if not coordinates: return [] import requests import utm from requests.exceptions import RequestException # Construcción más eficiente de parámetros locations = "|".join([f"{lat:.6f},{lon:.6f}" for lat, lon in coordinates]) try: response = requests.get( url="https://api.open-elevation.com/api/v1/lookup", params={'locations': locations}, timeout=20, verify=True ) response.raise_for_status() # Lanza excepción para códigos 4xx/5xx except RequestException as e: print(f"Error en la solicitud: {str(e)}") return [] try: data = response.json() except ValueError: print("Respuesta JSON inválida") return [] if "results" not in data or len(data["results"]) != len(coordinates): print("Formato de respuesta inesperado") return [] points = [] for result in data["results"]: try: # Conversión UTM con manejo de errores easting, northing, _, _ = utm.from_latlon( result["latitude"], result["longitude"] ) points.append(FreeCAD.Vector(round(easting), # Convertir metros a milímetros round(northing), round(result["elevation"])) * 1000) except Exception as e: print(f"Error procesando coordenadas: {str(e)}") continue return points def getElevationFromOE(coordinates): """Obtiene elevaciones de Open-Elevation API y devuelve vectores FreeCAD en coordenadas UTM.""" import certifi from requests.exceptions import RequestException if len(coordinates) == 0: return None from requests import get import utm locations_str="" total = len(coordinates) - 1 for i, point in enumerate(coordinates): locations_str += '{:.6f},{:.6f}'.format(point[0], point[1]) if i != total: locations_str += '|' query = 'https://api.open-elevation.com/api/v1/lookup?locations=' + locations_str try: r = get(query, timeout=20, verify=certifi.where()) # <-- Corrección aquí except RequestException as e: points = [] for i, point in enumerate(coordinates): c = utm.from_latlon(point[0], point[1]) points.append(FreeCAD.Vector(round(c[0], 0), round(c[1], 0), 0) * 1000) return points # Only get the json response in case of 200 or 201 points = [] if r.status_code == 200 or r.status_code == 201: results = r.json() for point in results["results"]: c = utm.from_latlon(point["latitude"], point["longitude"]) v = FreeCAD.Vector(round(c[0], 0), round(c[1], 0), round(point["elevation"], 0)) * 1000 points.append(v) return points def getSinglePointElevationFromBing(lat, lng): #http://dev.virtualearth.net/REST/v1/Elevation/List?points={lat1,long1,lat2,long2,latN,longnN}&heights={heights}&key={BingMapsAPIKey} import utm source = "http://dev.virtualearth.net/REST/v1/Elevation/List?points=" source += str(lat) + "," + str(lng) source += "&heights=sealevel" source += "&key=AmsPZA-zRt2iuIdQgvXZIxme2gWcgLaz7igOUy7VPB8OKjjEd373eCnj1KFv2CqX" import requests response = requests.get(source) ans = response.text s = json.loads(ans) print(s) res = s['resourceSets'][0]['resources'][0]['elevations'] for elevation in res: c = utm.from_latlon(lat, lng) v = FreeCAD.Vector( round(c[0] * 1000, 0), round(c[1] * 1000, 0), round(elevation * 1000, 0)) return v def getGridElevationFromBing(polygon, lat, lng, resolution = 1000): #http://dev.virtualearth.net/REST/v1/Elevation/Polyline?points=35.89431,-110.72522,35.89393,-110.72578,35.89374,-110.72606,35.89337,-110.72662 # &heights=ellipsoid&samples=10&key={BingMapsAPIKey} import utm import math import requests geo = utm.from_latlon(lat, lng) # result = (679434.3578335291, 4294023.585627955, 30, 'S') # EASTING, NORTHING, ZONE NUMBER, ZONE LETTER #StepsXX = int((polygon.Shape.BoundBox.XMax - polygon.Shape.BoundBox.XMin) / (resolution*1000)) points = [] yy = polygon.Shape.BoundBox.YMax while yy > polygon.Shape.BoundBox.YMin: xx = polygon.Shape.BoundBox.XMin while xx < polygon.Shape.BoundBox.XMax: StepsXX = int(math.ceil((polygon.Shape.BoundBox.XMax - xx) / resolution)) if StepsXX > 1000: StepsXX = 1000 xx1 = xx + 1000 * resolution else: xx1 = xx + StepsXX * resolution point1 = utm.to_latlon(xx / 1000, yy / 1000, geo[2], geo[3]) point2 = utm.to_latlon(xx1 / 1000, yy / 1000, geo[2], geo[3]) source = "http://dev.virtualearth.net/REST/v1/Elevation/Polyline?points=" source += "{lat1},{lng1}".format(lat1=point1[0], lng1=point1[1]) source += "," source += "{lat2},{lng2}".format(lat2=point2[0], lng2=point2[1]) source += "&heights=sealevel" source += "&samples={steps}".format(steps=StepsXX) source += "&key=AmsPZA-zRt2iuIdQgvXZIxme2gWcgLaz7igOUy7VPB8OKjjEd373eCnj1KFv2CqX" response = requests.get(source) ans = response.text # +# to do: error handling - wait and try again s = json.loads(ans) res = s['resourceSets'][0]['resources'][0]['elevations'] i = 0 for elevation in res: v = FreeCAD.Vector(xx + resolution * i, yy, round(elevation * 1000, 4)) points.append(v) i += 1 xx = xx1 + resolution # para no repetir un mismo punto yy -= resolution return points def getSinglePointElevation(lat, lon): source = "https://maps.googleapis.com/maps/api/elevation/json?locations=" source += str(lat) + "," + str(lon) source += "&key=AIzaSyB07X6lowYJ-iqyPmaFJvr-6zp1J63db8U" #print (source) #response = request.urlopen(source) #ans = response.read() import requests response = requests.get(source) ans = response.text # +# to do: error handling - wait and try again s = json.loads(ans) res = s['results'] from geopy.distance import geodesic for r in res: reference = (0.0, 0.0) v = FreeCAD.Vector( round(geodesic(reference, (0.0, r['location']['lng'])).m, 2), round(geodesic(reference, (r['location']['lat'], 0.0)).m, 2), round(r['elevation'] * 1000, 2) ) return v def _getSinglePointElevation(lat, lon): source = "https://maps.googleapis.com/maps/api/elevation/json?locations=" source += str(lat) + "," + str(lon) source += "&key=AIzaSyB07X6lowYJ-iqyPmaFJvr-6zp1J63db8U" #print (source) #response = request.urlopen(source) #ans = response.read() import requests response = requests.get(source) ans = response.text # +# to do: error handling - wait and try again s = json.loads(ans) res = s['results'] import pymap3d as pm for r in res: x, y, z = pm.geodetic2ecef(round(r['location']['lng'], 2), round(r['location']['lat'], 2), 0) v = FreeCAD.Vector(x,y,z) return v def getSinglePointElevation1(lat, lon): source = "https://maps.googleapis.com/maps/api/elevation/json?locations=" source += str(lat) + "," + str(lon) source += "&key=AIzaSyB07X6lowYJ-iqyPmaFJvr-6zp1J63db8U" #response = urllib.request.urlopen(source) #ans = response.read() import requests response = requests.get(source) ans = response.text # +# to do: error handling - wait and try again s = json.loads(ans) res = s['results'] for r in res: c = tm.fromGeographic(r['location']['lat'], r['location']['lng']) v = FreeCAD.Vector( round(c[0], 4), round(c[1], 4), round(r['elevation'] * 1000, 2) ) return v def getSinglePointElevationUtm(lat, lon): source = "https://maps.googleapis.com/maps/api/elevation/json?locations=" source += str(lat) + "," + str(lon) source += "&key=AIzaSyB07X6lowYJ-iqyPmaFJvr-6zp1J63db8U" print(source) #response = urllib.request.urlopen(source) #ans = response.read() import requests response = requests.get(source) ans = response.text # +# to do: error handling - wait and try again s = json.loads(ans) res = s['results'] print (res) import utm for r in res: c = utm.from_latlon(r['location']['lat'], r['location']['lng']) v = FreeCAD.Vector( round(c[0] * 1000, 4), round(c[1] * 1000, 4), round(r['elevation'] * 1000, 2)) print (v) return v def getElevationUTM(polygon, lat, lng, resolution = 10000): import utm geo = utm.from_latlon(lat, lng) # result = (679434.3578335291, 4294023.585627955, 30, 'S') # EASTING, NORTHING, ZONE NUMBER, ZONE LETTER StepsXX = int((polygon.Shape.BoundBox.XMax - polygon.Shape.BoundBox.XMin) / (resolution*1000)) points = [] yy = polygon.Shape.BoundBox.YMax while yy > polygon.Shape.BoundBox.YMin: # utm.to_latlon(EASTING, NORTHING, ZONE NUMBER, ZONE LETTER). # result = (LATITUDE, LONGITUDE) point1 = utm.to_latlon(polygon.Shape.BoundBox.XMin / 1000, yy / 1000, geo[2], geo[3]) point2 = utm.to_latlon(polygon.Shape.BoundBox.XMax / 1000, yy / 1000, geo[2], geo[3]) source = "https://maps.googleapis.com/maps/api/elevation/json?path=" source += "{a},{b}".format(a = point1[0], b = point1[1]) source += "|" source += "{a},{b}".format(a = point2[0], b = point2[1]) source += "&samples={a}".format(a = StepsXX) source += "&key=AIzaSyB07X6lowYJ-iqyPmaFJvr-6zp1J63db8U" import requests response = requests.get(source) ans = response.text # +# to do: error handling - wait and try again s = json.loads(ans) res = s['results'] for r in res: c = utm.from_latlon(r['location']['lat'], r['location']['lng']) v = FreeCAD.Vector( round(c[0] * 1000, 2), round(c[1] * 1000, 2), round(r['elevation'] * 1000, 2) ) points.append(v) yy -= (resolution*1000) FreeCAD.activeDocument().recompute() return points def getElevation1(polygon,resolution=10): StepsXX = int((polygon.Shape.BoundBox.XMax - polygon.Shape.BoundBox.XMin) / (resolution * 1000)) points = [] yy = polygon.Shape.BoundBox.YMax while yy > polygon.Shape.BoundBox.YMin: point1 = tm.toGeographic(polygon.Shape.BoundBox.XMin, yy) point2 = tm.toGeographic(polygon.Shape.BoundBox.XMax, yy) source = "https://maps.googleapis.com/maps/api/elevation/json?path=" source += "{a},{b}".format(a = point1[0], b = point1[1]) source += "|" source += "{a},{b}".format(a = point2[0], b = point2[1]) source += "&samples={a}".format(a = StepsXX) source += "&key=AIzaSyB07X6lowYJ-iqyPmaFJvr-6zp1J63db8U" try: #response = urllib.request.urlopen(source) #ans = response.read() import requests response = requests.get(source) ans = response.text # +# to do: error handling - wait and try again s = json.loads(ans) res = s['results'] except: continue #points = [] for r in res: c = tm.fromGeographic(r['location']['lat'], r['location']['lng']) v = FreeCAD.Vector( round(c[0], 2), round(c[1], 2), round(r['elevation'] * 1000, 2) ) points.append(v) FreeCAD.activeDocument().recompute() yy -= (resolution*1000) return points ## download the heights from google: def getElevation(lat, lon, b=50.35, le=11.17, size=40): #https://maps.googleapis.com/maps/api/elevation/json?path=36.578581,-118.291994|36.23998,-116.83171&samples=3&key=YOUR_API_KEY #https://maps.googleapis.com/maps/api/elevation/json?locations=39.7391536,-104.9847034&key=YOUR_API_KEY source = "https://maps.googleapis.com/maps/api/elevation/json?path=" source += str(b-size*0.001) + "," + str(le) + "|" + str(b+size*0.001) + "," + str(le) source += "&samples=" + str(100) source += "&key=AIzaSyB07X6lowYJ-iqyPmaFJvr-6zp1J63db8U" response = urllib.request.urlopen(source) ans = response.read() # +# to do: error handling - wait and try again s = json.loads(ans) res = s['results'] from geopy.distance import geodesic points = [] for r in res: reference = (0.0, 0.0) v = FreeCAD.Vector( round(geodesic(reference, (0.0, r['location']['lat'])).m, 2), round(geodesic(reference, (r['location']['lng'], 0.0)).m, 2), round(r['elevation'] * 1000, 2) - baseheight ) points.append(v) line = Draft.makeWire(points, closed=False, face=False, support=None) line.ViewObject.Visibility = False #FreeCAD.activeDocument().recompute() FreeCADGui.updateGui() return FreeCAD.activeDocument().ActiveObject class _ImportPointsTaskPanel: def __init__(self, obj = None): self.obj = None self.Boundary = None self.select = 0 self.filename = "" # form: self.form1 = FreeCADGui.PySideUic.loadUi(os.path.dirname(__file__) + "/PVPlantImportGrid.ui") self.form1.radio1.toggled.connect(lambda: self.mainToggle(self.form1.radio1)) self.form1.radio2.toggled.connect(lambda: self.mainToggle(self.form1.radio2)) self.form1.radio1.setChecked(True) # << --------------Poner al final para que no dispare antes de crear los componentes a los que va a llamar #self.form.buttonAdd.clicked.connect(self.add) self.form1.buttonDEM.clicked.connect(self.openFileDEM) self.form2 = FreeCADGui.PySideUic.loadUi(os.path.dirname(__file__) + "/PVPlantCreateTerrainMesh.ui") #self.form2.buttonAdd.clicked.connect(self.add) self.form2.buttonBoundary.clicked.connect(self.addBoundary) #self.form = [self.form1, self.form2] self.form = self.form1 ''' future: def retranslateUi(self, dialog): self.form1.setWindowTitle("Configuracion del Rack") self.labelModule.setText(QtGui.QApplication.translate("PVPlant", "Modulo:", None)) self.labelModuleLength.setText(QtGui.QApplication.translate("PVPlant", "Longitud:", None)) self.labelModuleWidth.setText(QtGui.QApplication.translate("PVPlant", "Ancho:", None)) self.labelModuleHeight.setText(QtGui.QApplication.translate("PVPlant", "Alto:", None)) self.labelModuleFrame.setText(QtGui.QApplication.translate("PVPlant", "Ancho del marco:", None)) self.labelModuleColor.setText(QtGui.QApplication.translate("PVPlant", "Color del modulo:", None)) self.labelModules.setText(QtGui.QApplication.translate("Arch", "Colocacion de los Modulos", None)) self.labelModuleOrientation.setText(QtGui.QApplication.translate("Arch", "Orientacion del modulo:", None)) self.labelModuleGapX.setText(QtGui.QApplication.translate("Arch", "Separacion Horizontal (mm):", None)) self.labelModuleGapY.setText(QtGui.QApplication.translate("Arch", "Separacion Vertical (mm):", None)) self.labelModuleRows.setText(QtGui.QApplication.translate("Arch", "Filas de modulos:", None)) self.labelModuleCols.setText(QtGui.QApplication.translate("Arch", "Columnas de modulos:", None)) self.labelRack.setText(QtGui.QApplication.translate("Arch", "Configuracion de la estructura", None)) self.labelRackType.setText(QtGui.QApplication.translate("Arch", "Tipo de estructura:", None)) self.labelLevel.setText(QtGui.QApplication.translate("Arch", "Nivel:", None)) self.labelOffset.setText(QtGui.QApplication.translate("Arch", "Offset", None)) ''' def add(self): sel = FreeCADGui.Selection.getSelection() if len(sel) > 0: self.obj = sel[0] self.lineEdit1.setText(self.obj.Label) def addBoundary(self): sel = FreeCADGui.Selection.getSelection() if len(sel) > 0: self.Boundary = sel[0] self.form2.editBoundary.setText(self.Boundary.Label) def openFileDEM(self): filters = "Esri ASC (*.asc);;CSV (*.csv);;All files (*.*)" filename = QtGui.QFileDialog.getOpenFileName(None, "Open DEM,", "", filters) self.filename = filename[0] self.form1.editDEM.setText(filename[0]) def mainToggle(self, radiobox): if radiobox is self.form1.radio1: self.select = 0 self.form1.gbLocalFile.setVisible(True) elif radiobox is self.form1.radio2: self.select = 1 self.form1.gbLocalFile.setVisible(True) def accept(self): from datetime import datetime starttime = datetime.now() site = PVPlantSite.get() try: PointGroups = FreeCAD.ActiveDocument.Point_Groups except: PointGroups = FreeCAD.ActiveDocument.addObject("App::DocumentObjectGroup", 'Point_Groups') PointGroups.Label = "Point Groups" PointGroup = FreeCAD.ActiveDocument.addObject('Points::Feature', "Point_Group") PointGroup.Label = "Land_Grid_Points" FreeCAD.ActiveDocument.Point_Groups.addObject(PointGroup) PointObject = PointGroup.Points.copy() if self.select == 0: # Google or bing or ... #for item in self.obj: #if self.groupbox.isChecked:break resol = FreeCAD.Units.Quantity(self.valueResolution.text()).Value Site = FreeCAD.ActiveDocument.Site pts = getGridElevationFromBing(self.obj, Site.Latitude, Site.Longitude, resol) PointObject.addPoints(pts) PointGroup.Points = PointObject else: if self.filename == "": return import Utils.importDEM as openDEM if self.select == 1: # DEM. import numpy as np root, extension = os.path.splitext(self.filename) if extension.lower() == ".asc": x, y, datavals, cellsize, nodata_value = openDEM.openEsri(self.filename) if self.Boundary: inc_x = self.Boundary.Shape.BoundBox.XLength * 0.05 inc_y = self.Boundary.Shape.BoundBox.YLength * 0.05 min_x = 0 max_x = 0 comp = (self.Boundary.Shape.BoundBox.XMin - inc_x) / 1000 for i in range(nx): if x[i] > comp: min_x = i - 1 break comp = (self.Boundary.Shape.BoundBox.XMax + inc_x) / 1000 for i in range(min_x, nx): if x[i] > comp: max_x = i break min_y = 0 max_y = 0 comp = (self.Boundary.Shape.BoundBox.YMax + inc_y) / 1000 for i in range(ny): if y[i] < comp: max_y = i break comp = (self.Boundary.Shape.BoundBox.YMin - inc_y) / 1000 for i in range(max_y, ny): if y[i] < comp: min_y = i break x = x[min_x:max_x] y = y[max_y:min_y] datavals = datavals[max_y:min_y, min_x:max_x] pts = [] if True: # faster but more memory 46s - 4,25 gb x, y = np.meshgrid(x, y) xx = x.flatten() yy = y.flatten() zz = datavals.flatten() x[:] = 0 y[:] = 0 datavals[:] = 0 pts = [] for i in range(0, len(xx)): pts.append(FreeCAD.Vector(xx[i], yy[i], zz[i]) * 1000) xx[:] = 0 yy[:] = 0 zz[:] = 0 else: # 51s 3,2 gb createmesh = True if createmesh: import Part, Draft lines=[] for j in range(len(y)): edges = [] for i in range(0, len(x) - 1): ed = Part.makeLine(FreeCAD.Vector(x[i], y[j], datavals[j][i]) * 1000, FreeCAD.Vector(x[i + 1], y[j], datavals[j][i + 1]) * 1000) edges.append(ed) #bspline = Draft.makeBSpline(pts) #bspline.ViewObject.hide() line = Part.Wire(edges) lines.append(line) ''' for i in range(0, len(bsplines), 100): p = Part.makeLoft(bsplines[i:i + 100], False, False, False) Part.show(p) ''' p = Part.makeLoft(lines, False, True, False) p = Part.Solid(p) Part.show(p) else: pts = [] for j in range(ny): for i in range(nx): pts.append(FreeCAD.Vector(x[i], y[j], datavals[j][i]) * 1000) elif extension.lower() == ".csv" or extension.lower() == ".txt": # x, y, z from gps pts = openDEM.interpolatePoints(openDEM.openCSV(self.filename)) PointObject.addPoints(pts) PointGroup.Points = PointObject FreeCAD.ActiveDocument.recompute() FreeCADGui.Control.closeDialog() print("tiempo: ", datetime.now() - starttime) def reject(self): FreeCADGui.Control.closeDialog() ## Comandos ----------------------------------------------------------------------------------------------------------- class CommandImportPoints: def GetResources(self): return {'Pixmap': str(os.path.join(DirIcons, "cloud.svg")), 'MenuText': QT_TRANSLATE_NOOP("PVPlant", "Importer Grid"), 'Accel': "B, U", 'ToolTip': QT_TRANSLATE_NOOP("PVPlant", "Creates a cloud of points.")} def IsActive(self): return not FreeCAD.ActiveDocument is None def Activated(self): self.TaskPanel = _ImportPointsTaskPanel() FreeCADGui.Control.showDialog(self.TaskPanel) if FreeCAD.GuiUp: class CommandPointsGroup: def GetCommands(self): return tuple(['ImportPoints' ]) def GetResources(self): return { 'MenuText': QT_TRANSLATE_NOOP("",'Cloud of Points'), 'ToolTip': QT_TRANSLATE_NOOP("",'Cloud of Points') } def IsActive(self): return not FreeCAD.ActiveDocument is None FreeCADGui.addCommand('ImportPoints', CommandImportPoints()) FreeCADGui.addCommand('PointsGroup', CommandPointsGroup())