GIS projections in Matlab and R

In spatial analysis it is very common to gather GIS layers from a variety of sources and more than often the layers may have different coordinate systems. Every GIS software know how to project each layer that comes with different projection information (The projection information is saved in a separate file with suffix *.prj) and can perform spatial analysis tasks on layers with different projection system. 

In Matlab and R however the spatial analysis operations work only if the layers have the same coordinate system. Luckily both languages provide functions to convert between projections. Matlab introduced some functions that help with the transformation in the latest update 2020b. So the following will not work in any previous Matlab version not even in 2020a.

Matlab

First let’s bring a couple of data into Matlab workspace. The first layer is a point layer with well locations and the second layer contains the Subregions of the C2Vsim model groundwater basins. For the first layer the *.prj file is missing therefore the coordinate is unknown. For the second file we know from the report that the coordinate system  is the EPSG:26910 – NAD83 / UTM zone 10N – Projected.
Mode info about EPSG see here 

% For the well layer
% Read the shapefile information
wells_info = shapeinfo('CVwelldata');
% Read the shapefile geometry
wells = shaperead('CVwelldata');

% For the Subregion layer
% Read the shapefile information
C2Vsim_basins_info = shapeinfo('C2Vsim_subBasins');
% Read the shapefile geometry
C2Vsim_basins = shaperead('C2Vsim_subBasins');

The wells_info.CoordinateReferenceSystem returns an empty field because the projection information is missing. On the other hand  C2Vsim_basins_info.CoordinateReferenceSystem returns the following structure

ans = 

  projcrs with properties:

                    Name: "NAD83 / UTM zone 10N"
           GeographicCRS: [1×1 geocrs]
        ProjectionMethod: "Transverse Mercator"
              LengthUnit: "meter"
    ProjectionParameters: [1×1 map.crs.ProjectionParameters]

Unfortunately I’m not aware of any automated process to identify the projection system but we can make guesses by examining the coordinate values. If we plot the well data set we get the following image

plot([wells.X]', [wells.Y]','.')
axis equal

We can see that the range of X and Y coordinates is between -124 – -114 and 33 – 42. If you work with California data this range corresponds to the WGS 84 system with EPSG:4326. This is the coordinate system that many web applications use such as google maps, earth engine etc. use. We can easily verify our assumptions in matlab using the webmap tool. If the coordinate system is indeed 4326 then they should be projected in the right place.

webmap % starts the webmap
wmmarker([wells.Y]', [wells.X]','Alpha',1, 'Color','green', 'IconScale', 0.2, ...
  'Icon','http://maps.google.com/mapfiles/kml/shapes/shaded_dot.png');

We can see that the well locations fall in the right place. Now we can transform this to any coordinate system using. Note that you may have to switch X and Y

[Wx,Wy] = projfwd(projcrs(3310),[wells.Y]', [wells.X]');

Next we have to convert the basin coordinates from EPSG:26910 to EPSG:3310. I haven’t found yet and way to do this transformation using one line so I first convert from EPSG:26910 to EPSG:4326 using projinv and then convert to EPSG:3310. Also because the projinv and projfwd accept only X,Y coordinates not structures we have to loop through the polygons

for ii = 1:length(C2Vsim_basins)
    [lat,lon] = projinv(projcrs(26910),C2Vsim_basins(ii,1).X, C2Vsim_basins(ii,1).Y);
    [C2Vsim_basins(ii,1).X3310, C2Vsim_basins(ii,1).Y3310] = projfwd(projcrs(3310),lat, lon);
end

Now that the two layers have their coordinates on the same system we can overlay them

plot(Wx,Wy,'.')
hold on
for ii = 1:length(C2Vsim_basins)
    plot(C2Vsim_basins(ii,1).X3310, C2Vsim_basins(ii,1).Y3310,'r')
end

 

Now that the data sets share the same coordinate system we can do spatial operations. For example we can isolate the wells that are located in the southern most sub-basin of the Central Valley.

in = inpolygon(Wx, Wy, C2Vsim_basins(21,1).X3310, C2Vsim_basins(21,1).Y3310);
plot(Wx(in), Wy(in),'.')
hold on
plot(C2Vsim_basins(21,1).X3310, C2Vsim_basins(21,1).Y3310,'r')

 

R

In R the spatial transformations are even easier. First load the library library("rgdal")

