Friday, December 10, 2010

matlab demo videos

I watched 14 of the demo videos yesterday:

http://www.mathworks.com/products/matlab/demos.html

It took me about 5 years to haphazardly learn most of the things that you could learn in an hour or so watching all of those videos.

Watch the videos.

You learn things like, did you know you can specify something like
mymatrix(10:end);

I've always done
mymatrix(10:length(mymatrix));

Also good:

http://www.mathworks.com/videos/matlab/optimizing-performance.html

I love matlab.

Tuesday, November 30, 2010

plotyy

plotyy - plot with two y-axes, independently control each axis

notes:
(1) specify the x-axis to be exactly the same for both
(2) define the handles for the axes and lines/objects AX(1) is left axis handle, AX(2) is right axis handle

long convoluted example:
    [AX,H1,H2] = plotyy(jda(event),sei(event), ...
        jda(event),dis(event));
    


    set(get(AX(1), 'Ylabel'), 'String', 'seismicity (m/s)', ...
        'FontName', 'Arial', 'FontSize', 12, 'FontWeight', 'bold')
    set(get(AX(2), 'Ylabel'), 'String', 'discharge (m^3/s)', ...
        'FontName', 'Arial', 'FontSize', 12, 'FontWeight', 'bold')
    set(AX(1), 'YLim', [0 .000001], ...
        'XLim', [axismat(1) axismat(4)], ...
        'XTick', axismat, ...
        'YTick', ...
        [0 .0000002 .0000004 .0000006 .0000008 .000001], ...
        'LineWidth',1.5, ...
        'FontName', 'Arial', 'FontWeight', 'bold', 'FontSize', 12) 
    set(AX(2), 'YLim', [0 5000], ...
        'XLim', [axismat(1) axismat(4)], ...
        'XTick', axismat, ...
        'YTick', [0 1000 2000 3000 4000 5000], 'LineWidth',1.5, ...
        'FontName', 'Arial', 'FontWeight', 'bold', 'FontSize', 12) 
    set(H1, 'Marker', 'o', 'LineStyle', 'none', 'MarkerSize', 3);
    set(H2, 'LineWidth', 2);
    axis ([AX(1) AX(2)], 'square');
    xlabel('julian day','FontName', 'Arial', ...
        'FontWeight', 'bold', 'FontSize', 12);

    set(gca, ...
    'LineWidth',1.5, ...
    'FontName', 'Arial', 'FontWeight', 'bold', 'FontSize', 12) 



Monday, November 29, 2010

How to put Excel data into Matlab and retain the 'cell' structure

