Create a NetCDF file

In this lesson, you’ll use the Ecological Marine Units (EMU) feature class of the Caribbean Sea to create a NetCDF file using Python in ArcGIS Pro. First, you’ll set up a new project from a project package. Second, you’ll start a new Python notebook. Third, you’ll create an empty NetCDF file with the EMU schema. Lastly, you’ll load the EMU values into the NetCDF file. The output EMU NetCDF file will have the required format to be loaded and visualized as a voxel layer in ArcGIS Pro.

Set up a new project

In this section, you’ll start a new ArcGIS Pro project using a project package. This package contains the country and a sample of EMU data. First, you'll process this EMU data using an ArcGIS Notebook.

  1. Download the project package and double-click it to open it.

    EMU project package

    A new ArcGIS Pro project starts with the Caribbean Ecological Marine Units local scene and the EMU_Caribbean_Sea 3D points feature class.

    Note:

    You can turn on the World Ocean Base layer and the WorldElevation3D/TopoBathy3D layer (and turn off the Countries layer) for better geographic reference. The less data on the map, the faster it will draw.

  2. On the ribbon, click the Insert tab and choose New Notebook.

    Insert a new notebook.

    A new notebook is added to the project. Notebook details can be edited in the Catalog pane.

  3. In the Catalog pane, expand the Notebooks folder, and right-click New Notebook.ipynb.
  4. Choose Rename and type Create EMU NetCDF file.ipynb. Press Enter.

    Rename the notebook.

A new ArcGIS Pro project with a local scene and an empty Python notebook are ready.

Create an empty NetCDF file

First, you’ll create a NetCDF file for the EMU in the Caribbean Sea. The NetCDF file will be built in a latitude, longitude, and depth grid. The file will also include global attributes, dimensions, and variables. These are the foundational attributes needed to map (latitude, longitude, and depth) and display the associated variables (temperature, salinity, and oxygen) in 3D.