The transformation between coordinate system is just one line of code

# Read the data sets
wells <- readOGR(dsn = '.', layer = 'CVwelldata')
C2Vsim_basins <- readOGR(dsn = '.', layer = 'C2Vsim_subBasins')

# Change the coordinate system to 3310
C2Vsim_basins_3310 <- spTransform(C2Vsim_basins, CRS("+init=epsg:3310"))
# Because the wells have no coordinate information first we have to assign the correct one
# we know that this is the 4326
crs(wells) <- CRS('+init=EPSG:4326')
# Now we can make the transformation
wells_3310 <- spTransform(wells, CRS("+init=epsg:3310"))

 

 

Customize Matlab color order

It’s very common to have to plot a number of individual lines on the same figure. In Matlab this can be as easy as plot(t, Data) which results in the following figure:

As you can see the default color order includes 7 colors which are repeated.  You can get information about the default colors using get(gca,'Colororder') which returns a matrix 7×3 with the rgb channels of the above colors.

There are multiple ways to define a custom color order. One I have found and works well for me is the following. Instead of just plotting, get a handle for the plot

h = plot(t, Data);

Then use the plot handle to change the color order

set(h, {'color'}, num2cell(colormat,2));

where colormat  is a matrix Nx3  where N is the number of unique colors.

Of course, the selection of colors is important and not an easy task. The approach I suggest is the following which makes it quite trivial.

Navigate to http://colorbrewer2.org/. Choose the number of classes and the nature of  your data. In my case I choose qualitative and 11 classes. Then from the dropdown menu choose RGB to display the RGB values of colors. Note that these are scaled between 0-255, while matlab requires scaling between 0-1.

A handy feature in this website is that you can copy/paste the values to an editor and convert it as matlab variable with a minimum effort. The result will look like the following:

colormat = [ ... 
    166,206,227; ... 
    31,120,180; ... 
    178,223,138; ... 
    51,160,44; ... 
    251,154,153; ... 
    227,26,28; ... 
    253,191,111; ... 
    255,127,0; ... 
    202,178,214; ... 
    106,61,154; ... 
    255,255,153... 
]/255;

Using the above color order matrix will make the plot to use a different color for each line

The http://colorbrewer2.org/ has some limits on the number of colors that you can define for each category. If you think you need  more colors then it’s better to question the type of plot you are trying to make.

Particle tracking in Houdini

In this post I’ll describe how one can do a simple particle tracking visualization in Houdini.

In general the velocity field will be given, as a product of a simulation for example. For this post we will create a random one

Generate random velocity field

First we create a geometry node and delete the default file node. Then we lay down three line nodes that were renamed as X,Y,Z discretization with the following properties:

Length:1, Points:5 (for start) and Direction {1,0,0}, {0,1,0}, {0,0,1} for the X,Y,Z line respectively.

Then copy the Z line to the points of X and then copy the result to the points of Y. So far we have a 3D grid of points.

Now we are going to assign random velocities on the nodes. Lay down an attribute wrangle node and copy the following code:

v@v = rand(@ptnum);
float w = fit(rand(@ptnum), 0, 1, 0.2, 0.8);
v@v = v@v*w;

The v@v  is read as “declare a vector with name v equal to”.The wrangle node is set to run over points and the @ptnum  is the point id which is used a seed to the rand function that generates a random.

The second line generates a random float value which is scaled between 0.2 and 0.8. Last the velocity vector is scaled using the w value.

We can visualize the velocity field using the visualize node. Connect the output of the attribute wrangle node to a new visualize node and change under the visualizers tab the type from color to marker the style to vector and choose v as attribute value. You may need to scale the vectors down a bit. (Note also that sometimes you need to go one level up to see the velocity field):

    

Last select all the nodes that are used for the random velocity field generation and create a subnet (shift+C) and rename it to Velocity field for example.

Generate Random initial particles

This is another process which most of the time is not random as the initial particle positions are defined by the problems. However here we will generate random initial particles within the space 0 1. This involves just on attribute wrangle node with the following code

int N = chi('Number_of_particles');
float seed = ch('Random_Seed');

for (int i = 0; i < N; i++){
    vector pos = rand( seed * (float(i)+1));
    addpoint(0, pos);
}

The first two lines generate GUI elements for this node. After inserting these two lines press the highlighted button which looks like sliders. This will generate two sliders below the code area.

