Part III Advanced areal weighted interpolation
Zhuotong Nan (南卓铜) ([email protected])
Motivation and objective
ArcView 9.2 comes with some basic interpolation approaches, for example, inverse distance weighted, Kringing, etc. In some cases, we want interpolate based on areal weighted. For example, Radar measurement can produce rainfall data over a study area. In order to use those Radar rainfall in hydrological model, e.g., VIC model, with a 1/8 degree spatial resolution, we need to do interpolation from a higher resolution, e.g., 4km by 4km. As shown in Figure 1, each grey cell has a measured value from Radar. We would like to generate an appropriate value for a VIC cell which might cover several Radar cells. In this part, we will illustrate an areal weighted approach.
Figure 1 Radar rainfall cells and 1/8 deg VIC cells
1. Create VIC grid polygons
Install Fishnet script with ArcView 9.2
ArcView 9.2 comes with a fishnet function (Data Management Tools → Feature Class → Create Fishnet) which is capable of generating a regular grid lines. However due to the limitation of ArcView license, it’s hard to implement conversion from polyline to polygon which is necessary to do areal weighted interpolation to be introduced in this section. Therefore, we will employ a free yet powerful script from ArcScripts website.
- Get this free script from http://arcscripts.esri.com/details.asp?dbid=12807. Unzip it to any location. We will see fishnet.dll in this location.
- Start ArcMap with a new document.
- Open Customize… dialog (Tools → Customize…), in Commands tab click Select from files… button to load fishnet.dll. You will see,
Figure 2 Customization dialog
Drag Create Fishnet to any place on toolbar. You will see on toolbar. Close Customize dialog.
Note, to uninstall the script, you can open Customize dialog again, and drag fishnet icon on toolbar to any place except toolbar.
Create 1/8 degree VIC grid polygons
VIC is a macroscale hydrologic model that solves full water and energy balances. It runs on basis of VIC grid. See http://ldas.gsfc.nasa.gov/LDAS8th/LDASspecs/LDASspecs.shtml for details of grid extent. We will generate VIC grid polygons using Fishnet utility download from ArcScripts in a scale of 1/8 degree.
- Open Create Fishnet utility (
).
- Configure parameters as following figure.
Figure 3 Parameters for creating 1/8 degree VIC grid
- Wait its completion. This might take a couple of minutes. A polygon shapefile “vicgrd8.shp” in geographic coordinate system (i.e., longitude and latitude, specified by GCS_North_American_1927 parameter) will be generated upon its completion.
Reproject VIC grid polygons to Albers equal-area projection system
- Run Feature project utility under Data Management Tools → Projections and Transformations → Feature → Project in ArcToolbox. Configure parameters as following,
Figure 4 Reproject to Albers Equal Area projection system
The produced vicgrd8_alb.shp now is in Albers Equal Area projection system. Note here we specified a transformation method to convert NAD_1927 to NAD_1983. Remember vicgrd8.shp is in NAD_1927 system and now vicgrd8_alb.shp is in NAD_1983 system.
- Close ArcMap without saving the current document.
2. Generate polygons for XMRG radar data file
- Get example data from http://www.nws.noaa.gov/oh/hrl/dmip/2/data/0703200407z.asc.gz. Unzip it (using WinRar, etc) to our current workspace (c:WorkspaceLab2).
- Convert the generated Grid in HRAP coordinate system to Stereogrpahic Polar coordinate system by using following equations. See this link http://www.nws.noaa.gov/oh/hrl/distmodel/hrap.htm for more details.
xster=hrapx*4762.5 – 401*4762.5
yster=hrapy*4762.5-1601*4762.5
Copy 0703200407z.asc to 0703200407z_ster.asc. Then open 0703200407z_ster.asc using any text editor, change its header to,
ncols 335
nrows 159
xllcorner -161925
yllcorner -6372225
cellsize 4762.5
Save 0703200407z_ster.asc.
- Convert 0703200407z_ster.asc to binary Grid file. In Command Line Window, run
ASCIIToRaster c:workspacelab2703200407z_ster.asc c:workspacelab2g0703200407z FLOAT
- Define Polar Stereographic Coordinate System for g0703200407z grid file.
In ArcCatalog, right click g0703200407z to open its Properties window. Navigate to Spatial Reference section, click Edit… to define its spatial reference. Choose Define the coordinate system interactively. Select Stereographic (Polar view) and configure it as shown below.
Figure 5 Parameters for Polar Ster projection system
- Convert raster to polygon
The conversion from raster to polygon only support grid file in integer. To do this, run following commands in Command Line Window.
Times_sa c:workspacelab2g0703200407z 100 c:workspacelab2g07032004_100
Int_sa c:workspacelab2g07032004_100 c:workspacelab2g07032004int
At this point, the produced g07032004int is already an integer grid file.
RasterToPolygon c:workspacelab2g07032004int c:workspacelab2g0703200407_poly NO_SIMPLIFY VALUE
A field “GRIDCODE” in the generated polygon (g0703200407_poly) is corresponding to cell values in g07032004int, which is 100 times real rainfall value.
- Project XMRG polygons to Albers projection system
Using Data Management Tools → Projections and Transformations → Feature → Project in ArcToolbox to project g0703200407_poly in Polar Stereographic coordinate system to g0703200407_alb in Albers Equal Area Coordinate System. This target coordinate system can be imported from vicgrd8_alb.shp.
3. Calculate area contribution to each VIC cell
Calculate area for each VIC cell
Using CalculateAreas command to calculate area for each VIC cell in Command Line Window.
CalculateAreas c:workspacelab2vicgrd8_alb.shp c:workspacelab2vicgrd8_alb_area.shp
The field “F_AREA” in vicgrd8_alb_area.shp indicates area in square meter.
Overlay
Intersect utility is used to calculate areal contribution from XMRG polygons to each VIC cell.
- First we will clip the whole VIC grid shapefile with XMRG polygons extent. In Command Line Window,
Workspace c:workspacelab2
Clip_analysis vicgrd8_alb_area g0703200407_alb c:workspacelab2vicgrd8_clp.shp
- In Arctoolbox, run Analysis Tools → Overlay → Intersect with input parameters specified as below.
Figure 6 Intersect input parameters
- Start ArcMap with a new document. Add intersect.shp. Save map document into c:workspacelab2partIII.mxd.
- Open intersect.shp attribute table.
Figure 7 Attribute table of intersect.shp
Here FID_vicgrd and F_AREA come from VIC cells, while FID_g07032 and GRIDCODE come from XMRG polygon file. FID_vicgrd uniquely represents each VIC cell whose area is show in F_AREA field. GRIDCODE is actually equal to 100 * rainfall.
Querying attribute, for example FID_vicgrd = 1400, will extract all small polygons consisting a VIC cell.
Calculate areal weighted rainfall
- Add a new field VIC_AREA with values come from F_AREA.
Open intersect.shp’s attribute table, Click Options → Add Field…
Figure 8 Add a new field
Click the header of VIC_AREA column to select this field, right click it then select Field Calculator…, Yes upon warning message. Input [F_AREA] in the lower textbox then click OK to start calculation.
Figure 9 Field calculator
- Calculate area by running CalculateAreas in Command Line Window.
CalculateAreas intersect c:workspacelab2inters_area
- Add a new field named “percentage” in double type to inters_area.shp.
- Calculate [percentage] = [F_AREA] * [GRIDCODE]/100 / [VIC_AREA].
Figure 10 Calculate Percentage field
- Export areal weighted rainfalls upon each VIC cell
Make inters_area layer visible, and make sure there is no selection in this layer. Otherwise use on the Tools toolbar to clear any existing selection. In ArcToolbox, run Analysis Tools → Statistics → Summary Statistics to calculate rainfalls cell by cell.
Figure 11 Summary statistics
Link interpolated rainfall data to VIC cells
- Make vicgrid8_clp.shp visible and be topmost. Link this layer’s attribute table to rainfalls.dbf by using context menu → Joins and Relates → Join. In Join Data window, enter inputs as shown below,
Figure 12 Join to rainfalls.dbf
- Export vicgrd8_clp to vicgrd8_rainfall.shp. The output shapefile will physically contain rainfall values (stored in field SUM_percen).
- Save partIII.mxd.
- Project vicgrd8_rainfall to rainfall_geo.shp.
- Add an alias “Rainfall” to field “SUM_Percen”. This makes the field name more straightforward.
Figure 13 Make field name self-explanatory
In order to use in combination with other raster, we can also convert this polygon shapefile to raster using FeatureToRaster command in Command Line Window.
4. Visualization
- Open partIII.mxd, switch all layers off except rainfall_geo. Set current data framework coordinate system to GCS_North_American_1983 using menu View→Data framework properties…→Coordinate System.
- Open rainfall_geo layer properties window, make symbology settings same as below. You can also choose a different color schema for this.
Figure 20 Symbology
- Export map for submission. See Figure 21 for reference.
Figure 21 Exported map
Conclusion
Although ArcView 9.2 has already included some basic interpolation approaches, it turns out in many cases hard to meet our specific requirements. This handout walks through how to interpolate Radar rainfall data to 1/8 degree VIC grid, illustrating necessary steps to implement areal weighted interpolation approach using ArcView. This handout also show us essential GIS operations such as project, clip, intersect, export, manipulate attribute table, link feature class to external table, etc, as well as making use of third party script.
(转载请保持信息的完整，并加以适当引用)