From bd6ca0ef78cacb87736a7e567b8da0ecd2524315 Mon Sep 17 00:00:00 2001 From: Javier Cladellas Date: Tue, 16 Jul 2024 09:03:16 +0200 Subject: [PATCH] add mercator projection, optimize and simple examples #30 --- src/notebooks/reader.ipynb | 577 ++++++++++++++++++++++++++++--------- 1 file changed, 447 insertions(+), 130 deletions(-) diff --git a/src/notebooks/reader.ipynb b/src/notebooks/reader.ipynb index ed94290..e4bdda2 100644 --- a/src/notebooks/reader.ipynb +++ b/src/notebooks/reader.ipynb @@ -22,7 +22,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "61ee8919-195b-4d1f-8ff1-a94603596143", "metadata": {}, "outputs": [], @@ -31,7 +31,7 @@ "\n", "# Setup the reader\n", "# We suppose that the necessary files are in the same directory as the notebook\n", - "case_file = \"strasbourg_sm_lod1/strasbourg_sm_lod1/City_Energy_Modeling.case\"\n", + "case_file = \"strasbourg_sm_lod1/City_Energy_Modeling.case\"\n", "# Setup the reader\n", "reader = vtk.vtkEnSightGoldBinaryReader()\n", "reader.SetCaseFileName(case_file)\n", @@ -47,42 +47,174 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "96eede0d-c63e-4187-94e0-f473df853c42", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Number of timesteps: 731\n", - "Number of cells: 212692\n", - "Number of scalar points?: 212692\n" - ] - } - ], + "outputs": [], "source": [ "print(\"Number of timesteps: \", timesteps)\n", "print(\"Number of cells: \", reader.GetOutput().GetBlock(0).GetNumberOfCells())\n", "print(\"Number of scalar points?: \", reader.GetOutput().GetBlock(0).GetCellData().GetArray(\"solar_shading_coeff\").GetNumberOfTuples())" ] }, + { + "cell_type": "markdown", + "id": "15d2113f", + "metadata": {}, + "source": [ + "Get the mesh" + ] + }, { "cell_type": "code", - "execution_count": 3, - "id": "a6e05806-9a22-488a-902d-495e19899620", + "execution_count": null, + "id": "7f7bc88b", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "output = reader.GetOutput().GetBlock(0)\n", + "mesh = np.empty((output.GetNumberOfCells(),3,3))\n", + "for i in range(output.GetNumberOfCells()):\n", + " #TODO: Slow, there should be a better way to convert to numpy\n", + " cell = output.GetCell(i)\n", + " cell_points = cell.GetPoints()\n", + " mesh[i] = np.array([cell_points.GetPoint(i) for i in range(cell_points.GetNumberOfPoints())])" + ] + }, + { + "cell_type": "markdown", + "id": "fb5f26e5", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "212692\n" - ] - } - ], "source": [ - "print(output.GetNumberOfCells())" + "# Inverse Mercator projection" + ] + }, + { + "cell_type": "markdown", + "id": "48772528", + "metadata": {}, + "source": [ + "### Getting the reference coordinates" + ] + }, + { + "cell_type": "markdown", + "id": "26ae5386", + "metadata": {}, + "source": [ + "Read the GIS file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1205003a", + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "\n", + "gis_path = \"gis.json\"\n", + "\n", + "with open(gis_path,\"r\") as f:\n", + " gis = json.load(f)\n", + "\n", + "reference_coordinates = gis[\"referenceCoordinates\"]\n", + "print(\"Reference coordinates: \", reference_coordinates)\n", + "\n", + "del gis" + ] + }, + { + "cell_type": "markdown", + "id": "1ffe1d19", + "metadata": {}, + "source": [ + "Inverse Mercator projection" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fbf0942e", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "\n", + "R_MAJOR = 6378137.0\n", + "R_MINOR=6356752.3142\n", + "ECCENT=np.sqrt(1-(R_MINOR/R_MAJOR)**2)\n", + "COM = 0.5 * ECCENT\n", + "SCALE_FACTOR = 1./np.cos(np.deg2rad(reference_coordinates[0]))\n", + "\n", + "\n", + "def xToLon(x):\n", + " return np.rad2deg(x*SCALE_FACTOR/R_MAJOR)\n", + "\n", + "def yToLat(y):\n", + " return np.rad2deg(np.pi/2 - 2 * np.arctan(np.exp( -y / R_MAJOR)))" + ] + }, + { + "cell_type": "markdown", + "id": "8f7b9065", + "metadata": {}, + "source": [ + "Validate" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8e40681", + "metadata": {}, + "outputs": [], + "source": [ + "latlon_mesh = np.empty((mesh.shape[0],mesh.shape[1],3))\n", + "\n", + "latlon_mesh[:,:,0] = xToLon(mesh[:,:,0]) + reference_coordinates[1]\n", + "latlon_mesh[:,:,1] = yToLat(mesh[:,:,1]) + reference_coordinates[0]\n", + "latlon_mesh[:,:,2] = mesh[:,:,2]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd5dcc29", + "metadata": {}, + "outputs": [], + "source": [ + "import plotly.graph_objects as go\n", + "\n", + "go.Figure(\n", + " data=[\n", + " go.Scattermapbox(\n", + " lon=latlon_mesh[:10000,0,0],\n", + " lat=latlon_mesh[:10000,0,1],\n", + " mode='markers',\n", + " marker=dict(size=5,color='red'),\n", + " )\n", + " ],\n", + " layout = go.Layout(\n", + " margin ={'l':0,'t':0,'b':0,'r':0},\n", + " mapbox = {\n", + " 'style': \"open-street-map\",\n", + " 'center': {'lon': reference_coordinates[1], 'lat': reference_coordinates[0]},\n", + " 'zoom': 13}\n", + " )\n", + ").show()" + ] + }, + { + "cell_type": "markdown", + "id": "69e8e41c", + "metadata": {}, + "source": [ + "Free the memory (Assuming we will convert the coordinates after)" ] }, { @@ -95,35 +227,110 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": null, "id": "db60d12d-b71a-41e7-a11f-14bfafe3d00a", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " timestep x_1 y_1 z_1 x_2 y_2 z_2 \\\n", - "0 0 -539.268982 600.049988 0.0 -538.495972 599.625000 0.0 \n", - "1 0 -541.656982 611.940979 0.0 -527.434998 605.570984 0.0 \n", - "2 0 -527.434998 605.570984 0.0 -538.495972 599.625000 0.0 \n", - "3 0 -487.437988 635.723999 0.0 -490.282013 625.955994 0.0 \n", - "4 0 -488.596985 642.520020 0.0 -487.437988 635.723999 0.0 \n", - "\n", - " x_3 y_3 z_3 scalar area \n", - "0 -538.987976 599.625000 0.0 0.0 0.104548 \n", - "1 -539.268982 600.049988 0.0 0.0 76.950973 \n", - "2 -539.268982 600.049988 0.0 0.0 4.648543 \n", - "3 -518.830994 625.955994 0.0 0.0 139.433298 \n", - "4 -518.830994 625.955994 0.0 0.0 112.334297 \n" - ] - } - ], + "outputs": [], + "source": [ + "# import numpy as np\n", + "# import pandas as pd\n", + "\n", + "# # Initialize an empty list to accumulate the arrays\n", + "# data = []\n", + "\n", + "# for t in range(timesteps):\n", + "# reader.SetTimeValue(time.GetValue(t))\n", + "# reader.Update()\n", + "# output = reader.GetOutput().GetBlock(0)\n", + "# cell_data = output.GetCellData()\n", + "# sc = cell_data.GetArray(\"solar_shading_coeff\")\n", + "\n", + "\n", + "# for i in range(output.GetNumberOfCells()):\n", + "# cell = output.GetCell(i)\n", + "# pts = cell.GetPoints()\n", + "# scalar = sc.GetTuple(i)\n", + "# np_pts = np.array([pts.GetPoint(j) for j in range(pts.GetNumberOfPoints())])\n", + " \n", + "# flattened_pts = np_pts.flatten()\n", + "# combined_data = np.append(flattened_pts, scalar[0])\n", + " \n", + "# # Append the timestep and combined data to the list\n", + "# data.append([t] + combined_data.tolist())\n", + "# # data.append(np.insert(combined_data, 0, t))\n", + "\n", + "# # Create column names for the DataFrame\n", + "# columns = ['timestep'] + ['x_1'] + ['y_1'] + ['z_1'] + ['x_2'] + ['y_2'] + ['z_2'] + ['x_3'] + ['y_3'] + ['z_3'] + ['scalar']\n", + "\n", + "# # Convert the list to a DataFrame\n", + "# df = pd.DataFrame(data, columns=columns)\n", + "\n", + "# # Function to calculate the area of a triangle given its vertices\n", + "# def calculate_triangle_area(p1, p2, p3):\n", + "# # Calculate the vectors for two sides of the triangle\n", + "# v1 = p2 - p1\n", + "# v2 = p3 - p1\n", + "# # Calculate the cross product of the vectors\n", + "# cross_product = np.cross(v1, v2)\n", + "# # Calculate the area of the triangle (0.5 * magnitude of the cross product)\n", + "# area = 0.5 * np.linalg.norm(cross_product)\n", + "# return area\n", + "\n", + "# # Calculate the area for each triangle and add it as a new column\n", + "# areas = []\n", + "# for index, row in df.iterrows():\n", + "# p1 = np.array([row['x_1'], row['y_1'], row['z_1']])\n", + "# p2 = np.array([row['x_2'], row['y_2'], row['z_2']])\n", + "# p3 = np.array([row['x_3'], row['y_3'], row['z_3']])\n", + "# area = calculate_triangle_area(p1, p2, p3)\n", + "# areas.append(area)\n", + "\n", + "# df['area'] = areas\n", + "\n", + "# # Print the resulting DataFrame\n", + "# print(df.head())\n" + ] + }, + { + "cell_type": "markdown", + "id": "d313c300", + "metadata": {}, + "source": [ + "**To DataFrame**\n", + "\n", + "Do this only if you REALLY need to have the data in tabular format. Otherwise, keep working with numpy.\n", + "\n", + "When having a tabular format, you are using approx. 11GB of RAM, when you could easily only use approx 1.25GB." + ] + }, + { + "cell_type": "markdown", + "id": "268438f3", + "metadata": {}, + "source": [ + "Flatten the mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a0dc9d3", + "metadata": {}, + "outputs": [], + "source": [ + "flat_mesh = mesh.reshape(-1,mesh.shape[1]*mesh.shape[2])\n", + "print(\"Flattened mesh shape: \",flat_mesh.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13d071e9", + "metadata": {}, + "outputs": [], "source": [ - "import numpy as np\n", "import pandas as pd\n", "\n", - "# Initialize an empty list to accumulate the arrays\n", "data = []\n", "\n", "for t in range(timesteps):\n", @@ -131,93 +338,203 @@ " reader.Update()\n", " output = reader.GetOutput().GetBlock(0)\n", " cell_data = output.GetCellData()\n", - " sc = cell_data.GetArray(\"solar_shading_coeff\")\n", - "\n", + " sc_field = cell_data.GetArray(\"solar_shading_coeff\")\n", "\n", + " sc_array = np.empty((output.GetNumberOfCells(),1))\n", " for i in range(output.GetNumberOfCells()):\n", - " cell = output.GetCell(i)\n", - " pts = cell.GetPoints()\n", - " scalar = sc.GetTuple(i)\n", - " np_pts = np.array([pts.GetPoint(j) for j in range(pts.GetNumberOfPoints())])\n", - " \n", - " flattened_pts = np_pts.flatten()\n", - " combined_data = np.append(flattened_pts, scalar[0])\n", - " \n", - " # Append the timestep and combined data to the list\n", - " data.append([t] + combined_data.tolist())\n", - " # data.append(np.insert(combined_data, 0, t))\n", - "\n", - "# Create column names for the DataFrame\n", - "columns = ['timestep'] + ['x_1'] + ['y_1'] + ['z_1'] + ['x_2'] + ['y_2'] + ['z_2'] + ['x_3'] + ['y_3'] + ['z_3'] + ['scalar']\n", - "\n", - "# Convert the list to a DataFrame\n", - "df = pd.DataFrame(data, columns=columns)\n", - "\n", - "# Function to calculate the area of a triangle given its vertices\n", - "def calculate_triangle_area(p1, p2, p3):\n", - " # Calculate the vectors for two sides of the triangle\n", - " v1 = p2 - p1\n", - " v2 = p3 - p1\n", - " # Calculate the cross product of the vectors\n", - " cross_product = np.cross(v1, v2)\n", - " # Calculate the area of the triangle (0.5 * magnitude of the cross product)\n", - " area = 0.5 * np.linalg.norm(cross_product)\n", - " return area\n", - "\n", - "# Calculate the area for each triangle and add it as a new column\n", - "areas = []\n", - "for index, row in df.iterrows():\n", - " p1 = np.array([row['x_1'], row['y_1'], row['z_1']])\n", - " p2 = np.array([row['x_2'], row['y_2'], row['z_2']])\n", - " p3 = np.array([row['x_3'], row['y_3'], row['z_3']])\n", - " area = calculate_triangle_area(p1, p2, p3)\n", - " areas.append(area)\n", - "\n", - "df['area'] = areas\n", - "\n", - "# Print the resulting DataFrame\n", - "print(df.head())" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "id": "2b590da1-99db-454c-acd3-b785e2c283b4", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " timestep point_1 point_2 point_3 point_4 point_5 \\\n", - "5104603 23 725.992004 -1084.459961 0.0 725.992004 -1084.459961 \n", - "5104604 23 739.231018 -1185.939941 0.0 737.825989 -1185.089966 \n", - "5104605 23 737.825989 -1185.089966 0.0 737.825989 -1185.089966 \n", - "5104606 23 763.109985 -1150.270020 0.0 760.265015 -1155.369995 \n", - "5104607 23 760.265015 -1155.369995 0.0 760.265015 -1155.369995 \n", - "\n", - " point_6 point_7 point_8 point_9 scalar area \n", - "5104603 5.2 726.484009 -1079.369995 5.2 1.0 13.295592 \n", - "5104604 0.0 739.231018 -1185.939941 6.5 1.0 5.336898 \n", - "5104605 6.5 739.231018 -1185.939941 6.5 1.0 5.336898 \n", - "5104606 0.0 763.109985 -1150.270020 3.0 1.0 8.759744 \n", - "5104607 3.0 763.109985 -1150.270020 3.0 1.0 8.759744 \n" - ] - } - ], - "source": [ - "print(df.tail())" - ] - }, - { - "cell_type": "code", - "execution_count": 13, + " sc_array[i] = sc_field.GetTuple(i)[0]\n", + "\n", + " assert sc_array.shape[0] == flat_mesh.shape[0], \"Solar shading coefficient and mesh do not have the same number of cells\"\n", + "\n", + " tmp_data = np.hstack(((np.ones(sc_array.shape[0])*t).reshape(-1,1),flat_mesh))\n", + " tmp_data = np.hstack((tmp_data,sc_array))\n", + " data.append(tmp_data)\n", + "\n", + "\n", + " if t % 50 == 0:\n", + " print(\"Processed timestep \", t)\n", + "\n", + "\n", + "del tmp_data\n", + "del sc_array" + ] + }, + { + "cell_type": "markdown", + "id": "7e04a14b", + "metadata": {}, + "source": [ + "Creating the DataFrame" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5b74c02", + "metadata": {}, + "outputs": [], + "source": [ + "columns=['timestep','x_1','y_1','z_1','x_2','y_2','z_2','x_3','y_3','z_3','scalar']\n", + "df = pd.DataFrame(np.vstack(data).reshape(-1,len(columns)),columns=columns)\n", + "del data" + ] + }, + { + "cell_type": "markdown", + "id": "0a4bdfd8", + "metadata": {}, + "source": [ + "**Compute the triangle area**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "418092da", + "metadata": {}, + "outputs": [], + "source": [ + "areas = 0.5*np.linalg.norm(\n", + " np.cross(\n", + " mesh[:,1]-mesh[:,0],\n", + " mesh[:,2]-mesh[:,0]\n", + " ), axis=1\n", + ")\n", + "\n", + "df[\"area\"] = np.array([areas]*timesteps).reshape(-1,1)" + ] + }, + { + "cell_type": "markdown", + "id": "2102474b", + "metadata": {}, + "source": [ + "Validate" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24f8b5dd", + "metadata": {}, + "outputs": [], + "source": [ + "df.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d5d849d", + "metadata": {}, + "outputs": [], + "source": [ + "df.tail(10)" + ] + }, + { + "cell_type": "markdown", + "id": "bd856f32", + "metadata": {}, + "source": [ + "To CSV" + ] + }, + { + "cell_type": "code", + "execution_count": null, "id": "7e470e53-0506-43a2-9d12-50deef85fbfa", "metadata": {}, "outputs": [], "source": [ "df.to_csv('initial_dataframe.csv', index=False)" ] + }, + { + "cell_type": "markdown", + "id": "d6d2956d", + "metadata": {}, + "source": [ + "## Using Numpy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9d2f870", + "metadata": {}, + "outputs": [], + "source": [ + "solar_shading_in_time = np.empty((mesh.shape[0],timesteps))\n", + "\n", + "for t in range(timesteps):\n", + " reader.SetTimeValue(time.GetValue(t))\n", + " reader.Update()\n", + " output = reader.GetOutput().GetBlock(0)\n", + " cell_data = output.GetCellData()\n", + " sc_field = cell_data.GetArray(\"solar_shading_coeff\")\n", + "\n", + " for i in range(output.GetNumberOfCells()):\n", + " solar_shading_in_time[i,t] = sc_field.GetTuple(i)[0]" + ] + }, + { + "cell_type": "markdown", + "id": "2a03de87", + "metadata": {}, + "source": [ + "## Some example statistics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b44a3be", + "metadata": {}, + "outputs": [], + "source": [ + "sample_size = 10000\n", + "latlon_mesh_sample = latlon_mesh[:sample_size]\n", + "solar_shading_sample = solar_shading_in_time[:sample_size]\n", + "t = 130\n", + "\n", + "#This takes all faces, so points in a 2d plot will be repeated." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cef2571b", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "go.Figure(\n", + " data=[\n", + " go.Scattermapbox(\n", + " lon=latlon_mesh_sample[:,0,0],\n", + " lat=latlon_mesh_sample[:,0,1],\n", + " mode='markers',\n", + " marker=dict(size=6,color=solar_shading_sample[:,t],colorscale='rdbu',cmin=0,cmax=1),\n", + " )\n", + " ],\n", + " layout = go.Layout(\n", + " margin ={'l':0,'t':0,'b':0,'r':0},\n", + " mapbox = {\n", + " 'style': \"open-street-map\",\n", + " 'center': {'lon': reference_coordinates[1], 'lat': reference_coordinates[0]},\n", + " 'zoom': 13},\n", + " )\n", + ").show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32c5c2e4", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -236,7 +553,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.3" + "version": "3.9.6" } }, "nbformat": 4,