With these two sliders we are able to change the number of particles and the random seed without writing code.

Velocity interpolation

Now its time to write the code for the velocity interpolation. Just to stay organized lay down a null node (which is doing absolutely nothing) and with the null node selected press Shift+C to create a subnet and rename the subnet as Particle tracer. Change also the input labels as follows:

and the network so far looks like the following:

Double click to enter the Particle tracer node and delete the null node.

The first input is expected to be a number of points with the initial particle positions. It is possible that these points form primitives which we don’t want. therefore we will make sure that the primitives are deleted. Lay down an attribute wrangle node, rename it to delete primitives and set its Run Over option to Primitives. There is only one line of code for this wrangle node:

removeprim(0, @primnum, 0);

Now that we are sure that there are no primitives, we will create one primitive per particle, a polyline of the trajectory of the particle within the velocity field. Add another attribute wrangle node, rename it to create polylines and leave the Run Over option to Points. Paste the following two lines which, for each point add a primitive and vertex. At this point the primitives cannot be displayed. we need at least two vertices to display a line.

int iprim = addprim(0, "polyline");
addvertex(0, iprim, @ptnum);

Setup solver node

Finally its time to do the interesting part. Lay down a solver node and double click to enter inside the solver node. What makes solver special is that it gives us access to the geometry at the previous time step. Connect the Prev_Frame node with a new attribute wrangle node. Rename the new node to find next position, make sure to change the Run Over at primitives and connect the second input of particle tracer to the second input of the new wrangle node. The network should look similar to this:

Since the find next position node will do all the work, there is going to be a lot of code here.

First we are going to create a slider that will help us choose the number of points to be taken into account during velocity interpolation and a slider to define at what distance the points will be considered. Here we also define a scale velocity slider that will be used at a later time. The if statement ensures that even if one uses wrong number of interpolation nodes the code will still work. however the interpolation method will be similar to nearest neighbor interpolation.

int Nint = chi('N_points_velocity_interpolation');
float search_radius = ch('Search_radius');
float scale_vel = ch('Velocity_scale');

if (Nint <= 0)
    Nint = 1;

Next step is to find how many vertices the primitive has.

int Nverts = primvertexcount(0, @primnum);

There should be at least one vertex and the following lines get the point index and the position of that last vertex in the primitive

int pt_last = primpoint(0, @primnum, Nverts-1);
vector pos_last = point(0,'P', pt_last);

Now we can search for Nint  nearest points within a sphere with radius search_radius

int closept[] = pcfind(1, "P", @P, search_radius, Nint);

The array closept contains the ids of the velocity points that satisfy the above criteria. After we initialize some help variables we loop through the velocity points. (Note that the code will use Inverse Distance Weighting Interpolation IDW)

For each velocity point get the velocity vector and calculate the distance of the current trajectory point with the velocity interpolation point.

