Interpolate values along a line from a surface in ArcGIS

This is kinda stupid and I’m sure there are much more elegant ways to do this. But for a quick and dirty job it’s good enough to me.

Anyways, what I want to do is to draw a line on a surface (could be anything, DEM, velocity field from CFD, or precipitation map) and extract the values along this line then output as a table. It can be easily done using ArcGIS 3D Analyst toolbar then you get something like this:

Capture

By clicking right mouse this data can be saved in a file and job done.

However, I just found that if you need to output data along the same line for different surfaces, this toolbar thing could become very annoying and you may need to draw multiple lines manually, which make them not exactly the same.

So here is another way to do this. First create a polyline, then in “Editor” tool, split it to equal length segments. Second, in 3D Analyst Tools -> Functional Surface -> Add Surface Information, select the surface and polyline just created, and choose Z_MEAN as output property. Run it, then you get the mean value of each segment in the attribute table.

Final trick, in attribute table, select some features, then right click the very left column and hold “shift” at the same time, there’s a copy option to allow you paste the table in Excel.

Advertisements

ArcGIS shape file plot in Matlab

ArcGIS shapes such as points and polygon, or even raster files can be plotted in Matlab. Here is an example of reading and plotting HUC4 polygons in Matlab.

First if read in the shape file that’s generated in ArcGIS

%read in shapefile
huc = shaperead('C:\HUC4.shp','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);

then we get a struct like this, in which most of the columns are attribute table variables, with the first few columns showing the geometrical information about each polygon.

Capture

Sometimes the vertices are too many for each polygon, which makes the plotting process extremely slow. Here I am using a tool called DecimatePoly to reduce the total number of vertices 10 times less without making noticeable changes on the shape.

Next we will define the color of the polygon by one of the variables in the attribute table using makesymbolspec:

faceColors = makesymbolspec('Polygon',{'INDEX',[1 lenS],'FaceColor',color_att});
%lenS is the number of polygon
%color_att is a 3-column matrix of RGB triplets with length of LenS

Next is to create a map frame to plot on:

ax = usamap(latlim,lonlim); %plot frame using usamap as template 
axis off; framem on; gridm on; mlabel on; plabel on;
setm(gca,'MLabelLocation',10) %interval for meridians labeling
setm(gca,'PLabelLocation',5)  %interval for parallels labeling

Then plot the polygons

geoshow(ax, huc, 'SymbolSpec', faceColors);

When plotting points on the map, use scatterm

scatterm(lat,lon,s,c,'filled'); % s is vector defines the size of the dots
                                % c is RGB defines the color of the dots

untitled

Python in ArcGIS 1

1. Python editor IDLE will be installed with ArcGIS. Right click .py file can find the IDLE option in the menu.

2. print

strName = "Tian"
print strName
print (strName)
print 10*strName
print 'Tian'
print "Tian"

3. in Python Shell, type into variable name gives the value of the variable. (kinda like matlab)

4. type(strName) gives the type of variable.

5. vars() gives all the variables exist.

6. list contains more variables

lstMonths=['May','Jun','July']
lstMonths[0]   #shows the 1st element of the list "May"
len(lstMonths) #length of the list, like matlab
lstMonths[1:]  #all the elements except the 1st one
lstMonths[:-1] #all the element except the last one
lstYears = range (2000,2014,2) #doesn't contain the last one, 2 is the increment
lstYears.index(2008) #shows the index of this element
lstYears.count(2008) #counts how many 2008 in the list

7. Loops

lstMonths=['May','Jun','July']
for i in lstMonths:
    print i

# or
for i in range(0, len(lstMonths)-1):
    print(lstMonths[i])

8. Modules
First, create a module, save it as “testmodule.py”

#test module

#function do the square
def valueSquared(myVar):
 fltSquared = myVar*myVar
 return (fltSquared)

Second, use this module:

import testmodule
dir (testmodule) #see what function in this module
help (testmodule.valueSquared) #show the comment above this function
testmodule.valueSquared(2.0) #will return 4.0

Another way to import the module:

from testmodule import*
valueSquared(2.0) #then you don't need to put the module name (prefix) all the time

9. OS module
OS module operates files:

import os
os.getcwd() #returns the current working dir
os.chdir('C:/temp') #changes dir
os.listdir('C:/temp') #list files

10. glob module

import glob
glob.glob('C:/temp/*.py') #glob can find the files using wildcard

Delineate basins based on flow directions

1. Download the flow direction file (FDR) from here (Wu, et. al. 2012) then import it into ARCGIS.
Capture

2. Using “Basin” tool to identify the basins and then convert the generated raster to polygons.
Capture

3. Select the desired basin and output the polygon as a single file (Mekong Basin here).
Capture

4. Use “extract by mask” to extract the flow direction for the desired basin.
Capture

5. Export the flow direction as ascii