Multidimensional data is collected over space and time, depth or height.

  1. In the notebook’s empty cell, type the following commands to load the required Python modules:
    import os
    import arcpy
    import netCDF4
    import numpy as np
    import datetime as dt

    These modules are the libraries that contain the functions used in the script. For example, the arcpy module has the ArcGIS geoprocessing, data management, and mapping functions. The netCDF4 module allows us to read and write NetCDF files.

  2. On the Notebook ribbon, click Run.

    Run the first cell of the notebook.

    Once the first cell runs, the braces will show the number 1. A new cell is added below the first.

    Note:

    You can add new cells from the notebook ribbon by clicking the Insert tab and choosing Insert Cell Below.

  3. In the new cell, paste the following commands related to the output NetCDF file and the regular 3D grid and click Run:
    # Output file
    output_nc = "EMU_Caribbean_Voxel.nc"
    
    # Grid
    lowerleft_x = -88.625
    lowerleft_y = 8.375
    cellsize = 0.25
    
    # Number of cells
    no_rows = 58
    no_cols = 117
    no_depths = 102
    
    # Values
    lat_vals = np.arange(lowerleft_y, lowerleft_y + cellsize*no_rows, cellsize)
    lon_vals = np.arange(lowerleft_x, lowerleft_x + cellsize*no_cols, cellsize)
    z_vals = np.concatenate([np.arange(-5500, -2000, 100), np.arange(-2000, -500, 50),
                             np.arange(-500, -100, 25), np.arange(-100, 5, 5)])
    
    # Print vectors
    print('Latitude')
    print(lat_vals)
    print('Longitude')
    print(lon_vals)
    print('Depth')
    print(z_vals)

    Each variable definition is broken up with associated comment (#). Output_nc specifies the name of the output NetCDF file. Lowerleft_x and lowerleft_y specify the coordinates for the NetCDF file. Cell size represents the resolution of 1/4 degree. No_rows specifies the number of rows, which is affiliated with latitude; in this case, 58 rows at 0.25 cell size means 14.5 degrees of latitude. No_cols specifies the number of columns, which is affiliated with longitude; in this case, 117 columns at 0.25 cell size means 29.25 degrees of longitude. This means that the output will be a little more than twice as wide (longitude) as it is long (latitude). No_depths specifies the number of depth values that will be included in the NetCDF file. This represents our water column space in 3D. In the next variable definition (Values), you specify how these can all be used together to produce a regularly spaced horizontal grid and irregularly spaced vertical (depth) grid.

    The latitude, longitude, and depth vectors are printed on the screen.

  4. In the empty cell, paste the following code using the Dataset function of the netCDF4 Python module to create an empty NetCDF file. Provide global attributes to add documentation to the file and run the cell.
    # Create File
    nc_file = netCDF4.Dataset(output_nc, 'w', format="NETCDF4", set_auto_mask=False)
    
    # Global Attributes
    nc_file.title = 'Ecological Marine Units'
    nc_file.summary = 'Ecological Marine Units (EMUs) for the Caribbean Sea.'
    nc_file.keywords = 'Global, EMU, USGS, NOAA, Esri, 3D, Ocean, Marine, Ecology'
    nc_file.license = ('This work is licensed under a Creative Commons Attribution 4.0 '
                       'International License.')
    nc_file.references = ('Please see refrences Cited in this publication: ' 
                        'http://www.aag.org/galleries/default-file/AAG_Marine_Ecosyst_bklt72.pdf')
    nc_file.source = 'Esri, USGS,  NatureServe, GEO, NOAA, NASA, NIWA'
    nc_file.Conventions = 'CF-1.6'
    nc_file.institution = 'Esri'
    nc_file.history = ('{0} Creation of EMU Caribbean '
                       'netcdf file.').format(dt.datetime.now().strftime("%Y-%m-%d"))

    Note:

    It is fundamental to document a NetCDF file so other people can understand it better. You can follow the Climate and Forecast (CF) Metadata Conventions to guide you on standard attributes and variable names.

  5. In a new cell, paste the latitude, longitude, and depth dimensions.
    # Create Dimensions
    lat_dim = nc_file.createDimension('latitude', no_rows)
    lon_dim = nc_file.createDimension('longitude', no_cols)
    top_dim = nc_file.createDimension('depth', no_depths)
  6. Add the following code for the mapping and gridding variables: crs, latitude, longitude, and depth. Load the lat_vals, lon_vals, and z_vals vectors into their respective (in other words, lat_var, lon_var, and depth_var) variables and run the cell.
    # Create the map and grid variables
    crs_var = nc_file.createVariable('crs', 'l', (), fill_value=-9999, zlib=True)
    crs_var.standard_name = 'crs'
    crs_var.grid_mapping_name = 'latitude_longitude'
    crs_var.crs_wkt = arcpy.SpatialReference(4326).exportToString()
    
    lat_var = nc_file.createVariable('latitude', 'f4', ('latitude'), fill_value=-9999)
    lat_var.units = 'degrees_north'
    lat_var.standard_name = 'latitude'
    lat_var.axis = 'Y'
    
    lon_var = nc_file.createVariable('longitude', 'f4', ('longitude'), fill_value=-9999)
    lon_var.units = 'degrees_east'
    lon_var.standard_name = 'longitude'
    lon_var.axis = 'X'
    
    depth_var = nc_file.createVariable('depth', 'l', ('depth'), fill_value=-9999)
    depth_var.short_name = 'depth'
    depth_var.standard_name = 'depth'
    depth_var.positive = 'up'
    depth_var.units = 'meters'
    depth_var.axis = 'Z'
    
    # Load data
    lat_var[:] = lat_vals
    lon_var[:] = lon_vals
    depth_var[:] = z_vals

    Notice that the crs variable does not have any dimensions or variables; it provides the mapping spatial reference information.

  7. In a new cell, paste three variables for the EMU data: sea water temperature, salinity, dissolved oxygen, and cluster37, and run the cell.
    # Create EMU Variables
    temp_var = nc_file.createVariable('temp', 'f4',
                                      ('depth', 'latitude', 'longitude'), fill_value=-9999)
    temp_var.short_name = 'temp'
    temp_var.standard_name = 'sea_water_temperature'
    temp_var.units = 'degree_Celsius'
    temp_var.grid_mapping = 'crs'
    
    sali_var = nc_file.createVariable('salinity', 'f4',
                                      ('depth', 'latitude', 'longitude'), fill_value=-9999)
    sali_var.short_name = 'salinity'
    sali_var.long_name = 'sea_water_salinity'
    sali_var.units = ''
    sali_var.grid_mapping = 'crs'
    
    dissO2_var = nc_file.createVariable('dissO2', 'f4',
                                      ('depth', 'latitude', 'longitude'), fill_value=-9999)
    dissO2_var.short_name = 'dissO2'
    dissO2_var.long_name = 'mole_concentration_of_dissolved_molecular_oxygen_in_sea_water'
    dissO2_var.units = 'microliters/liter'
    dissO2_var.grid_mapping = 'crs'
    
    cluster37_var = nc_file.createVariable('cluster37', 'l',
                                           ('depth', 'latitude', 'longitude'), fill_value=-9999)
    cluster37_var.short_name = 'cluster37'
    cluster37_var.long_name = 'cluster37_id'
    cluster37_var.units = ''
    cluster37_var.grid_mapping = 'crs'

Sea water temperature, salinity, and dissolved oxygen are key variables used for the characterization of the ocean and its environmental condition. The cluster variable holds the classification of a grid cell within the EMU dataset. Notice that the first three variables have a floating number data type (in other words, ‘f4’) while the cluster37 variable has an integer number data type (in other words, ‘l’).

Note: Visit the Ecological Marine Units (EMU) web app in ArcGIS Living Atlas of the World or read the American Association of Geographers publication to learn more about the EMU dataset. The EMU dataset contains more variables than the three listed in this exercise. You can download the full EMU NetCDF file from ArcGIS Online.

Load data into the NetCDF file

In this section, you’ll load the values from the EMU point feature class into the NetCDF file.

  1. In the Contents pane, right-click the EMU_Caribbean_Sea feature class and choose Attribute Table.

    Open the attribute table.

    The temp, salinity, dissO2, and Cluster37 fields for each 3D point (coordinates X, Y, and Z) contain the information that will be loaded into the NetCDF file variables.

    View Ecological Marine Unit attributes.

  2. Close the attribute table.
  3. In the last cell of the notebook, assign the EMU_Caribbean_Sea feature class into the fc variable. Create a Search Cursor with the fc variable and the fields: temp, salinity, dissO2, Cluster37, and Shape@. Do not run the cell yet.
    # Load EMU data
    fc = os.path.abspath(r'.\EMU.gdb\EMU_Caribbean_Sea')
    cursor = arcpy.da.SearchCursor(fc, ['temp', 'salinity', 'dissO2', 'Cluster37', 'Shape@'])
  4. In the Contents pane, double-click the EMU_Caribbean_Sea layer. In the Layer Properties window, click the Source tab.

    Find the source path to the geodatabase.

  5. Copy the source path from the Database row, and paste it into the notebook as the fc path variable. Make sure to include the file name EMU_Caribbean_Sea at the end.

    The complete variable fath for fc

  6. In the same cell, use a for loop to iterate over the rows in the cursor. Get the EMU variable values, replace null values with -9999, get the array indices for the 3D point, and store the variable values in the NetCDF array.
    # For each record in the table
    for row in cursor:
        # Get values
        temp, salinity, dissO2, cluster37 = row[:4]
        coords = row[4].getPart()
        # Replace null values with -9999
        temp = float(temp) if temp else -9999
        salinity = float(salinity) if salinity else -9999
        dissO2 = float(dissO2) if dissO2 else -9999
        cluster37 = float(cluster37) if cluster37 else -9999
        # Get indices in array of current 3D point
        idx = (np.abs(lon_vals - coords.X)).argmin()
        idy = (np.abs(lat_vals - coords.Y)).argmin()
        idz = np.where(z_vals==int(coords.Z))[0]
        # Store variables values
        temp_var[idz, idy, idx] = temp
        sali_var[idz, idy, idx] = salinity
        dissO2_var[idz, idy, idx] = dissO2
        cluster37_var[idz, idy, idx] = cluster37
    del(cursor)

    For loop that will iterate over the rows in the cursor

  7. Run the cell.

    The cell may take some time to run.

  8. In a new cell, use nc_file.close() to close the NetCDF file and run the cell.
    # Close NetCDF file
    nc_file.close()

    Note:

    You can download the final Python notebook from ArcGIS Online.

  9. On the ribbon, on the Notebook tab, click Save.

In this lesson, you created an empty NetCDF file. Then, you loaded values into the NetCDF file from a point feature class using Python Notebooks in ArcGIS Pro. In the next lesson, you’ll load this NetCDF file as a multidimensional voxel layer in a local scene.


Add a NetCDF file as a voxel layer

In the previous module, you created a NetCDF file using Python notebooks in ArcGIS Pro and the EMU point feature class. In this module, you’ll load the EMU NetCDF file in a local scene as a multidimensional voxel layer. The voxel layer can be used to visualize EMU variables such as sea water temperature, salinity, or dissolved oxygen in 3D.

Add a multidimensional voxel layer

In this section, you’ll add the EMU NetCDF file as a voxel layer.

  1. Click the Caribbean Ecological Marine Units tab to open the local scene.
  2. On the ribbon, on the Map tab, click the Add Data arrow. In the menu, choose Multidimensional Voxel Layer.

    Add multidimensional voxel layer.

  3. In the Add Voxel Layer window, click Browse and choose the EMU_Caribbean_Voxel.nc file as the Input Data Source.

    Once you choose the file, the available attribute fields will populate the Select Variables section.

  4. For Select Variables, in the Default Variable column, set dissO2 as the default.

    Add Voxel Layer window

  5. Click OK.

    Each variable defined in the NetCDF file is listed in the Add Voxel Layer window. You can decide which variables to include in the voxel layer by checking the variables on or off. You can choose which variable to set as the default variable. When the voxel layer gets added to the scene, this will be the variable shown.

    A variable has a defined data type that determines how the variable can be symbolized. For example, a continuous data type could be numeric values such as sea water temperature, salinity, or dissolved oxygen. The continuous variables can be visualized using stretch values. In cases in which the data type is discrete like cluster37, you can symbolize the variable using unique values.

    The EMU NetCDF file appears as a voxel layer in ArcGIS Pro. It's currently drawn with default colors and elevation values. To more accurately visualize this data, you'll change the elevation and styling.

    EMU voxel layer

    This layer is currently drawing on top of the ocean. To correctly place it at depth, you'll set the vertical exaggeration.

  6. In the Contents pane, click the EMU_Caribbean_Sea_Voxel layer to select it.
  7. On the ribbon, in the Voxel Layer contextual menu, click the Appearance tab. In the Elevation group, set Offset to -200,000 m and Vertical Exaggeration to 3,000 m.

    Elevation parameters

    Next, for visualization purposes, you'll apply a style developed specifically for showing oceanographic data. The oceanographic style includes a variety of color schemes designed to visualize characteristics such as oxygen level in the ocean. The color range and label need to match to make it easy for your audience to interpret the values in the voxel layer.

  8. Download the cmocean.stylx style package.
  9. On the ribbon, click the Insert tab. In the Styles group, click Add and choose Add Style.

    Add Style

  10. In the Add a style file window, choose the cmocean.stylx file and click OK.

    Now you can apply the style to the layer.

  11. On the ribbon, click the Voxel Layer contextual Appearance tab. In the Drawing group, click Symbology.

    The Symbology pane appears.

  12. In the Symbology pane, for Color Scheme, expand the menu and check Show names and Show all.

    Show all color schemes and their names.

    The ocean styles you downloaded are shown at the top of the list.

  13. Find and select the oxy color scheme.

    Choose color scheme.

    The map redraws to show the new color scheme. A low concentration of dissolved oxygen is shown in reds and maroons, while a high concentration is shown in yellows.

    Dissolved O2 voxel layer

    The middle values are shown in gray. To emphasize high and low values, you'll use transparency to hide these intermediate values.

  14. In the Symbology pane, for Label and Color range, set the Minimum value to 3.00 and Maximum value to 5.00.

    Set the Label and Color range values.

  15. Check the Transparency function box.

    Define the transparent range by adding control points.

  16. Click Add control point twice.

    Two additional points are added to the symbol range.

  17. Click the first control point and set Value to 100 percent and Position to 3 percent.

    First control point

  18. For the second control point, set Value to 100 percent and Position to 80 percent.

The voxel layer shows areas of high (yellow) and low (red) dissolved oxygen in the Caribbean Sea.

Finished symbology

In this lesson, you added the EMU NetCDF file as a multidimensional voxel layer in a local scene in ArcGIS Pro. In the next lesson, you’ll visually analyze the voxel layer by defining surfaces and making basic changes to the voxel layer appearance such as selecting a variable and changing stretch parameters.


Visually analyze the EMUs in 3D

In the previous module, you loaded the EMU NetCDF file as a voxel layer in ArcGIS Pro. In this module, you’ll create surfaces such as isosurfaces and sections to visually analyze variables.

Create an isosurface

In this section, you will define an isosurface at a specific value to visually analyze the distribution of dissolved oxygen. An isosurface describes a surface within the voxel layer along a specified value. For example, an isosurface of dissolved oxygen at 5 microliters/liter shows all the places in the Caribbean Sea with that value which can happen at different depths.

  1. In the Contents pane, under EMU_Caribbean_Voxel, click mole_concentration_of_dissolved_molecular_oxygen_in_sea_water two times.

    The field becomes editable.

  2. Name the variable Concentration of Dissolved Oxygen.

    Rename the dissolved oxygen surface variable.

  3. Click Surfaces and expand it.

    The voxel will be removed from the map since you haven't created a surface yet.

  4. Under Surfaces, right-click Isosurfaces and select Create Isosurface.

    Create an isosurface.

    A new isosurface is created within the voxel layer, and the Voxel Exploration pane appears.

  5. In the Voxel Exploration pane, for Name, type Concentration of Dissolved Oxygen 5. For Value, type 5.

    Voxel Exploration pane

    An isosurface at 5 microliters/liter of dissolved oxygen is drawn on the map.

  6. For Color, click the default color and choose a dark purple color.

    When drawn with other layers, such as the basemap, the dark purple will make the threshold value more prominent.

  7. In the Contents pane, turn on the World Ocean Base and WorldElevation3D/TopoBathy3D layers.

    Concentration of dissolved oxygen at a value of 5 microliters/liter

  8. On the ribbon, click the Voxel Layer Appearance tab. In the Lighting group, for Diffuse, choose 70 percent and for Specular, choose 20 percent.

Voxel layer symbology

By changing the Diffuse and Specular parameters, the 3D effect of the isosurface is better represented on the local scene.

Create a section

In the previous section, you defined an isosurface to visualize a specific value within a variable to visualize the same variable conditions at different places and depths within the voxel layer. In this section, you'll visually explore variables by defining sections that span across the voxel layers. You'll define a section, then dynamically change its position, tilt, and orientation. The new section will be locked to be visible together with other variables.

  1. On the ribbon, on the Voxel Layer Appearance tab, in the Variable group, choose sea_water_temperature.

    Visualize the sea water temperature variable.

  2. In the Drawing group, click Symbology. In the Symbology pane, for Color scheme, check Show Names and choose Thermal.

    Choose thermal symbology.

  3. In the Contents pane, uncheck the Countries layer. Make sure the EMU_Caribbean_Voxel layer is selected with Surfaces chosen.
  4. On the ribbon, click the Voxel Layer Data contextual tab. In the Explore group, choose Slice and Section.

    The Slice and Section toolbar appears at the bottom of the Local Scene.

  5. In the Slice and Section toolbar, select Horizontal Section.

    Horizontal slice

  6. On the map, click anywhere within the volume to create the section.

    A new horizontal section is added to the Local Scene, and the Push and Pull tool is automatically selected.

  7. Click and drag the tool to move the horizontal slice up and down to explore the sea water temperature at different depths.

    Navigate through the horizontal slice.

  8. On the Slice and Section toolbar, choose the Tilt tool. Click into the scene and move the pointer to the left to tilt the section.

    Side tilt profile

  9. Tilt the section until it is vertical.

    Section tilted vertically.

    You can see the temperature profile from the surface to the deepest parts of the Caribbean Sea.

  10. On the Slice and Section toolbar, select the Orient tool. Click into the scene and move the pointer to the right to rotate the sections.

    Side section

    Tip:
    To better navigate the scene, click the Map tab and use the tools in the Navigate group.

  11. Find an interesting spot, such as a transversal section from Cuba to Panama to lock to your scene.
  12. In the Contents pane, right-click Section 1. Choose Lock Section.

    Lock the section for later visualization.

    The locked section moves into the Locked Sections group of the voxel layer. The section is now a snapshot at the defined location and can be viewed even if a different variable is selected. This allows you to explore different variables together.

  13. On the ribbon, on the Voxel LayerAppearance tab, in the Variable group, choose the mole_concentration_of_dissolved_molecular_oxygen_in_sea_water variable.

    The Concentration of Dissolved Oxygen 5 isosurface and the sea water temperature section appear on the Local Scene.

    Visual of a section and isosurface

    The image above shows a Sea Water Temperature cross section from Cuba to Panama. It also shows an isosurface of Dissolved Oxygen for values of 5 microliters/liter.

In this lesson, you created a NetCDF file using Python notebooks in ArcGIS Pro. You downloaded the EMU point feature class and loaded its values into the new EMU NetCDF file. You loaded and visually analyzed this NetCDF file in 3D as a multidimensional voxel layer in a local scene.