Reading numbers from a string in MATLAB

I always have needs to find numbers from a formatted string using Matlab. For example, find the lats and longs for all the grid cells of a river basin based on a bunch of input files, something like (forcings_35.25_-100.75). Then in Matlab you have two ways to read these numbers out.

  • using ‘sscanf’. It’s kinda stupid that you can’t find the two numbers by
  • sscanf('forcings_35.25_-100.75','%s_%f_%f')

    instead, you will get an 1d array that contains 22 numbers, with each number representing a single character in that string. What you need to do is using

  • sscanf('forcings_35.25_-100.75','forcings_%f_%f')

    I guess the reason behind this is that sscanf can only output one numerical matrix that based on the first character it processes. So if you do

  • sscanf('35.25_-100.75_forcings','%f_%f_%s')

    you will get the first two number you want, followed by 8 numbers representing the string “forcings”.

  • Another way is to use “textscan”. I feel safer using this method.
  • A = textscan('forcings_35.25_-100.75','%s %f %f', 'Delimiter', '_')

    This way, you will end up getting a three cell array representing the three parts separated by “_” in the string. Then simply using cell2mat() to convert the numbers to regular array. To convert “forcings” to string, using

  • char(A{1})
Advertisements

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.

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 directly read in USGS flow by given site number

Python packages are very powerful. Here is an example for hydrological data analysis. This small piece of code reads USGS flow by a given site number and store it as pandas DataFrame. You need to install “ulmo” first.

Original post is from here

import pandas as pd
import numpy as np
import ulmo

def importusgssite(siteno):
    sitename = {}
    sitename = ulmo.usgs.nwis.get_site_data(siteno, service="daily", period="all")
    sitename = pd.DataFrame(sitename['00060:00003']['values'])
    sitename['dates'] = pd.to_datetime(pd.Series(sitename['datetime']))
    sitename.set_index(['dates'],inplace=True)
    sitename[siteno] = sitename['value'].astype(float)
    sitename[str(siteno)+'qual'] = sitename['qualifiers']
    sitename = sitename.drop(['datetime','qualifiers','value'],axis=1)
    sitename = sitename.replace('-999999',np.NAN)
    sitename = sitename.dropna()
    #sitename['mon']=sitename.index.month
    return sitename

d = importusgssite('12472800');
d.plot(style='r', linewidth=1.0)

flowseries.png

Plotting in Python

Finally start to use python for plotting.

  • Plotting in notebook

ipython notebook is good and plots can be lined up with scripts. But if I need to run the same script on cluster without Xming, it will crash. Simple way is to add

import matplotlib as mpl
mpl.use('Agg')

when importing libraries.

  • Seaborn

Seaborn is a great package for color scheme picking. Color schemes are called “palettes” in seaborn. There are many pre-set palettes for you to choose. If you want to use “colorblind”, just simply add

import seaborn as sns
with sns.color_palette("colorblind"):

before the plotting script.

Customizing the palette is easy too:

colors = ["windows blue", "amber", "greyish", "faded green", "dusty purple"]
with sns.xkcd_palette(colors):

or

flatui = ["#9b59b6", "#3498db", "#95a5a6", "#e74c3c", "#34495e", "#2ecc71"]
sns.set_palette(flatui)
with sns.color_palette():

 

Tecplot in cmd

It’s kinda weird to use command lines in Windows (maybe a lot of people do it, who knows) but the reason I’m using it is that in Windows based Matlab, the “system” function will just call cmd. So if I need to run something inside a Matlab loop, I should learn how to use cmd.

Here is the command line I used to open a data file and run a python script in Tecplot:

C:\ /WAIT tec360 -b -datafile C:\something.plt -p C:\some\python.py

“/WAIT” is to hold the program until it’s done
“-b” is batch mode, include this otherwise you will see the GUI
“-datafile” is reading file
“-p” is executing python script
Note that the python script must be put in the python path otherwise Tecplot won’t find it.