|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "markdown", |
| 5 | + "id": "f6a6711f", |
| 6 | + "metadata": {}, |
| 7 | + "source": [ |
| 8 | + "# Example using Geography Module\n", |
| 9 | + "\n", |
| 10 | + "FAModel was originally designed to represent a floating array model in arbitrary space but it also has the ability to represent a model in an actual, global space (e.g., a specific floating offshore wind lease area).\n", |
| 11 | + "\n", |
| 12 | + "This file will walk through the options available to a user who wishes to use geography-related design inputs" |
| 13 | + ] |
| 14 | + }, |
| 15 | + { |
| 16 | + "cell_type": "markdown", |
| 17 | + "id": "a145fb42", |
| 18 | + "metadata": {}, |
| 19 | + "source": [ |
| 20 | + "### CRS Setup\n", |
| 21 | + "\n", |
| 22 | + "The first step towards using geography-related tools is to set the proper coordindate-reference system (CRS). There are many CRSs that set up coordinates for the world in either a geographic coordinate system (i.e., latitude and longitude) or a projected coordinate system (i.e., UTM). The main goal is to represent the 3-D coordinates of a globe into a 2-D plane. Each coordinate reference system is going to distort the global map in some way and you're not going to be able to preserve all properties (like distances or areas). That's why Greenland appears so big on some maps. Imagine peeling an orange and then laying the peels out on a table into a rectangle...it's hard to do.\n", |
| 23 | + "\n", |
| 24 | + "There is a Python package called 'pyproj' which we use to store CRS information. As a start, we can set a 'global' CRS based on the conventional lat/long coordinate system, which is signified by a specific EPSG code: 4326" |
| 25 | + ] |
| 26 | + }, |
| 27 | + { |
| 28 | + "cell_type": "code", |
| 29 | + "execution_count": 112, |
| 30 | + "id": "3532056e", |
| 31 | + "metadata": {}, |
| 32 | + "outputs": [], |
| 33 | + "source": [ |
| 34 | + "from pyproj import CRS\n", |
| 35 | + "\n", |
| 36 | + "latlong_crs = CRS.from_epsg(4326)" |
| 37 | + ] |
| 38 | + }, |
| 39 | + { |
| 40 | + "cell_type": "markdown", |
| 41 | + "id": "ab905e5e", |
| 42 | + "metadata": {}, |
| 43 | + "source": [ |
| 44 | + "We can also create other CRSs based on our interests. We can create a UTM CRS (a common 2D projection) if we have a specific location we want to focus on. (Think of UTM CRSs as if you cut the Earth into 60 equidistance slices from the north pole to just above the equator and from the south pole to just below the equator). If you have a specific area of interest, the UTM projection will represent that area very well, but not exactly." |
| 45 | + ] |
| 46 | + }, |
| 47 | + { |
| 48 | + "cell_type": "code", |
| 49 | + "execution_count": 113, |
| 50 | + "id": "ae65b853", |
| 51 | + "metadata": {}, |
| 52 | + "outputs": [], |
| 53 | + "source": [ |
| 54 | + "from pyproj.aoi import AreaOfInterest\n", |
| 55 | + "from pyproj.database import query_utm_crs_info\n", |
| 56 | + "import famodel.geography as geo\n", |
| 57 | + "\n", |
| 58 | + "target_crs = geo.getTargetCRS([-125, -124], [40, 41])" |
| 59 | + ] |
| 60 | + }, |
| 61 | + { |
| 62 | + "cell_type": "markdown", |
| 63 | + "id": "424e638d", |
| 64 | + "metadata": {}, |
| 65 | + "source": [ |
| 66 | + "More appropriately, we can create 'custom' CRSs that create a reference system about a specific longitude and latitude (i.e., the centroid of an offshore wind lease area). This will result in the least amount of distortion in the results." |
| 67 | + ] |
| 68 | + }, |
| 69 | + { |
| 70 | + "cell_type": "code", |
| 71 | + "execution_count": 114, |
| 72 | + "id": "c7395829", |
| 73 | + "metadata": {}, |
| 74 | + "outputs": [], |
| 75 | + "source": [ |
| 76 | + "custom_crs = geo.getCustomCRS(-124.5, 40.5)" |
| 77 | + ] |
| 78 | + }, |
| 79 | + { |
| 80 | + "cell_type": "markdown", |
| 81 | + "id": "08465007", |
| 82 | + "metadata": {}, |
| 83 | + "source": [ |
| 84 | + "### Boundaries\n", |
| 85 | + "\n", |
| 86 | + "Next, we can set up specific boundaries in space (like a lease area, or any area that can be defined by latitude and longitude coordinates). We will use the Humboldt offshore wind lease area as an example.\n", |
| 87 | + "\n", |
| 88 | + "We use the Python package 'geopandas' to store geography-related information, in the conventional pandas format (the main difference between pandas and geopandas is that geopandas uses the 'shapely' python package to store coordinates and shapes).\n", |
| 89 | + "\n", |
| 90 | + "We can create a geopandas dataframe using 'shape files' (.shp) to store the boundary information contained in the shape file," |
| 91 | + ] |
| 92 | + }, |
| 93 | + { |
| 94 | + "cell_type": "code", |
| 95 | + "execution_count": 115, |
| 96 | + "id": "4e2c6338", |
| 97 | + "metadata": {}, |
| 98 | + "outputs": [], |
| 99 | + "source": [ |
| 100 | + "import geopandas as gpd\n", |
| 101 | + "\n", |
| 102 | + "lease_areas = gpd.read_file('../geography/Wind_Lease_Outlines_2_2023.shp')" |
| 103 | + ] |
| 104 | + }, |
| 105 | + { |
| 106 | + "cell_type": "markdown", |
| 107 | + "id": "cff4b8e8", |
| 108 | + "metadata": {}, |
| 109 | + "source": [ |
| 110 | + "and then we can extract certain lease area information from that geopandas dataframe (for futher development, you will have to take a peek inside the geodataframe to determine what exact part of the shapefile you want to extract)." |
| 111 | + ] |
| 112 | + }, |
| 113 | + { |
| 114 | + "cell_type": "code", |
| 115 | + "execution_count": 116, |
| 116 | + "id": "fb432da4", |
| 117 | + "metadata": {}, |
| 118 | + "outputs": [], |
| 119 | + "source": [ |
| 120 | + "lease_area = lease_areas.loc[lease_areas['LEASE_NUMB']=='OCS-P0562 - Provisional']" |
| 121 | + ] |
| 122 | + }, |
| 123 | + { |
| 124 | + "cell_type": "markdown", |
| 125 | + "id": "1676e305", |
| 126 | + "metadata": {}, |
| 127 | + "source": [ |
| 128 | + "You can then extract the latitudes and longitudes of the boundary of that shape," |
| 129 | + ] |
| 130 | + }, |
| 131 | + { |
| 132 | + "cell_type": "code", |
| 133 | + "execution_count": 117, |
| 134 | + "id": "031f5ff7", |
| 135 | + "metadata": {}, |
| 136 | + "outputs": [], |
| 137 | + "source": [ |
| 138 | + "area_longs, area_lats = lease_area.geometry.unary_union.exterior.coords.xy" |
| 139 | + ] |
| 140 | + }, |
| 141 | + { |
| 142 | + "cell_type": "markdown", |
| 143 | + "id": "b7b36192", |
| 144 | + "metadata": {}, |
| 145 | + "source": [ |
| 146 | + "or even get the centroid of the lease area." |
| 147 | + ] |
| 148 | + }, |
| 149 | + { |
| 150 | + "cell_type": "code", |
| 151 | + "execution_count": 118, |
| 152 | + "id": "ef7a7865", |
| 153 | + "metadata": {}, |
| 154 | + "outputs": [ |
| 155 | + { |
| 156 | + "name": "stderr", |
| 157 | + "output_type": "stream", |
| 158 | + "text": [ |
| 159 | + "C:\\Users\\shousner\\AppData\\Local\\Temp\\1\\ipykernel_23620\\3218981639.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n", |
| 160 | + "\n", |
| 161 | + " centroid = ( lease_area.geometry.centroid.values.x[0], lease_area.geometry.centroid.values.y[0] )\n" |
| 162 | + ] |
| 163 | + } |
| 164 | + ], |
| 165 | + "source": [ |
| 166 | + "centroid = ( lease_area.geometry.centroid.values.x[0], lease_area.geometry.centroid.values.y[0] )" |
| 167 | + ] |
| 168 | + }, |
| 169 | + { |
| 170 | + "cell_type": "markdown", |
| 171 | + "id": "21e07e5d", |
| 172 | + "metadata": {}, |
| 173 | + "source": [ |
| 174 | + "All of this is done in the geography module's function ```getLeaseCoords(lease_name)```." |
| 175 | + ] |
| 176 | + }, |
| 177 | + { |
| 178 | + "cell_type": "markdown", |
| 179 | + "id": "e403eb0d", |
| 180 | + "metadata": {}, |
| 181 | + "source": [ |
| 182 | + "Next, we need to convert these latitudes and longitudes to meters, since that is what FAModel is most familiar with. This is where the custom CRS is useful. We can convert the latitudes and longitudes into meters relative to a specific reference point. We run the following converter function to get the boundaries in units of meters.\n", |
| 183 | + "\n", |
| 184 | + "This function uses the geopandas.to_crs() function to convert the data and then ensure that the new coordinates are relative to the input centroid (the input centroid should be the same point used to create the custom CRS).\n", |
| 185 | + "\n", |
| 186 | + "This process can also be reversed using the 'convertMeters2LatLong()' function" |
| 187 | + ] |
| 188 | + }, |
| 189 | + { |
| 190 | + "cell_type": "code", |
| 191 | + "execution_count": 119, |
| 192 | + "id": "616eaecc", |
| 193 | + "metadata": {}, |
| 194 | + "outputs": [ |
| 195 | + { |
| 196 | + "name": "stderr", |
| 197 | + "output_type": "stream", |
| 198 | + "text": [ |
| 199 | + "c:\\code\\floatingarraydesign\\famodel\\famodel\\geography.py:146: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n", |
| 200 | + "\n", |
| 201 | + " centroid = ( lease_area.geometry.centroid.values.x[0], lease_area.geometry.centroid.values.y[0] )\n" |
| 202 | + ] |
| 203 | + } |
| 204 | + ], |
| 205 | + "source": [ |
| 206 | + "lease_name = 'Humboldt_SW'\n", |
| 207 | + "lease_longs, lease_lats, centroid_latlong = geo.getLeaseCoords(lease_name)\n", |
| 208 | + "lease_xs, lease_ys, centroid_m = geo.convertLatLong2Meters(lease_longs, lease_lats, centroid_latlong, latlong_crs, custom_crs, return_centroid=True)" |
| 209 | + ] |
| 210 | + }, |
| 211 | + { |
| 212 | + "cell_type": "markdown", |
| 213 | + "id": "1bcd0ce9", |
| 214 | + "metadata": {}, |
| 215 | + "source": [ |
| 216 | + "### Bathymetry\n", |
| 217 | + "\n", |
| 218 | + "Once the boundaries are set up, the next important part is to set up the bathymetry based on an input GEBCO file. This can all be done by running something like the following " |
| 219 | + ] |
| 220 | + }, |
| 221 | + { |
| 222 | + "cell_type": "code", |
| 223 | + "execution_count": 120, |
| 224 | + "id": "c3b69bac", |
| 225 | + "metadata": {}, |
| 226 | + "outputs": [], |
| 227 | + "source": [ |
| 228 | + "gebco_file = '../geography/gebco_humboldt_2023_n41.3196_s40.3857_w-125.2881_e-123.9642.asc'\n", |
| 229 | + "bath_longs, bath_lats, bath_depths, ncols, nrows = geo.getMapBathymetry(gebco_file)\n", |
| 230 | + "bath_xs, bath_ys, bath_depths = geo.convertBathymetry2Meters(bath_longs, bath_lats, bath_depths, centroid_latlong, centroid_m, latlong_crs, custom_crs, ncols, nrows)" |
| 231 | + ] |
| 232 | + }, |
| 233 | + { |
| 234 | + "cell_type": "markdown", |
| 235 | + "id": "5f03af75", |
| 236 | + "metadata": {}, |
| 237 | + "source": [ |
| 238 | + "This section reads in a GEBCO file and extracts a list of latitudes and longitudes that define a \"rectangle\" of coordinates, where each point within that rectangle has a depth. However, in the 'meters' reference system, this will be a slightly distorted rectangle due to the differences in CRSs. This function also returns the discretization in the number of rows and columns.\n", |
| 239 | + "\n", |
| 240 | + "Those latitude and longitude values are then read into the next function, but only the maximum and minimum values are used to generate a perfect rectangle in the custom_CRS \"inside\" of the distorted geographic rectangle. It can set a new discretization between those maximum and minimum points, creates a python meshgrid, and then converts this new square back into the geographic coordinate system to make a slightly smaller distorted rectangle. It only does this so that it can reference the geographic coordinate system's depths and assign the right depths to the new \"custom CRS\" square (using the getDepthFromBathymetry function). The result is a list of x coordinates (2D array), y coordinates (2D array), and water depths (3D array) at each x/y point." |
| 241 | + ] |
| 242 | + }, |
| 243 | + { |
| 244 | + "cell_type": "markdown", |
| 245 | + "id": "e9ddef7b", |
| 246 | + "metadata": {}, |
| 247 | + "source": [ |
| 248 | + "For use with FAModel/MoorPy, you can then write this information to a MoorPy/MoorDyn-readable text file:" |
| 249 | + ] |
| 250 | + }, |
| 251 | + { |
| 252 | + "cell_type": "code", |
| 253 | + "execution_count": 121, |
| 254 | + "id": "fb333bc3", |
| 255 | + "metadata": {}, |
| 256 | + "outputs": [], |
| 257 | + "source": [ |
| 258 | + "geo.writeBathymetryFile('bathymetry_humboldt.txt', bath_xs, bath_ys, bath_depths)" |
| 259 | + ] |
| 260 | + }, |
| 261 | + { |
| 262 | + "cell_type": "markdown", |
| 263 | + "id": "557d4f13", |
| 264 | + "metadata": {}, |
| 265 | + "source": [ |
| 266 | + "All of the above processes can be summarized in:" |
| 267 | + ] |
| 268 | + }, |
| 269 | + { |
| 270 | + "cell_type": "code", |
| 271 | + "execution_count": 122, |
| 272 | + "id": "518a9b21", |
| 273 | + "metadata": {}, |
| 274 | + "outputs": [ |
| 275 | + { |
| 276 | + "name": "stderr", |
| 277 | + "output_type": "stream", |
| 278 | + "text": [ |
| 279 | + "c:\\code\\floatingarraydesign\\famodel\\famodel\\geography.py:146: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n", |
| 280 | + "\n", |
| 281 | + " centroid = ( lease_area.geometry.centroid.values.x[0], lease_area.geometry.centroid.values.y[0] )\n" |
| 282 | + ] |
| 283 | + }, |
| 284 | + { |
| 285 | + "name": "stdout", |
| 286 | + "output_type": "stream", |
| 287 | + "text": [ |
| 288 | + "dict_keys(['lease_longs', 'lease_lats', 'lease_centroid', 'lease_xs', 'lease_ys', 'bath_longs', 'bath_lats', 'bath_xs', 'bath_ys', 'bath_depths'])\n" |
| 289 | + ] |
| 290 | + } |
| 291 | + ], |
| 292 | + "source": [ |
| 293 | + "info = geo.getLeaseAndBathymetryInfo(lease_name, gebco_file, bath_ncols=100, bath_nrows=100)\n", |
| 294 | + "print(info.keys())" |
| 295 | + ] |
| 296 | + }, |
| 297 | + { |
| 298 | + "cell_type": "markdown", |
| 299 | + "id": "b276baab", |
| 300 | + "metadata": {}, |
| 301 | + "source": [ |
| 302 | + "### Soil\n", |
| 303 | + "\n", |
| 304 | + "Similar processes can be done for extracting and setting up soil information\n", |
| 305 | + "\n", |
| 306 | + "```getSoilGrid()``` reads an input soil shape file and creates a rectangular grid of x coordinates (2D), y coordinates (2D) and soil names (3D string array), based on the shape file while also converting lat/long information to meters. It works by creating shape objects based on the data of the shapefile and using data within the shape file to extract the soil type name (e.g., 'mud', or 'rock')\n", |
| 307 | + "\n", |
| 308 | + "```getSoilType()``` returns the soil type name (string) that an input lat/long is above" |
| 309 | + ] |
| 310 | + }, |
| 311 | + { |
| 312 | + "cell_type": "markdown", |
| 313 | + "id": "ffd18aec", |
| 314 | + "metadata": {}, |
| 315 | + "source": [ |
| 316 | + "These are just the immediate needs for the related floating array modeling and design work. There are likely many more applications using geography-related tools that can be done." |
| 317 | + ] |
| 318 | + } |
| 319 | + ], |
| 320 | + "metadata": { |
| 321 | + "kernelspec": { |
| 322 | + "display_name": "famodel-env", |
| 323 | + "language": "python", |
| 324 | + "name": "python3" |
| 325 | + }, |
| 326 | + "language_info": { |
| 327 | + "codemirror_mode": { |
| 328 | + "name": "ipython", |
| 329 | + "version": 3 |
| 330 | + }, |
| 331 | + "file_extension": ".py", |
| 332 | + "mimetype": "text/x-python", |
| 333 | + "name": "python", |
| 334 | + "nbconvert_exporter": "python", |
| 335 | + "pygments_lexer": "ipython3", |
| 336 | + "version": "3.11.5" |
| 337 | + } |
| 338 | + }, |
| 339 | + "nbformat": 4, |
| 340 | + "nbformat_minor": 5 |
| 341 | +} |
0 commit comments