1. Make your cell structure by using curly brackets:
A={1}
2. Copy the excel data to the clipboard
3. Then open up "A" and right click on the top left cell-- choose 'paste excel data'
4. Save your variable A
5. Use it in Matlab:
B = cell2mat(A(row # to start on : row # to end on , column # );

Monday, October 25, 2010

calculate river slope from coarse DEM

Here is one way to calculate the slope of a stream (and long profile) from a coarse DEM.

(1) If the DEM is huge, clip it to only the area of interest, to speed up the processing. You can do this by creating a polygon outlining the area, adding it to the map, and then using extract by mask (Spatial Analyst, Extraction, Extract by Mask).

(2) It may help to visualize things by making a Hillshade.

(3) Under Spatial Analyst, Hydrology, there are three steps.

(a) Fill

(b) Flow direction

(c) Flow accumulation

(4) Now you want to extract the stream by choosing an upstream point, then selecting all points with a flow accumulation greater than that value, using the conditional to turn the map to 0s and 1s to delineate the stream, and then turn that stream into a vector.

Map algebra, Single Output Map
Conditional

Hydrology
Stream to feature (may be optional and affect choice in Step 5)

(5) Select equally spaced points along this line either by using Hawth’s tools / Geospatial Modeling Environment or this previous post.

(6) Extract Values to Points, Export data

(7) You can open the .dbf in Excel or Gnumeric to inspect the slope.

There’s also this:
Tarboton, D. G., R. L. Bras, and I. Rodriguez-Iturbe. 1991. On the Extraction of Channel Networks from Digital Elevation Data. Hydrological Processes. 5: 81-100.

Monday, October 18, 2010

macro in perl

This is one way to write a macro in perl.

Useful if you want to use perl to manipulate files in the directory, and then run sac commands on them. The sac is run all at once after the if/elsif statements are completed.

macro: A sequence of commands stored for later use in an application

# count number of files in directory

$dir = '.';
@files = <$dir/*>;
$count = @files;

print "numfiles: $count \n";

#######################

# to run sac commands, write to sac macro and run
open(scratch1,'>sacmacro');

# set initial counter
$x = 1;

foreach $fn (<*BHZ*.SAC>){

if ($x==1){

$fnout=$fn . '.merge';
print scratch1 "r $fn \n";
print "read first file: $fn \n";
$x++;

} elsif ($x < $count) {

print scratch1 "merge $fn \n";
print "merge $fn \n";
$x++;

} elsif ($x == $count) {

print scratch1 "merge $fn \n";
print scratch1 "w $fnout \n";
print "merge $fn \n";
}

} # ends the foreach

close(scratch1);

# now run the sac macro
# within the ` `, until EOF, all lines are executed in sac

`sac << EOF
m sacmacro
q
EOF`;

Wednesday, October 13, 2010

Clip a Raster to a polygon - GIS

Make a new polygon with the shape you want.
Toolbox:
Data Management Tools>Raster>Raster Processing>Clip
Your raster in the "input raster"
Your polygon is the Output Extent
(check the "Use Input Features for Clipping Geometry (Optional)" box)

signal processing

spectrogram

h = figure;
Nwin=2048;
Twin=Nwin*dt;

[dum,f,t,p] = spectrogram(y,Nwin,0,[],1/dt,'yaxis');

surf(t,f,10*log10(abs(p)),'EdgeColor','none');

periodogram

making an envelope

envelope

downsample

upsample

diff

Tuesday, October 12, 2010

ftp-ing things

some things that i forget:

lcd: change local directory

prompt: turn off interactive prompting (if you want to mget everything but not be asked every time)

mget: get things and put them in the local directory

if not using the command line i like fileziilla

Tuesday, October 5, 2010

symbolic/soft links

in unix,

ln -s [directory or file]

to point to a directory or file and have it show up in the directory listing but not have the files copied there.

useful if you have files in multiple places (data folders, etc.) but you want them to all show up in one place.

Friday, October 1, 2010

mac vs. pc, matlab version

There are a bunch of reasons why a matlab script written for a pc might not work on mac and vice versa. Here are some examples (some (most) content contributed from loopy):

1. dir function to display directory listing (like ls in unix)

On a pc, you'll get this:

>> dir

.          ..         pathdef.m


. stands for "this directory" and .. stands for "directory above"

As far as I know, mac does not do this, but just starts with the regular file names. A lot of times I want to cycle though files in a directory, so I'll create an array with all of the file names, populating it with dir. On my pc, I start at n=3 to skip . and ..    but on a mac it should start at n=1. Or I should write some better code that just skips . and ..

2. for some reason my mac version never runs out of memory, whereas on the pc some of my scripts always run out of memory. same computer hardware, but the software handles memory differently.

3. folder / \ are different, and sometimes it cares, sometimes it doesn't.

4. mac can't write to excel files - there might be some individually written codes out there to handle this. [note: ugh who would write to .xls?]

5. writing pdfs in pc usually makes them sideways, or 90 degrees off, the mac version.

If there are any others out there please post.

Thursday, September 30, 2010

Create points from a line (ArcGIS)

I need two text files, one with the x coordinates and another with the y coordinates of a stream (to put into matlab). I have a polyline shapefile of the stream in ArcGIS. Don't know if this is the easiest way, but it works :)

1. In ArcCatalog, make a new point file with the same coordinate system as your line file.
2. In ArcGIS, start editing your POINT file, and select the line you want to make poins for (so that it is highlighted green).
3. Under the editor pull down (in the editor tool bar) select to 'divide'.
Choose the second option: 'place points separated by every __ units' (the units are your map units... you can check this under dataframe properties).
4. Stop editing.
5. Open the 'Add XY Coordinates' tool. (ArcToolbox>Data Management Tools>Features)
Choose your point shapefile and say OK.
6. Once it's complete, open the attribute table for your point shapefile.
7. Under attribute table options, choose 'export' and save your file as a .dbf file.
8. Open this file in excel (excel is the default program to open .dbf files anyway).
9. For my matlab application I need two files, one with just x coordinates, and one with just y coordinates, so I made two copies of this file and deleted all but the POINT_X column in one and all but the POINT_Y column in the other. Then for each I deleted the title row, changed it to scientific notation and saved the file as a text(ms dos) file.
10. Wahoo! maybe there is an easier way to do this, but it works :)

some chinese translation tools

OT = off topic

(1) old standby babel fish - type in english or paste chinese characters, get unintelligible response

http://babelfish.yahoo.com/translate_txt

(2) draw the character with your mouse and get the closest matches. okay for characters without a hundred strokes. can't cut and paste the result easily, but you get the pin yin and tone.

http://www.chinese-tools.com/tools/mouse.html

(3) chinese word processor and dictionary. extremely useful if you know pin yin. might only be 30 days free.

http://www.njstar.com/cms/njstar-chinese-word-processor-download

Monday, September 27, 2010

vi stuff

some vi stuff

show line numbers:
:set number

copy or move block of lines or line:

:line,linecotarget
:line,linemtarget
:linemtarget

Go to a line:

type line number then G

error in transfer

If you get this error when applying transfer in sac:  

check_channel; decimation blockette with no associated filter,
skipping to next response now
No transfer function applied! 
ERROR 2118: No transfer function applied.


Then go into the RESP file and make sure that the blocks for each stages are in the following order (cut and paste if necessary):

FIR response
Decimation
Channel Gain

Friday, September 24, 2010

perl

some perl stuff

http://esdynamics.geowissenschaften.uni-tuebingen.de/~esd/tutorials/perl/start.html

http://www.perl.com/pub/2008/04/23/a-beginners-introduction-to-perl-510.html

http://perldoc.perl.org/index-tutorials.html


#!/usr/bin/perl -w

foreach $fn (<*BHZ*.SAC>){
   foreach (`sachdr -l $fn`){
      if (/kstnm/){
          @a=split;
          $sta=$a[2];
          }
      if (/stla/){
          @b=split;
          $stla=$b[2];
          }
      if (/stlo/){
          @c=split;
          $stlo=$c[2];
          }
       }
print "$sta $stla $stlo\n";
};

x windows

Xming for xwindows.

Run Xlaunch for multiple windows.

Then in SSH under tunneling settings, check Enable X11.

Wednesday, September 15, 2010

station lat/lon to 3 decimal places

On the IRIS site the seismic station locations are only given to 2 decimal places lat/lon, which is +/- l km. They don't want to real the precise position of the instruments.

If you need more info, location to 3 decimal places is available in the .sac files.

(from rob)

In the sac header files you can find the latitude and longitude to at least 3 decimal places. You can access this either via sac:
read (file)
lh stla stlo

or a couple matlab options based on which sac reader you're using
with load_sac:
[header, data] = load_sac('file');
header.stla, header.stlo

(In the above two cases stla is short for station latitude and stlo is short for station longitude)

or with readsac
[time, amp, header] = readsac('file');
header(17)
header(18)

in this case header(17) is the latitude and header(18) is longitude.

italian insar site

This link has been stuck in my inbox for a long time, courtesy of JP. Might as well archive it here. We're all about databases.

http://webgis.irea.cnr.it/webgis.html

IREA-CNR InSAR WEB GIS

(wow. acronyms.)

Thursday, August 26, 2010

matlab + geospatial

futzing around on matlab help site and ran into this:

Geospatial Data Import and Access
http://www.mathworks.com/access/helpdesk/help/toolbox/map/ref/f3-12193.html#f3-32594

which includes things like:

usgs24kdem
Read USGS 7.5-minute (30-m or 10-m) Digital Elevation Models

Really?? I guess you need Mapping Toolbox and I'm not sure if that comes standard.

sac remove instrument response

*InstRm: remove instrument signal from seismic signal
[lhsu@thought leslie]$ more InstRm

* read signal
r 2008.005.00.00.16.1150.XQ.ME30.01.BHZ.M.SAC

*possibly transfer first, but remove mean and trend
rmean
rtrend

*taper ends
taper w 0.1

*transfer function to remove instrument noise
transfer from evalresp to vel freqlimit 0.0009 0.001 15 17

*convert to m from nm (sac default)
mul 1e-9

*write out the resultant velocity vector (binary)
w Day005_vel

*fft
fft

*write out the power spectrum (binary)
wsp Day005_sac

Thursday, August 19, 2010

more gis tools

Geospatial Modelling Environment
http://www.spatialecology.com/gme/gmecommands.htm

Yesterday Jon downloaded something from the predecessor of these tools, and the interface and functionality looked really good. You can do a bunch of stuff with points, lines, polygons, that I think is more intuitive than arcGIS.

Also, linked to that page is a whole list of open source GIS. Who knew??

Quantum GIS
GRASS GIS
uDIG GIS
SAGA GIS
OSSIM Image Processing
The R Project
PostGIS
Wikipedia GIS Page

Brain about to explode.

Wednesday, August 18, 2010

mapping software alternatives

Up until last week, I didn't know about any viable alternatives to arcGIS.

But here are a couple (one still under development...)

GeoMapApp

Open Earth Framework

Both are being developed at academic institutions, the first at Lamont-Doherty and the second at UCSD.

I haven't used them very much yet, but I hope to check them out more in the near future. 

Thursday, July 29, 2010

sea monkey

New html editor. Open source, free.

http://www.seamonkey-project.org/

Completing my very perplexing office suite of Google docs, TeXnicCenter, gnumeric, and gimp.

trace of a river profile

A somewhat less lazy explanation of how to extract the trace of a river profile. [lazy explanation]

(1) Download the relevant file from the National Hydrography Dataset. (I recommend using the ftp site, which unfortunately only works on Internet Explorer.)

(2) Open the hydrology layer in arc.

(3) Determine the source latitude and longitude of your river of interest. You may be able to find this under the Geographic Names Information System entry for your river.

(4) Go to that lat/lon in the arc file and locate the source.

(5) Select the feature using Navigation [from this tutorial]

a. Click on Analysis in the Utility Network Analyst toolbar, click Options, and then Results.

b. Click on the Selection radio button, Apply, then OK.

c. Perform navigation. The navigation path will select and display NHDFFlowlines.

i. Click on View, Toolbars, and then check the Utility Network Analyst box.
ii. HYDRO_NET will appear automatically in the Network box.
iii. click the Add Edge Flag tool. The flag marks the start of the navigation path.
iv. click cursor at appropriate point on flowline (the source).
v. Select Trace Task, Trace Downstream
vi. Click on the Solve icon - it should highlight the river downstream of the flag

(6) Save the new shapefile. Right-click on NHDFlowline in the Table of Contents, scroll to Data, click on Export Data. Export selected entries, name the new shapefile.

Monday, July 26, 2010

Quaternary faults, reproject

Quaternary Faults in California, maintained by California Geological Survey:
http://www.consrv.ca.gov/cgs/information/publications/Pages/QuaternaryFaults_ver2.aspx

(1) Download .e00 files
(2) Import from Interchange File

(3) Data Management Tools -> Projections and Transformations -> Feature -> Project
Fill in:
Input Dataset
Input Coordinate System (automatic)
Output Dataset (put what you want, specify folder you want)
Output Coordinate System (for me: NAD_1983_UTM_Zone_10N)
Geographic Transformation (for me: (guessed) NAD_1927_To_NAD_1983_NADCON)

Aside:
wth is NADCON
http://www.ngs.noaa.gov/TOOLS/Nadcon/Nadcon.shtml

dealing with .e00 files v.2

Converting .e00 files. An easier way than here:

In ArcCATALOG

View -> Toolbars -> ArcView 8x Tools -> Conversion Tools
Import from Interchange File

data frames, neatlines, graticules

In ArcMap sometimes you want to add lat/lon ticks around your map, add a scale bar, etc.

This stuff is all filed under Page Layout and Map Composition. You see most of this stuff under the Layout View, not the Data View. (Buttons in lower left corner of map window toggle these.)

data frame: most fundamental element in a map document and in ArcMap user interface.

neatline: border line commonly drawn around the extent of a map to enclose the map, keeping all o the information pertaining to that map in one "neat" box

graticule: the grid of intersecting lines, esp. of latitude and longitude on which a map is drawn. vs. gridicule.

Accessing the Grids and Graticules wizard:
  1. Click the View menu and click Data Frame Properties.
  2. Click the Grids tab.
  3. Click the New Grid button.
Adding grids and graticules (reference systems)

Working with data frames
Adding data frames
Customizing data frames (you can rotate them)
How to do inset maps (extent rectangles)

Thursday, July 22, 2010

gnumeric

I am now using gnumeric as my spreadsheet software on my work computer.

http://projects.gnome.org/gnumeric/

You can download a Microsoft Windows build if you want. [.exe]

Not too many variations on a spreadsheet program. It's not that I actually reject Excel and Open Office Calc, it's just that they won't successfully install on my !#$^(!&#)^ computer.

spatial hydrology data

I've been trying to download river spatial data (for arc) for a while now, and finally have come up with some sort of resolution. But the problem seems to be information overload. I still haven't figured out if there's a way to, say, query "Truckee River" and have the trace of only that river returned.

Instead, here are some convoluted ways to get data of Quads or defined map areas.

(1) USGS maintains the National Hydrography Dataset.

(2) The Data Availability page describes ways to obtain the data.

(3) It appears that USGS is in the middle of transitioning to a new map viewer/data downloader interface.

New one: http://viewer.nationalmap.gov/viewer/

It is easy to download hydrology data by quadrangle, if you click on the map, it tells you the quadrangle and the available data. (Using the download data tool.)

(4) One way to figure out all of the quadrangles covered by a particular river is to use the Geographic Names Information System.

(5) The page for the river will list each quadrangle, which you can then locate in the Map Viewer.

(6) You send in a request for a data and it may take several hours to be returned. (Still waiting for some I ordered a couple hours ago.) It is not very efficient.

(7) Another way is to download pre-defined subregions. These datasets are ~250 MB zipped and ~750 MB unzipped. It is information overkill, but you can also download it right away. The old map viewer seems to be better in identifying the 4-digit identifier (first 4 of the 8 digits).

FTP for predefined regions:
ftp://nhdftp.usgs.gov/SubRegions/

Only works in Internet Explorer! Not Firefox or Chrome! Also has a README in a .doc format! gah.

I'm not sure how much closer this gets me to the trace of a single river, but possibly this information will be useful to someone somewhere.

Here's their tutorial.
http://nhd.usgs.gov/tutorials.html

Note:

The last section of the tutorial, on Navigating, actually lays out how to save the trace of a river. First you need to click on Analysis in the Utility Network Analyst toolbar, and click Options, Results, and check Selection.

Then "perform navigation," set an Edge Flag at the source of the river (this is located in the Geographic Names Information System entry, or at least the correct quadrangle). Then use the Utility Network Analyst to navigate to all points downstream.

You will be able to save it to the map and export the data as a shapefile.

Read the last page of the tutorial linked above.

Wednesday, July 14, 2010

Export Table to ASCII

Under the Spatial Statistics Tools, under Utilities, there is a tool called Export Feature Attribute to ASCII.

This tool will write out a .dbf Attribute Table to ASCII, thereby negating my previous bitching about it. You can choose which columns to write out, and you can choose the delimiter.

Great, Way to make me feel dumb, arc. I can still complain that this tool is hidden where it is - why not just have a button for Export .dbf, or Save As? Why?

Oh well, now we know.

Tuesday, July 13, 2010

refreshing and labeling

In ArcCatalog, if you've just created a new layer in ArcMap, it may not appear until you refresh the list. To refresh the list, use F5.

Labels

You can label features (ugh, an arc word...) with their attributes (ugh, another arc word). There is a toolbar for it.

You might need to be in an open Edit Session for this to work.

View -> Toolbars -> Labeling

The first icon is the Label Manager.

You can label with more than one attribute field. Click on the Expression button under Text String. You can Append more than one field name in the label. (e.g. you can label a fluvial terrace with both its unique id number and also the terrace grouping Qt1, Qt2, etc.)

Edit the size and color and etc. to your liking.

The end.

Text in ArcGIS is nontrivial (labels vs. annotation vs. graphic text), but labeling is relatively straight-forward.

Overview of working with text

Friday, July 9, 2010

linear referencing and hatch marks

Linear referencing: point location along the line as an alternative to expressing the location using an xy coordinate. Useful, for example, for measuring distance along a stream.

An Overview of Linear Referencing

Here are steps outlining a way to put hatch marks (with distance labels) on a polyline (e.g. following a stream).

1. When you create the polyline shapefile in arc catalog, you need to check the box for "Coordinates will contain M values. Used to store route data."

2. Make the polyline (start edit session, click click click, double click to end, stop edit session)

3. In the arc toolbox -> linear referencing tools -> create routes  [routes?]

a. Input line features: the polyline layer
b. Route identifier field: (you can use id or if you created another field like "route." If there's only one line, that's ok.)
c.  Output route feature class (it fills in a default)
d.  Coordinate Priority: where is the start of the polyline? In my case, lower_left.
e.  click OK - this will create another later eg. channel_CreateRoutes

4. I'm still fuzzy on why you need an entirely new layer, but anyway... it is on this new layer that you can specify hatch marks.

This page nicely outlines the 19 (!!!) steps necessary to display hatch marks on that layer:

Displaying hatches

  1. Right-click the layer in the table of contents that you want to hatch and click Properties.
  2. Click the Hatches tab.
  3. Check Hatch features in this layer.
  4. Type an appropriate hatch interval.
  5. Click Hatch Def(1).
  6. The Hatch Definition view becomes active.
  7. Type an appropriate line length.
  8. Right-click the hatch class and click Add Hatch Definition.
  9. The Hatch Definition view becomes available.
  10. Type an appropriate multiple of the hatch interval.
  11. Type an appropriate line length.
  12. Check Label these hatches.
  13. Right-click the hatch class and click Add End Hatch Definition.
  14. The End Hatch Definition view becomes active.
  15. Type an appropriate end hatch tolerance if necessary.
  16. Type an appropriate line length.
  17. Check Label these hatches.
  18. Click OK.
  19. Poke your eye out with a pen.

5. Peripherally related: It may also be useful to you to divide a line into segments: Splitting line features

6. Next up: putting something useful along the route: Event Tables

in the beginning...

Sometimes you have to start from the beginning, there can be useful information in the beginning.

LaTeX 
"haha, if you're not smart enough to figure this out, too bad for you."

http://en.wikibooks.org/wiki/LaTeX/Introduction
http://en.wikibooks.org/wiki/LaTeX/Importing_Graphics 

Arc
"haha, we don't give a fck about user experience, too bad for you."

http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=An_overview_of_the_Geodatabase

Wednesday, July 7, 2010

jabref and TeXnicCenter and BibTeX

BibTeX, I hate you for making me feel stupid. But I love you because you are going to make my life a lot easier. right? Right?! RIGHT??!!

1. JabRef databases are saved in .bib format, which are just .txt files and can be opened with WordPad or vi or whatever. They might have a bunch of extra crap in them like full abstracts, the notes you typed into JabRef, whatever. Export the entries for your paper into a separate .bib file: blah.bib.

2. Using the BibTex keys in JabRef, cite the stuff in your .tex file (in TeXnicCenter on a dumb PC if you are like me)

\citet and \citep are for textual and parenthetical citations
Reference sheet for natbib usage

\citep{Dietrich2003, Gauer2004, Stock2006}

3. You have downloaded all of the necessary style files or whatever, for example from the AGU website, and put the agu08.bst file into your working directory. In your working directory, you also have your .tex file and the blah.bib file. In the .tex file you have referenced the .bst and .bib you are using.

\bibliographystyle{agu08}
\bibliography{blah}

4. Build your document in TeXnicCenter. In the profiles for building files (e.g. for LaTeX -> PDF) there is a box to check for running BibTeX (mine was already checked by default). This first Build creates a .bbl file in your working directory.

5. Cut and paste the contents of the .bbl file into your .tex file AND comment out the two lines above in #3.

6. Build the document again, it will now use the \bibitems in the .tex file and everything will work great and you will feel great too. (You might need to build more than two times for all cross references.)

Have fun.

[9/30/2010 ps. It would be very helpful to know if this explanation is useful to you, or if you need/want more detail, or if it worked for you or not. There seems to be the most demand for this information and I really hate how there is not much info out there. I will monitor comments.]

dealing with .e00 files

.e00 files are from an older version of arc. a simpler time. but now it is complicated to open them.

update 7/26/2010:
Easier way:

In ArcCATALOG

View -> Toolbars -> ArcView 8x Tools -> Conversion Tools
Import from Interchange File

**** old stuff below ***

I ended up importing some .e00 files that were in a .tar.Z file (1996 USGS Open report) in a very roundabout manner, and later found some other more straightforward procedures (but can't figure out how to get those to work with my version of Arc). I'll just post them all here.

1. Official explanation on arc site:  Importing an ArcInfo interchange file

2. Mac person's explanation:
there is a arctool box feature for importing arc/info exchange (.e00) files. Toobox-> Coverage -> Conversion -> To Coverage -> Import exchange files. I was able to unzip the tar file using the Mac Archive Utility.

3. My convoluted method on dumb PC (involving much sighing and swearing):

1. download tar.Z to PC
2. ftp to unix account
3. uncompress and untar in unix

uncompress *.Z
tar -xvf *.tar (I may have left off one of the flags to get it to work)

4. ftp back to PC
5. use Import71.exe (wtf, really) to import the .e00 files to the newer version (see below)

my [source]

Importing .e00 files using Import71.exe, ArcView 9:

a. Before you start, note that this import tool does not deal well with folder paths that have space in their names (e.g., like "My Documents"). Therefore, prior to using this tool, you should create a new folder for that file in your C: drive. E.g., if I was importing a TxDot county file for Williamson County called urban246.e00, I should create a folder called williamson in the c:\temp\ folder, and copy my urban246.e00 file to that folder.

b. Locate the Import71.exe file - it should be in the C:\arcgis\bin\ folder (or use the Windows search tool to find it)

c. Run the Import71.exe file by double-clicking on it - a dialog box should appear

d. Browse to the location of your .e00 file and select that file and its directory path as the input file (e.g., c:\temp\barbara\williamson\urban246.e00)

e. For your output file, simply specify the name you want for the newly imported file (e.g., txdot) - no directory path is needed, the file will import to the same location as the .e00 file

f. Be patient - the import process may take 10 minutes or so. You will know when it is done because a small box saying Import Complete will appear over the Import71 dialog box.

g. After you see the Import Complete notice, start ArcView again and check that the coverage displays properly. A "coverage" is an older GIS data format.

contours

Simple:

Toolbox -> Spatial Analyst Tools -> Surface -> Contour

Specify input raster, output name, contour interval.

Creating contours in spatial analyst

x, y, z table for a bunch of points

This is how I got an xyz table of a bunch of points. (There is probably a better way out there.)

1. Using arc catalog, create a shapefile of points. Then add columns to the attribute table, x, y
2. In arcmap, start an edit session and click away on the points.
3. End and save edit session, I think.
4. Right click on shapefile in Table of Contents column and open Attribute Table. All of the points should be there.
5. Right click on the x column
6. Click Calculate Geometry
7. Choose X coordinate of Point, Use coordinate system of the data source (UTM meters)
8. Repeat with x column and Y coordinate of the Point
9. Use Toolbox -> Spatial analyst Tools -> Extraction -> Extract Values From Points
Input point features (point shapefile)
Input raster (DEM)
Output point features (choose your name)
Click OK
10. This should create a new column in the table. Now you can save the table and export the .dbf if you wish, and do whatever you want with it.

Overview of Extraction tools

aspect-slope map plus other stuff

While hunting around for some guidance in producing a slope map (a matter of a few mouse clicks), I came across this page [aspect-slope map], which is useful because it gave an example of how to implement the following things:

Styles (predefined colors, symbols, map elements, etc.)
Slope map (with link to how it works)
Aspect map
Reclassify: redefine the default number and width of display bins
Map Algebra

Thursday, June 24, 2010

lazy about tex equations

It's only been a 15 year hiatus but i am now seriously working on a .tex file.

Although I have no problem with brute-force typing in equations, the lazy part of me wishes there were something that would take the MathType output and generate the equation code. Just for efficiency's sake, you know? There are a lot of equations.

I found this thing, which helps. You have to click in the symbols, but for a beginner like me, it's a good way to learn the syntax.

http://www.codecogs.com/latex/eqneditor.php

.dbf files

ArcGIS output may be saved as a .dbf file, with no immediately apparent way to save as a .txt, or even to copy and paste the cells to a simple spreadsheet program (why? why?).

I had no idea what to do about this so I just downloaded a small program that converts the .dbf to .csv.

http://www.fileguru.com/DBF-To-CSV/download

I did this a couple days ago so now I have no recollection of how it worked, I just remember that it worked, and I copied the data into a google docs spreadsheet. Open Office won't even install on my computer, how lame is that?

If anyone else knows a more intelligent way to do this, please share.

Somebody else talks about .dbf and .xls:
http://gisatvassar.blogspot.com/2008/02/excel-dbf-files-and-arcgis-92.html

Note, 7/14/2010, This is the way to do it:
http://matlabor.blogspot.com/2010/07/export-table-to-ascii.html

Monday, June 21, 2010

Elevation statistics of traced polygons

I wanted to find the mean, min, max, stddev elevation of a bunch of traced terraces, but this can be done with any input files of a zone dataset (i.e. shapefile of polygons) and a value raster (i.e. raster of elevation values).

1. Make the zone dataset, in my case, make a shapefile of polygons [make polygon layer]

2. Add field(s) to the Attribute Table to define the zones. I have a separate ID for each polygon and also the field "Name" for the groupings Qt1, Qt2, etc. To add fields, open the shapefile in ArcCatalog and in the Fields tab, just click in the table to add the Field Name and Data Type (must be integer or text).

Then edit the field to define your zones. First add the zone dataset to the map.
Then, Editor toolbar -> Start edit session -> right-click on your shapefile layer in the leftbar -> Open Attribute Table -> click in the table and fill in.

Stop edit session and Save changes.

3. Zonal analysis by attribute
Spatial Analyst toolbar -> Zonal Statistics
Fill in stuff in the pop up window:
Zone dataset: shapefile of polygons
Zone field: ID or Name or whatever your defining field is
Value raster: elevation dataset
Check or uncheck Ignore NoData and Join output table to zone layer (I like to join)
Chart statistic: I checked Mean but got min, max, range, mean, std, sum in the output table
Specify output table

4. Click OK, the output table pops up, and you're Done.

Friday, June 18, 2010

make polygons

To make some @#$*%$@#% polygons in arcgis (e.g. trace fluvial terrace deposits):

1. Open Arc Catalog
2. In the Location bar, navigate to the folder where you want to put the polygons (shapefile)
3. File -> New -> Shapefile
4. Name it and choose Feature Type -> Polygon
5. If you want, edit the Coordinate System
Select -> Projected Coordinate Systems ->
(UTM NAD 1983 10 N for local stuff)

In ArcCatalog you can right-click a shapefile and edit Properties - e.g add Fields and specify the DataType (text, long integer, double, etc.). But I think the shapefile can't be open in ArcMap if you want to edit it.

6. Add the (empty) layer to the ArcMap window by dragging or using the AddData button
7. Start an edit session (Edit Toolbar)
8. Task dropdown menu: Create New Feature
9. Target dropdown menu: choose your shapefile - if the target drop down is empty, it might be a bug, restart ArcMap and ArcCatalog

10. Use the Sketch button (pencil icon), click and double-click to make vertices and close polygon
11. Stop edit session, Save

12. *Get mean elevation of each polygon from an elevation raster (TBD)

obtain and view seismic data from IRIS

This is deteriorating from a pedagogical tool to a storage place for cryptic notes. But, better than nothing.

One way to get some seismic data from IRIS and look at it

1. Use the GMAP (google map) interface on the IRIS site
http://www.iris.washington.edu/gmap/

The instructions tell you what to enter in the URL
lat/lon example:
http://www.iris.edu/gmap/?minlat=46&maxlat=49&minlon=-125&maxlon=-117

2. A map will load in your browser, Click on a station, click on the "more information" link

3. Look for BHZ data (Broadband H-something z-orientation )

a. cursory look
   4. Click on "Explore with QUACK"
   5. Click on "NEIC PDF PLOT"
   6. Click on "by Year" and look at monthly variation

b. download the actual data
   7. Click on "Make a batch request for data (breq_fast)
   8. Network, Station, and Channel should already be filled in for you, fill in other required (starred) data - Start and End time (GMT - HHMMSS), and contact info
   9. In some amount of minutes, you will be emailed with a link to the data, download it to linux desktop

10. In the directory with the .seed file, unpack the .seed data with rdseed (make sure your .cshrc and .login are set up for this)

Use all defaults except:
Options: d for dump
R for raw

you might need to edit Start/End/Sample Buffer Length depending on the length of the data and what you want

You will get a bunch of .sac files

12. Use pql to view data

13. Use sac to view data

read: r *.SAC
plot: p1
bandpass: bp
ylim
etc.

14. Use matlab (yay. hopefully there will be some posts about this at some point)

Wednesday, June 16, 2010

hillshade from seamless data

It's ArcGIS time. First, let's reflect on how Arc is the anti-google in terms of user-friendliness. It's a constant source of exasperation and I am here to document my exasperation and hopefully avoid some in the future.

Starting from the beginning, here's an example if you want to download high resolution topographic data (1/9 arc second, where 1 arc second at the equator is ~31 m)  from the USGS Seamless Server, reproject it to UTM, and then make a hillshade.

1. Download stuff from seamless - http://seamless.usgs.gov/

Click on the "view and download united states data" image [link].

a. On the right hand side of the screen:
In the Display Tab, under Layer Extents, specify NED 1/9 arc second (National Elevation Dataset). This shows where that data exists.
In the Download Tab, click 1/9 arc sec, assuming that's what you want (default is 1 arc).

b. On the left hand side of the screen, under the Downloads section, define your area:
The way I like to do it is with the 2nd button, Define download area by coordinates, because then you know exactly what the boundaries are, if you need to mosaic it later or whatever. I like to get my lat/lon coordinates from Google Earth.

c. Fill in the window that pops up and then click "Add Area", another window will pop up telling you your request is in a queue and then the file will be downloaded to your default location. If the area you specify is large (I don't know how large), it will be divided into separate files.

Note: Sometimes one or more of the files will stall during downloading, and since the naming convention is incomprehensible (an 8-digit string, seemingly random), it is hard to figure out what part you are missing. Oh well.

2. Unzip all files in the directory you want

3. If there are multiple files, make a mosaic from the separate files

In ArcMap,

a. Make sure the Spatial Analyst extension turned on (Tools -> Extensions -> Spatial Analyst)
(!$#%!#$)%!#$%& that this is not on by default)

b. Click on red tool box (or bring up the toolbar with View ->Toolbars ->Tools)

Click Data management Tools -> Raster -> Raster Dataset -> Mosaic to New Raster

if you don't specify an extension, it will default to GRID

4. Convert from lat/lon to UTM

In Data Management Tools,
Projections and Transformations -> Raster -> Project Raster -> UTM

Input: make sure to erase the default output name .img extension, and make under 13 characters (bug in GIS)

usual convention for new projection: name_prj

Select -> Projected Coordinate Systems

NAD 1983 Zone 10N (for local CA stuff)

Scroll down to Resampling technique-> change to CUBIC (the default method tends to create striping)

5. Restart ARC because it is possible that it thinks it is still in lat/lon

6. Make hillshade

In tools:
Spatial Analyst Tools -> Surface -> Hillshade

You should have a nice hillshade of the topo data.

Tuesday, June 15, 2010

LaTeX on windows

Matlabor is back, but not only about matlab. It's a place to put things that I would otherwise forget.

Today we're going to document how to do LaTeX on windows. Why? Let's just assume the following
(1) I had ubuntu working but when I tried to upgrade it yesterday, I ended up three hours later very demoralized with no more ubuntu on my computer (the fault of incomplete wubi uninstall, and my incomplete understanding of computers)
(2) I've always wanted to be cool enough for LaTeX
(3) I'm just not cool enough for a mac

Here are some steps
(following http://doams.wordpress.com/2009/10/11/how-to-setup-latex-on-windows-7/)

1. install MikTeX
http://miktex.org/2.8/setup

2. install TeXnicCenter
http://www.texniccenter.org

When you start TeXnicCenter for the first time you will be asked where your MiKTeX executables are located (in my case it is in …\program files\MiKTeX2.8\miktex\bin).

Then you will be asked what program you want t use for displaying Post Script files. Skip that part.


3. download journal-specific templates
e.g.
http://www.agu.org/pubs/authors/manuscript_tools/journals/latex/index.shtml

4. add the new class and style files to the path

I got errors if I didn't exactly follow their directory naming conventions:

Make this directory and put the .cls and .sty files in here:
C:\Local TeX Files\tex\latex\misc

Add  C:\Local TeX Files
to MikTeX, following the directions here:
http://docs.miktex.org/manual/localadditions.html

5. Open TeXnicCenter, edit template, and ctrl-shift-F5
I like LaTeX -> PDF

Getting the bibliography to work will be for another day, but I'm hoping that everything will work nicely with JabRef.

Wednesday, March 3, 2010

profile & tic & toc

it seems like i have exhausted my knowledge of matlab functions. i am still waiting for guest blogs for patch and ode45. i won't hold my breath. or perhaps i will.

here's something pc told me about once, i haven't used it much because my programs aren't that intensive. (not intensive for the computer, but very intensive for me to produce)

profile

The profile function helps you debug and optimize M-files by tracking their execution time. For each M-function, M-subfunction, or MEX-function in the file, profile records information about execution time, number of calls, parent functions, child functions, code line hit count, and code line execution time. Some people use profile simply to see the child functions; see also depfun for that purpose. To open the Profiler graphical user interface, use the profile viewer syntax. By default, Profiler time is CPU time. The total time reported by the Profiler is not the same as the time reported using the tic and toc functions or the time you would observe using a stopwatch.

Hm, reminds me of that horrible song by Ke$ha, TiK ToK. horrible horrible song.

anyway. remember to debug and optimize!

Wednesday, February 24, 2010

textread and ignoring columns

okay, say you have a .csv file (comma separated variable, excel likes those) and for come reason you don't want to use csvread. i found textread to be more flexible.

roka collected some data that has 22 headerlines, and then 4 columns of data. i only care about the first column of data. so this is what i used to read in the file:

[kilograms] = textread(namestr,'%f %*f %*f %*f', 'delimiter', ',', 'headerlines',22);

namestr is set to the name of the file, e.g. '090512.csv'
%f indicates that the data are floating point numbers
the *s indicate that i am skipping the last three columns of floating point data
there are several options for the delimiter, in my case it is commas ','
there are 22 headerlines before the data starts. 

more data i/o, woo.

anyway, i just went to make the links to textread and matlab says that it recommends textscan over textread. oh well. 

Wednesday, February 17, 2010

break your plot. breakplot

there are several options on file exchange for putting a break in your plot axis if you have a region of uninteresting data.

for example if you have 12 subplots and 11 of them have the same y-bounds, but one of them has a very high value, you can put a break in that plot and standardize the rest of the axes bounds. got it?

breakplot

breakxaxis (you can modify it so it's the y-axis)

i haven't actually used these yet, so that's all i have to say about that. but i hope to in the near future.

Friday, February 12, 2010

fread: read binary data

probably the moment that i felt like the biggest matlab stud (although there are so many i cannot even count, haha JOKE), is when i finally decoded the binary data from the laser scanner my flume. the folks at SAFL had written the program to store the voluminous data in binary format, and ever since dec. 2006 i knew that one day i would have to figure out how to read it.

fread reads data from a binary file.

numsamples = fread(fid, 1, 'long','ieee-le');

first, you probably need to use fseek to move around in the file (to skip header info, etc.):

status = fseek(fid, 3000, 'bof');

i love the words bof, cof, and eof (beginning, current position, and end of file)

anyway, then you can read all the data and write it out to a matlab matrix

for i = 1:numsamples
    mirrorrev(i) = fread(fid, 1, 'double','ieee-le');
    drumrev(i) = fread(fid, 1, 'double','ieee-le');
    distance(i) = fread(fid, 1, 'single','ieee-le');
    laserstatus(i) = fread(fid, 1, 'uint16','ieee-le');
    ambient(i) = fread(fid, 1, 'int8','ieee-le');
    amplitude(i) = fread(fid, 1, 'int8','ieee-le');
end

data I/O is very important!

Tuesday, February 9, 2010

suplabel

this is one that i downloaded from matlab central file exchange: suplabel, places one x or y axis label on a bunch of subplots. Because sometimes redundancy is bad.

http://www.mathworks.com/matlabcentral/fileexchange/7772-suplabel

I changed the default supAxes to [.12 .12 .85 .85 ] so they wouldn't be chopped off. Download it to your path and off you go.

Places text as a title, xlabel, or ylabel on a group of subplots. Returns a handle to the label and a handle to the axis.

[ax,h]=suplabel(text,whichLabel,supAxes)

returns handles to both the axis and the label.

ax=suplabel(text,whichLabel,supAxes)

returns a handle to the axis only. suplabel(text) with one input argument assumes whichLabel='x'

whichLabel is any of 'x', 'y', or 't', specifying whether the text is to be the xlable, ylabel, or title respectively.

supAxes is an optional argument specifying the Position of the "super" axes surrounding the subplots. supAxes defaults to [.075 .075 .85 .85] specify supAxes if labels get chopped or overlay subplots

EXAMPLE:
subplot(2,2,1);ylabel('ylabel1');title('title1')
subplot(2,2,2);ylabel('ylabel2');title('title2')
subplot(2,2,3);ylabel('ylabel3');xlabel('xlabel3')
subplot(2,2,4);ylabel('ylabel4');xlabel('xlabel4')
[ax,h1]=suplabel('super X label', 'x');
[ax,h2]=suplabel('super Y label','y');
[ax,h3]=suplabel('super Title' ,'t');
set(h3,'FontSize',30)



the image doesn't match the example code. so kill me.

Monday, February 8, 2010

codecs and videos

you can make videos with matlab, usually i use avifile and addframe (usually in a loop) to make animations of my plots.

matlab will automatically use some kind of compression and this is where you might run into a problem where you aren't able to actually play the .avi on your computer. you might not have the correct codec (encoder/decoder).

These are the options given for the 'compression' parameter of avifile:
'Indeo3'
'Indeo5'
'Cinepak'
'MSVC'
'RLE'
'None'

Most of those are pretty old.

Also, if you want to enrage an open-source person you can start talking about proprietary codecs. I know very little about this stuff, but here's a list of open source codecs.

http://en.wikipedia.org/wiki/List_of_open_source_codecs

What I know is:

1. MATLAB defaults to Indeo5 on Windows, but this isn't part of the base Windows install now, so I had problems playing my own MATLAB generated .avi movies

2. 'Cinepak' seemed to work, but it loses some sharpness, especially in text

3. 'None' is fine but the files are much larger (duh)

4. Picasa will upload the default Indeo videos and play them, but degrade the quality, and mess up the timing

5. Another option is to write out .avis using 'None' compression, then using a third party software to compress.

6. Also some people have written stuff on the matlab file exchange such as mpgwrite, i haven't tried it yet though.

Man, MATLAB is so cool.

Saturday, February 6, 2010

making nice graphs

here's a much more comprehensive link on making nice graphs

http://blogs.mathworks.com/loren/2007/12/11/making-pretty-graphs/

covered (in well-written code):
  • line properties (both functional and aesthetic)
  • legends and labels
  • font and axes properties

Wednesday, February 3, 2010

serial port I/O and pause

i employed the help of yesdogs' friend to write this .m file to collect data from an instrument connected to my laptop's serial port:
http://seismo.berkeley.edu/~lhsu/help/get9830.m

instead of prostituting oneself for matlab code (figuratively, not literally), one could start here:

Serial Port I/O introduction
http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_external/f105659.html

okay i'll paste the .m below. there's a lot to talk about. i don't understand it all. i think i had a version 2 of this that collected data for a set time, but i can't find it.

i think my favorite part is pause(0.25)
pause halts execution temporarily. as written above, it halts things for 0.25 sec. so data is collected at 4 Hz.

note: i have not yet used this data, but i wrote an abstract for a conference in 2011 that requires it.


% get9830.m 
%
% An example of how to get data from an instrument connected to a serial port.
% This was written specifically for an Interface 9830 Digital Indicator
% (connected to a load cell). Therefore some variables - espectially in the
% s = serial(....) line - may need to be adjusted for different systems.
% 
% pause(0.25) sets the frequency of data logging, units of seconds so 0.25
% is 4 Hz.
%
% finalized on 16 june 2008
% written by terry
%
% open serial port on COM1
s=serial('COM1','BaudRate',9600,'DataBits',8,'Parity','none','StopBits',1,'Terminator','CR');

fopen(s);

s.Status

disp('get9830.m');
disp('ctrl-C to end');
input('press enter to begin');
i = 1;

while(i)
    % send transmit on
    %fprintf(s,'data');
    fwrite(s,17,'uint8');  %send an XON
    fwrite(s,13,'uint8');  %send a  CR


    % get first 8 bytes
    %data = fread(s, s.BytesAvailable, 'uint8');
    if(s.BytesAvailable >= 8)
        data = fread(s, 8, 'uint8');
    
        %convert data into decimal
        reading = bitshift(data(3),24) + bitshift(data(4),16) + bitshift(data(5),8) + data(6);  %combine data bytes
        

        %handle negatives
        if(data(3)>127)
            reading = bitcmp(reading,32); % complement the bits
            reading = -1*(reading+1);
        end;

        %place decimal point
        reading = reading*10^-(5-data(7));
        
        %save in array
        output(i) = reading;
        pause(0.25); % I think the indicator makes 4 measurements per second
        i=i+1;
        
        % press ctrl-c to end
        
    end;
end;

fclose(s);
delete(s);
clear s;

Sunday, January 31, 2010

switch, otherwise, disp

more things i should've known and/or always forget:

switch is supposed to be better than a list of if, elseifs. i can think of some cases where i should've used this.

a very clear and not flashy 3 minute video explaining this:
http://blogs.mathworks.com/videos/2008/01/02/matlab-basics-switch-case-vs-if-elseif/

otherwise is the default value of switch (also explained in video)

the video shows the use of disp, which is a very simple command to display text that i have forgotten in the past and spent a long time trying to remember. it is particularly good if you want to print things out to the screen to know how your program is running (e.g. error messages, transient values, loop number, etc.).

disp('SIGH');

Saturday, January 30, 2010

text as TeX or un-TeX

continuing the tradition of writing about stuff that i haven't really implemented in the past...

text properties

important: the default is to interpret strings as TEX instructions, and that is why the underscores in my titles always mess it up by creating subscripts. i need to set the text interpreter to none. then i need to be cool and start having some LaTeX equations in my other plots.

from the Text Properties page:

Interpreter
latex | {tex} | none

Interpret TEX instructions. This property controls whether MATLAB interprets certain characters in the String property as TEX instructions (default) or displays all characters literally. The options are:

latex — Supports a basic subset of the LATEX markup language.

tex — Supports a subset of plain TEX markup language. See the String property for a list of supported TEX instructions.

none — Displays literal characters.

Latex Interpreter

To enable the LaTEX interpreter for text objects, set the Interpreter property to latex. For example, the following statement displays an equation in a figure at the point [.5 .5], and enlarges the font to 16 points.

text('Interpreter','latex',...
'String','$$\int_0^x\!\int_y dF(u,v)$$',...
'Position',[.5 .5],...
'FontSize',16)



Wednesday, January 27, 2010

setting defaults: startup.m

this might have been useful to learn 7 years ago: make a startup.m file and put it in your path, you can set all sorts of defaults like line widths, colors, etc.

http://www.mathworks.com/support/tech-notes/1200/1205.html

see:
Section 5: Setting default property values

If I ever get around to it I'll set font to Arial, size to 14, LineWidth to 2. Or something like that.

Tuesday, January 26, 2010

octave: matlab's cousin

since we all like free things, and open source, here is the link to matlab's cousin, octave:

http://www.gnu.org/software/octave/

i've never used it, but i'd feel like a stud if i knew how.

GNU Octave is a high-level language, primarily intended for numerical computations. It provides a convenient command line interface for solving linear and nonlinear problems numerically, and for performing other numerical experiments using a language that is mostly compatible with Matlab. It may also be used as a batch-oriented language.

...There were still some people who said that we should just be using Fortran instead, because it is the computer language of engineering [if you were born in the 1800s], but every time we had tried that, the students spent far too much time trying to figure out why their Fortran code failed and not enough time learning about chemical engineering. We believed that with an interactive environment like Octave, most students would be able to pick up the basics quickly, and begin using it confidently in just a few hours.

***

Somebody needs to kill FORTRAN, okay? I guess if I actually knew that language, I could help the world be free of it. But I don't, I just like to make disparaging remarks about it. In ~10 years I'm going to be one of those old folks who still tries to hold on to the technology/software/languarges that I learned in my youth. I already feel sorry for myself now.

Monday, January 25, 2010

anti-illustrator

adobe illustrator has its purposes, but let those purposes NOT be to tweak matlab output so that it looks more aesthetically pleasing.  matlab output edited by illustrator one time is tolerable, but if you have to do it two, three, four times before you get your final plot, you're likely to poke your eyes out with a blunt object.

don't be a fool!

not that i've mastered these commands yet, but setting font, font size, line width, and tick properties goes a LONG way. i haven't used illustrator to "fix" a plot since i've started using stuff like this:

xlb=xlabel('strength (kPa)','FontName','arial');
ylb=ylabel('erosion','FontName','arial');

set([xlb ylb],'FontSize',14,'FontName','arial','LineWidth',1.5 )

axis square;

set(gca,'FontSize',14,'FontName','arial', 'LineWidth',1.5, 'Ticklength', [0.025 0.025])

gca is important.

Friday, January 22, 2010

axis of evil

i love me some square axes:

axis square;

i am also always setting axes bounds:

axis([xmin xmax ymin ymax]);

but there are tons of other axis options i have not checked out fully.

axis 'auto yz';

axis auto
axis manual
axis tight
axis fill
axis ij
axis xy
axis equal
axis image 
axis square
axis vis3d
axis normal
axis off
axis on


aesthetics people, aesthetics...

and another reason that matlab beats the pants off of excel.

Thursday, January 21, 2010

colors

colors are important if you want your MATLAB plots to be aesthetically pleasing. which you should. because it makes life more enjoyable. i mean, see how the colors of the header image and the titles and links sort of match each other on this page, that is NOT by accident. and it pleases me.

another advantage of specifying your own colors is that you can make a bunch of personal colors that convert to the correct gray-shade. that way you can make one set of images for presentations (colorful!) and easily convert it to grey shade for papers. ingenious!

here's a webpage with the RGB codes for many colors:
http://en.wikipedia.org/wiki/Web_colors

in matlab colors are specified with RGB values between 0 and 1, not 0 and 255. so simply specify the code given on that page, e.g. i,j,k, like so: [i/255 j/255 k/255]

example:
plot(.6,1.5,'o', 'MarkerEdgeColor', [0 0 0], 'MarkerFaceColor', [70/255 130/255 180/255], 'MarkerSize',11, 'LineWidth', 1.5);

steel blue, mmmmmm.

btw, "coulours" is one of my favorite songs by The Prodigy.

Wednesday, January 20, 2010

continue

this command reminds me of jay-z's song on to the next one.

continue: passes control to the next iteration of a for or while loop.

if [something that makes you want to skip the rest of the iteration go to the next one]
continue;
end

Tuesday, January 19, 2010

position and PaperPositionMode

There are lots of figure properties, including one i use to specify the size of screen figure output and also .jpg or .eps size. this is helpful if you care about the aspect ratio of your output, which you should.

figurehandle = figure('position', [x1 y1 x2 y2],  'PaperPositionMode', 'auto');

i usually start from the bottom left corner of the screen (x1 = 0, y1 = 0) with units in pixels. The 'auto' makes the printed output (.jpg, .eps, etc.) the same aspect ratio/size.

h1 = figure('position', [0 0 800 1100], 'PaperPositionMode', 'auto');

many other options! figure properties page.

Monday, January 18, 2010

datestr(now) is now

a very simple command to make a timestamp:

>> datestr(now)
ans =
18-Jan-2010 15:31:58

i like to write this string to my output files, especially if it is a file that i append to. then if you change something in the program on some date, you can look at lines before that date to see the difference. (or, if you make 11 versions of the same dumb program, you can always look at the output file and know which version it's from.)

i use these lines a lot to make the first 3 columns of output files. i guess it's overkill to output the date string for every single row/iteration, but whatever.

today = datestr(now);

fprintf(fid, '%s , %s, %d, ', today, filename, i);

the first line goes at the top of my program, and the second one is inside a loop for i.

Sunday, January 17, 2010

fit.m, an alternative to excel

it's the year 2010 and i finally figured out the fit.m function. why is this awesome? because it takes the place of the tedious excel fits i've used for over a decade. how many hours of my life have i spent clicking the chart button, scatter plot button, fit trendline, power fit, show equation, show R^2? very inefficient! so here is a way that i like much better*:

x and y are your vectors holding your x_i and y_i

A simple linear fit (polynomial order 1):
[cfun, gof, output] = fit(x,y,'poly1')

Finding the exponent of a power law fit:
xlog = log10(x);
ylog = log10(y);
[cfun, gof, output] = fit(xlog,ylog,'poly1')

cfun holds the fit parameters cfun(x) = p1*x + p2, which you can access by cfun.p1 or cfun.p2
gof holds goodness of fit parameters including rsquare and rmse (root mean square error)
output holds other parameters like numobs (# of observations), residuals, and more

i left off the ending semi-colons so the structures will appear automatically

or, you can output these to a text file like so:

fid = fopen('C:\path\filename.txt', 'a+');
 % open file for writing, append mode (will create file if it does not exist yet)

fprintf(fid, '%8.3f, %8.3f, %8.3f %8.3e, %d\n', cfun.p1, cfun.p2, gof.rsquare, gof.rmse, output.numobs); 
% the different formats used are f = fixed-point 8 digits, 3 decimals, e = exponential notation, d = integer, \n = newline

status = fclose('all');
% close the file properly so you can open it with other programs

this is especially useful when employed in a loop, say, if you want to look at a large data set and then certain subsets of that data, or if you want to look at y versus several different x variables.

here is a list of other fit types.
note/question: with the method above, i get the same power law fit parameters as using excel. there is a power fit option in matlab (instead of taking log10 and doing a polynomial fit) but i don't get the same answers.

now you can
(1) open the text file with wordpad (or other) to inspect
(2) use the file as input for matlab plot
(3) *sigh* open the file in excel as a comma delimited, save to .xls and plot

post questions or improvements to comments.

* i must admit that i probably wasn't using excel to its full capabilities by doing things the tedious way, but i just never learned any cool shortcuts in excel. that says something about me and something about excel.