Embedded Files
This tutorial will show you how to quantify areal distortions for any type of projection using Nicolas Tissot's ellipse of distortion--also known as an indicatrix. While Tissot's indicatrix was originally designed to quantify both areal and conformal distortions, this tutorial will focus solely on quantifying areal distortions.
This tutorial was created using ArcGIS version 3.3.
Before starting this tutorial, you will need to download and install a dataset following these instructions:
Create a folder called indicatrix somewhere in your course working directory (e.g. C:\Users\yourname\Documents\ArcGIS_tutorials\indicatrix\).
Download the data for this exercise then extract the contents of indicatrix.zip into your newly created indicatrix folder.
Open the map document
Open the map document
Open the Tissot.aprx map document.
Identifying a layer's coordinate system
Identifying a layer's coordinate system
To identify a data layer's coordinate system, bring up it's Properties window by double-clicking on the layer name in the Contents pane.
Select Source from the left-hand tabs and expand Spatial Reference.
The layer's property indicates that the World polygons are stored in a Miller Cylindrical projected coordinate system (PCS). It also indicates that the PCS is built on top of a WGS 1984 geographic coordinate system (GCS). Recall that all PCS' are built off of a GCS ,hence why you will always see a GCS referenced alongside a PCS.
In the next steps, we will build Tissot indicatrix features to evaluate this coordinate system's spatial properties.
Create a point grid
Create a point grid
We will first need to create a point layer about which tissot circles will be generated. This point layer can be created manually (this may be useful if you want to measure distortions at specific locations), or it can be created in a semi-automated fashion following a gridded pattern. In this exercise, we will opt for the latter approach using ArcGIS' built-in fishnet tool.
If the geoprocessing pane is not open in your document, click on the Tools icon in the Analysis tab.
In the geoprocessing pane type fishnet in the search window.
One of the first results should be the Create Fishnet tool. Click on this tool to bring up its window.
The fishnet tool is designed to create a grid of polygons or polylines along with a point layer depicting the grid cells center. Our interest is in the point layer (we will discard the polygon or polyline geometries later).
In the Create Fishnet pane, name the output grid fishnet. Save it as a shapefile.
Next, you will need to set the extent of this new layer. We will use the World layer's extent.
Under Template Extent, click on the Extent of a Layer pulldown menu and select World.
Once selected, the World layer's extent boundaries should appear under the Template Extent rubric.
In the previous step, setting the Template Extent to the World layer will ensure that the points are confined to the feature's extent. However, we would like the points to be evenly distribute north and south of the equator. Since the World layer's extent does not go beyond the existing features (i.e. there is no land mass around the north pole), we will manually modify the extents such that both northern and southern hemispheres are evenly represented. This will require bumping the upper boundary to 14338757 m to match the -14338757 m of the lower boundary.
Change the upper boundary limit to 14338757 (or the absolute value of whatever lower boundary value the software generated from your layer).
The fishnet tool allows us to define the grid intervals in one of two ways: by cell height and width, or by number of rows and columns. We'll opt for the latter.
Jump to the Number of Rows and Number of Columns fields and type 6 for the former and 10 for the latter.
Make sure that the Create Label Points is checked (this is the layer we are after).
You can leave the Geometry Type set to Polyline (the geometry type does not matter since this layer will be removed from our project once created).
Click Run.
You should now have two new layers in your Contents pane: a fishnet layer and a fishnet_label layer.
Remove the fishnet polyline layer from the Contents pane by right-clicking it and selecting Remove from the pull-down option. You can also delete it from your project folder since it will not be used in this exercise.
Create indicatrix ellipses
Create indicatrix ellipses
When converting spatial features from a geographic coordinate system (GSC) to a projected coordinate system (PCS) one or more spatial properties may be distorted in the transformation. These properties include area, shape, distance and direction .
If you want to assess a projection’s spatial distortion across your study region, you can generate Tissot indicatrix (TI) ellipses. The idea is to project a small circle (i.e. small enough so that the distortion remains relatively uniform across the circle’s extent) and to measure its distorted shape on the projected map.
Navigate back to the Geoprocessing pane.
If the Create Fishnet pane is still active, click on its back arrow to navigate back to the search window.
In the search field, type Buffer.
Select Buffer (analysis Tools) from the suggested tools list.
Set fishnet_label as the Input Features.
Name the output polygon Indicatrix. You can save it as a shapefile or a geodatase file (but make sure to save it somewhere in your indicatrix/. folder).
Set the Linear distance to 300 US Survey miles (this will be the Tissot circle's radius). You may need to tweak this radius value to ensure that the resulting circles are big enough to properly compute area values, but not too big so as to have many overlapping circles (though this may not be possible for poleward regions).
Set the Method to Geodesic (this is important!).
Click Run.
By adopting a geodesic approach to generate a circle about each point, we have ArcGIS generate the circles on an unprojected coordinate system (i.e. on the ellipsoid representation of the earth surface). This guarantees that the circles are "true" to form and not distorted by any projection used in the layer or map frame.
A true indicatrix quantifies the distortion in area and angles at a point location (or an infinitesimally small area). However, doing so in an ArcGIS environment would generate circles not viewable to the user. The trick is to pick a radius just big enough for the user to make out any spatial distortion while minimizing the influence a non-uniform spatial distortion distribution (inherent in most projections) can have on the outer bounds of the indicatrix circles.
Your output should look something like this. Note the varying degree in shape distortion as one moves poleward.
Checking that these are indeed circles
Checking that these are indeed circles
So how can we be sure that the tool generated conformal circles? One way to address this question visually is to look at the indicatrix layer on a 3D globe. We can add a new 3D map to our current ArcGIS project.
Click on the Insert tab.
Click on New Map and select New Global Scene.
This will add a new map tab to your current document. You can toggle back and forth between the Map frame and the 3D Scene frame by selecting the appropriate tab. In doing so, you will notice that each map has its own unique set of layers (i.e. they do not share the same Contents pane). You will therefore need to add the indicatrix layer to the 3D scene's Contents pane. You can add the layer via the Add Item button under the Insert tab, or you can copy the layer from the Map frame to the 3D scene window. We'll opt for the latter approach in the next step.
Switch to the Map frame, right-click the indicatrix layer and select Copy.
Switch back the the 3D Scene frame, right-click Scene, then select Paste.
Feel free to spin the 3D globe around. Do the tissot circles appear to be true to form?
This is how we expect to view the tissot circles if the projection preserves area and shape. Given that we know that the tissot circles are conformal, observing their shape on a projected surface will inform us on that projection's spatial properties.
Calculating areal distortion
Calculating areal distortion
In this next step, we will quantify the areal distortion afforded by the Miller projection. We will do so by calculating both the geodesic area of the circle (this will be the "true" area) and the projected area as defined by the Miller projection, then compute their difference as a percentage.
Switch to the Map frame.
Right-click the indicatrix layer and select Attribute Table.
NOTE:
If you saved the indicatrix layer in a geodatabase, you will see a Shape_Length and a Shape_Area field. These are automatically updated and their units are those of the Miller coordinate system. We will not use these columns, instead, we will compute the area fields from scratch.
In the Attribute table, click on the Add Field icon in the upper-left hand corner of the table pane.
This will bring up the Fields pane. This is where you can add, delete and define new fields.
Name the field Area_mil (this will be the field that will store the projected area measurements).
Make it a Float data type.
Click on "Click here to add a new field" at the bottom of the Fields pane.
Name the second new field Area_geo and make it a Float data type.
In the Fields tab, click on the Save button.
Note that it is important to save any changes made to the attribute table fields. Not doing so makes these fields unavailable in the attribute data table.
Navigate back to the data table.
Right-click on the Area_mil field and select Calculate Geometry.
Note: If you don't see the Area_mil field, you probably did not save changes as outlined in the previous step.
Select Indicatrix for Input Features if this is not already set.
In the Target Field, assign the property Area to the Area_mil field.
Note that this tool allows for multiple field calculations. We'll have it calculate the geodesic area too.
Select Area_geo from the next pull-down field and assign the Area (geodesic) property.
Set the area unit to square US survey miles (this will be applied to both area measurements).
Finally, we can specify the coordinate system to use in calculating the planar (projected) area. If this field is not explicitly defined, it will adopt the feature layer's CS. Just to be safe, we'll select World from the pull-down menu. This will set the coordinate system to this layer's Miller Cylindrical coordinate system. Note that the internal name for this PCS (World_Miller_Cylindrical) is a bit different than that shown in the Properties window.
Click OK.
Next, we'll compute the percent difference between both area measurements.
Note: If the calculations were only applied to just a handful of rows, you probably had indicatrix features selected. To ensure that the calculations apply to all records, clear any selection you might have in your map view and rerun the geoprocess.
Following the steps outlined earlier in this section, add a new field called diff_mil and set its data type to Float. (Don't forget to save these changes before returning to the data table).
Right-click on the newly created diff_mil field in the indicatrix attribute table and select Calculate Field.
Make sure that the indicatrix layer and the diff_mil field are populating the Input Table and Field Name fields.
You can use either a Python syntax or an Arcade syntax when defining an expression. Since we've used Python syntax in previous tutorials, we'll switch to Arcade syntax for a change (but using Python would be perfectly fine).
Change the expression type to Arcade.
In the Expression box, type the expression:
($feature.Area_mil - $feature.Area_geo) / $feature.Area_geo * 100
You can double-click on the field names from the list of Fields to have them automatically populate the expression box. This is helpful in minimizing typos.
Also, make sure that none of the indicatrix features are selected in your map before running the geoprocess.
Click OK to run the geoprocess.
Feel free to symbolize the indicatrix features using the diff_mil field. Refer to the Symbolizing Features tutorial for guidance. What would be an appropriate color scheme to adopt?
A positive diff_mil value indicates that the projected coordinate system over-estimates the true area value. A negative value indicates that the projected coordinate system under-estimates the true area value.
In the adjoining figure, a manually defined classification scheme with four breaks is used to symbolize the circles. Note the increasing distortion in area of the indicatrix circles as one moves poleward. The maximum distortion is 1462% of the true area value!
In addition to the observed distortion in shape (note the flattened circles close to the poles), the error in area increases as one approaches the poles. This is expected given the nature of this cylindrical coordinate system.
A note about Mercator projections
A note about Mercator projections
As of version 2.4.2 of ArcGIS Pro, some projected coordinate systems, such as Mercator, do not offer the planar area option in the Calculate Geometry tool (See BUG-000122213 ). If you want to measure the aerial distortion associated with a Mercator projection, you can adopt this workaround:
Use the Calculate Field tool instead of the Calculate Geometry tool to compute the attribute value.
Type the following Python expression in the Calculate Field tool instead of the Calculate Geometry tool:
!Shape!.getArea("planar", "squaremiles")
This should get you the planar area for the desired projection. Other units that can be computed are: "SQUAREKILOMETERS", "ACRES" , "HECTARES" , "SQUAREMETERS".
This wraps up this in-class lab activity.
Page updated
Report abuse