float w[];
float sum_w = 0;
vector vel[];
foreach (int pt; closept){
    int out;
    vector interp_vel = pointattrib(1, 'v', pt, out);
    float d = distance(pos_last, point(1, 'P', pt));

According to IDW if the distance is 0 set directly the velocity equal to the velocity node. When this condition is true the code removes any other weight and velocity, appends only the values of this velocity node and exits the loop.

if (d == 0){
    resize(w, 0);
    resize(vel, 0);
    sum_w = 1;
    append(vel, interp_vel);
    append(w, 1);
    break;

if the distance is not zero, the weight of this velocity point is set equal to the inverse distance

else{
    append(w, 1/d);
    sum_w += 1/d; 
    append(vel, interp_vel);
}

Now we can calculate the velocity of the point

vector my_vel = {0,0,0};
for (int i = 0; i < len(w); ++i){
    my_vel.x = w[i]*vel[i].x;
    my_vel.y = w[i]*vel[i].y;
    my_vel.z = w[i]*vel[i].z;
}
my_vel = my_vel/sum_w;

Finally we can compute the position of the next point in the particle trajectory. Note that we multiply the calculated velocity with the scalar that has the effect of speeding up or slowing down the velocity.

vector new_pos = pos_last + my_vel*scale_vel;

Although the following is not needed, will help us during visualization. The first line sets the calculated velocity as point attribute of the previous point and calculates the normal of the previous point so that it points in the new position.

setpointattrib(0, 'v', pt_last, my_vel, 'set');
vector normal_pt = new_pos - pos_last;
setpointattrib(0, 'N', pt_last, normal_pt, 'set');

So far we have found the new position but we have to add it to the primitive.

int new_pt = addpoint(0, Nverts);
setpointattrib(0, 'P', new_pt, new_pos, 'set');
addvertex(0, @primnum, new_pt);

That concludes the code if this wrangle node. Here is the full code without interruptions

//Number of points involved in the velocity interpolation
int Nint = chi('N_points_velocity_interpolation');
float search_radius = ch('Search_radius');
float scale_vel = ch('Velocity_scale');

if (Nint <= 0)
    Nint = 1;

// find how many vertices exist in the primitive
int Nverts = primvertexcount(0, @primnum);
if (Nverts >= 1){
    
    // get the id of the last one
    int pt_last = primpoint(0, @primnum, Nverts-1);
    //i@temp = pt_last;
    
    // position of last vertex
    vector pos_last = point(0,'P', pt_last);
    
    // Find the position of the last point in the primitive
    int closept[] = pcfind(1, "P", pos_last, search_radius, Nint);
    //i@Nfind = len(closept);
    
    // Calculate the distace from each velocity point
    
    float w[];
    float sum_w = 0;
    vector vel[];
    foreach (int pt; closept){
        int out;
        // get the velocity of the pt point
        vector interp_vel = pointattrib(1, 'v', pt, out);
        // find the distance between the velocity point and the point in question
        float d = distance(pos_last, point(1, 'P', pt));
        
        // If the distance is zero then assign this velocity
        // ignore the other values and exit the loop
        if (d == 0){
            resize(w, 0);
            resize(vel, 0);
            sum_w = 1;
            append(vel, interp_vel);
            append(w, 1);
            break;
        }
        else{
            append(w, 1/d);
            sum_w += 1/d; 
            append(vel, interp_vel);
        }
    
    }
    
    // Next compute velocity as weighted average
    vector my_vel = {0,0,0};
    for (int i = 0; i < len(w); ++i){
        my_vel.x = w[i]*vel[i].x;
        my_vel.y = w[i]*vel[i].y;
        my_vel.z = w[i]*vel[i].z;
    }
    
    my_vel = my_vel/sum_w;
    // assign the velocity as attribute to the previous point
    setpointattrib(0, 'v', pt_last, my_vel, 'set');
    
    // find the next position
    vector new_pos = pos_last + my_vel*scale_vel;
    
    // calculate the normal vector for the last position
    vector normal_pt = new_pos - pos_last;
    setpointattrib(0, 'N', pt_last, normal_pt, 'set');
    
    // add the point in the primitive with id the next number
    int new_pt = addpoint(0, Nverts);
    setpointattrib(0, 'P', new_pt, new_pos, 'set');
    addvertex(0, @primnum, new_pt);
}

Make sure that the blue flag is set to the wrangle node and exit the solver node. Pressing play will show the particle trajectories as lines.

Particle trajectory appearance

The computation of particle tracking has finished in the previous section. Here we will work just on the appearance. First we will add color to the points according to velocity. To this end add another attribute wrangle node, rename it to set color and add the following lines. Press the button to generate the gui elements. This snippet calculates the velocity magnitude and assigns a color between green and red.

float min_vel = ch('Minimum_value');
float max_vel = ch('Maximum_value');
vector red = {0.9, 0, 0};
vector green = {0, 0.9, 0};
float l = length(@v);

@Cd = lerp(green, red, fit(l, min_vel, max_vel, 0.0, 1.0));

Convert polylines to tubes

Depending on the velocity field the polyline may be jagged. Therefore we lay down 2 nodes. A smooth node and a resample node. Both of these nodes have their own parameters that define the smoothing amount resampling amount. Its best to play around with those values to find what looks better.

It is very common to display the particle trajectories as tubes. To achieve this in Houdini is trivial.

First create a cirlce node and scale it down to something that makes sense. In my example I set the Uniform Scale parameter to 0.004.

Now we are going to loop through the primitives (at the moment the number of primitives is equal to the initial number of particles). However we will do that using the For-Each Primitive node. This is going to lay down two nodes. A for each begin and a for each end. This is going to loop through the primitives and execute whatever nodes are in between. Add the copy to points and skin nodes as shown in the figure:

Last we add a normal node which corrects the face normals so that the polygon faces are lighted correctly