Linux - SoftwareThis forum is for Software issues.
Having a problem installing a new program? Want to know which application is best for the job? Post your question in this forum.

Notices

Welcome to LinuxQuestions.org, a friendly and active Linux Community.

You are currently viewing LQ as a guest. By joining our community you will have the ability to post topics, receive our newsletter, use the advanced search, subscribe to threads and access many other special features. Registration is quick, simple and absolutely free. Join our community today!

Note that registered members see fewer ads, and ContentLink is completely disabled once you log in.

If you have any problems with the registration process or your account login, please contact us. If you need to reset your password, click here.

Having a problem logging in? Please visit this page to clear all LQ-related cookies.

Introduction to Linux - A Hands on Guide

This guide was created as an overview of the Linux Operating System, geared toward new users as an exploration tour and getting started guide, with exercises at the end of each chapter.
For more advanced trainees it can be a desktop reference, and a collection of the base knowledge needed to proceed with system and network administration. This book contains many real life examples derived from the author's experience as a Linux system and network administrator, trainer and consultant. They hope these examples will help you to get a better understanding of the Linux system and that you feel encouraged to try out things on your own.

Click Here to receive this Complete Guide absolutely free.

I need to create a vector image from a gridded dataset. The image is a map of the world (or even a limited region) and data are regularly spaced topography data. I can easily create a raster image using Matlab or Ocean Data View, but my requirement is to create a vector image and I don't know which linux software (if any) has the capability to read/import gridded data from a file and plot them on a vector graphic. Any idea?

Assuming you mean contour plots, gnuplot might fit the bill. You can plot to a vector graphics file (SVG is my favourite).

Another alternative is to generate an SVG directly from the data e.g. using awk or Perl/Python script. First build a net (no need to do e.g. Delaunay triangulation, if the data is reasonably regular), then draw the curve(s) within each cell. Follow the net to connect the curves. It is not too difficult, if you first work it out with a smaller dataset.

Would you like to describe what you are trying to image in a bit more detail? Contours? Projections? Silhouettes?

Last edited by Nominal Animal; 09-10-2011 at 02:38 AM.

I need to create a vector image from a gridded dataset. The image is a map of the world (or even a limited region) and data are regularly spaced topography data. I can easily create a raster image using Matlab or Ocean Data View, but my requirement is to create a vector image and I don't know which linux software (if any) has the capability to read/import gridded data from a file and plot them on a vector graphic. Any idea?

Thank you.

Have a look at GRASS. This is a very capable GIS which should be able to handle your data.
Sorry if I have misunderstood your requirement.

Thank you Nominal Animal for your reply. I will take a look at Gnuplot. I like the idea of building the XML code using awk or some other scripting language, but I'm afraid it's far too difficult, since there are a lot of data and usually I don't plot them as contours, rather as indexed images where every pixel (that represents one cell of my grid) picks a color from the current colormap based on its value.

Quote:

Originally Posted by Nominal Animal

Would you like to describe what you are trying to image in a bit more detail? Contours? Projections? Silhouettes?

Please see an example in attachment. Here the size of the data matrix is 2820x960.

While GRASS and other GIS are obviously the best choice if they have the desired functionality -- no need to reinvent the wheel, and all that --, I wrote one possible approach in C99 from scratch.

(The reason I did this is twofold: One, it is not that difficult, and I wanted to show a real life example. Two, I wanted an algorithm that could produce the contour curves in 3D (4D data), for example for the full Earth elevation model.)

One function does the actual work. It takes a triangle mesh, and the samples at each vertex as parameters, and returns the contour curve(s). Let's look at the header file first.

The function will work just fine with 3D data (4D if you count the "elevation" or sample at each point); just define USE_3D . The between() inline function is responsible for producing a point on the contour curve when the curve intersects an edge. The two points are the edge endpoints, and t=0..1 describes where the contour line is in relation to the two point samples. Let me explain the isocurve() function interface in more detail further below; for now, let's look at the source code:

Code:

#include <stdlib.h>
#include <errno.h>
#include "isocurves.h"
struct edge {
point_t c; /* Intersected point */
size_t p[2]; /* Point indices identifying this edge */
size_t e; /* Next edge along the contour curve */
signed char s; /* Status: 3 = start, 2 = tracing, 1 = unused, 0 = used. */
};
path_t *isocurves(const size_t points,
const point_t point[],
const double value[],
const size_t tris,
const tri_t tri[],
const double level)
{
size_t edges_max = 0;
size_t edges = 0;
struct edge *edge = NULL;
path_t *root = NULL, *path;
size_t t, e, i, n;
/* We need at least one triangle and three points. */
if (points < (size_t)3 || tris < (size_t)1) {
errno = EINVAL;
return NULL;
}
/* Check that all triangle points are valid. */
t = tris; while (t--)
if (tri[t].p[0] >= points ||
tri[t].p[1] >= points ||
tri[t].p[2] >= points ||
tri[t].p[0] == tri[t].p[1] ||
tri[t].p[0] == tri[t].p[2] ||
tri[t].p[1] == tri[t].p[2]) {
errno = EDOM;
return NULL;
}
/* Find triangles bisected by the isocurve. */
t = tris; while (t--) {
const size_t p[3] = { tri[t].p[0], tri[t].p[1], tri[t].p[2] };
const int type = ( (value[p[0]] >= level) ? 1 : 0 )
| ( (value[p[1]] >= level) ? 2 : 0 )
| ( (value[p[2]] >= level) ? 4 : 0 );
size_t curr_p[2], curr_e;
size_t next_p[2], next_e;
/* Nothing in this triangle? */
if (type < 1 || type > 6)
continue;
/* We know there is a contour line bisecting this triangle,
* since the three points straddle the contour level.
* Orient the line segment so that the upper half is on the right.
*/
/* Determine which edge starts the contour line segment
* passing through this triangle. For efficiency,
* pick the points so that curr_p[0] < curr_p[1].
*/
switch (type) {
case 1: case 5:
if (p[0] < p[1]) {
curr_p[0] = p[0];
curr_p[1] = p[1];
} else {
curr_p[0] = p[1];
curr_p[1] = p[0];
}
break;
case 2: case 3:
if (p[1] < p[2]) {
curr_p[0] = p[1];
curr_p[1] = p[2];
} else {
curr_p[0] = p[2];
curr_p[1] = p[1];
}
break;
case 4: case 6:
if (p[0] < p[2]) {
curr_p[0] = p[0];
curr_p[1] = p[2];
} else {
curr_p[0] = p[2];
curr_p[1] = p[0];
}
break;
}
/* Determine which ends starts the contour line segment
* passing through this triangle. For efficiency,
* pick the points so that next_p[0] < next_p[1].
*/
switch (type) {
case 1: case 3:
if (p[0] < p[2]) {
next_p[0] = p[0];
next_p[1] = p[2];
} else {
next_p[0] = p[2];
next_p[1] = p[0];
}
break;
case 2: case 6:
if (p[0] < p[1]) {
next_p[0] = p[0];
next_p[1] = p[1];
} else {
next_p[0] = p[1];
next_p[1] = p[0];
}
break;
case 4: case 5:
if (p[1] < p[2]) {
next_p[0] = p[1];
next_p[1] = p[2];
} else {
next_p[0] = p[2];
next_p[1] = p[1];
}
break;
}
/* Increase edge buffer allocation, if necessary.
* we may add up to two edges at once.
*/
if (edges + (size_t)2 >= edges_max) {
const size_t max = edges + (size_t)4096;
struct edge *temp;
temp = realloc(edge, max * sizeof(struct edge));
if (!temp) {
if (edge)
free(edge);
errno = ENOMEM;
return NULL;
}
edge = temp;
edges_max = max;
}
/* Find the index of the edge where the bisecting line segment ends.
* If we create a new edge, calculate the point where the contour
* curve crosses the edge, too. (Mark the new edge active, "s=1", for later.)
*/
next_e = 0;
while (next_e < edges && (edge[next_e].p[0] != next_p[0] || edge[next_e].p[1] != next_p[1]))
next_e++;
if (next_e >= edges) {
/* Create new edge. */
next_e = edges++;
edge[next_e].p[0] = next_p[0];
edge[next_e].p[1] = next_p[1];
edge[next_e].e = next_e;
edge[next_e].c = between(point[next_p[0]], point[next_p[1]], (level - value[next_p[0]]) / (value[next_p[1]] - value[next_p[0]]));
edge[next_e].s = 1;
}
/* Find the index of the edge where the bisecting line segment starts.
* If we create a new edge, calculate the point where the contour
* curve crosses the edge, too. (Mark the new edge active, "s=1", for later.)
*/
curr_e = 0;
while (curr_e < edges && (edge[curr_e].p[0] != curr_p[0] || edge[curr_e].p[1] != curr_p[1]))
curr_e++;
if (curr_e < edges) {
/* Found; link. (This should always work.) */
if (edge[curr_e].e == curr_e)
edge[curr_e].e = next_e;
} else {
/* Create new edge. */
curr_e = edges++;
edge[curr_e].p[0] = curr_p[0];
edge[curr_e].p[1] = curr_p[1];
edge[curr_e].e = next_e;
edge[curr_e].c = between(point[curr_p[0]], point[curr_p[1]], (level - value[curr_p[0]]) / (value[curr_p[1]] - value[curr_p[0]]));
edge[curr_e].s = 1;
}
}
/* Scan thorough all active edges. .s has values:
* 0=used, 1=unused/active, 2=tracing, 3=starting point.
*/
e = edges; while (e--) if (edge[e].s == 1) {
/* Mark starting point. */
edge[e].s = 3;
/* Trace the curve; we need to know the length. */
n = 2; /* Extra points to allocate */
i = edge[e].e;
while (edge[i].s == 1) {
n++;
edge[i].s = 2;
i = edge[i].e;
}
/* Ignore points. */
if (n < (size_t)3) {
edge[e].s = 0;
continue;
}
/* Allocate a new sub-path. */
path = malloc(sizeof(struct path) + n * sizeof(struct point));
if (!path) {
/* Oops. Out of memory. */
if (edge)
free(edge);
while (root) {
path = root;
root = root->next;
free(path);
}
errno = ENOMEM;
return NULL;
}
/* Prepend to the singly-linked list of sub-paths. */
path->next = root;
root = path;
/* Retrace the curve, but this time copy and mark it too. */
n = 0;
i = e;
edge[i].s = 2;
while (edge[i].s == 2) {
path->point[n] = edge[i].c;
n++;
edge[i].s = 0;
i = edge[i].e;
}
/* If it is a closed curve, duplicate the starting point. */
if (i == e) {
path->point[n] = edge[e].c;
n++;
}
path->points = n;
}
/* The edges are no longer needed. */
if (edge)
free(edge);
/* Done. */
errno = 0;
return root;
}

The logic of the function is surprisingly simple, even if my implementation does need a lot of work before it is ready for real use.

Determine which triangles are bisected by the contour curve.
If a triangle has one vertex sample on the other side of the comparison level than the other vertex samples for that triangle, then it is bisected by the contour curve.

The contour curve can cross each edge only once.
(The edge may be, and usually is, shared by two triangles, but the curve only crosses the edge once.)
The result is that each triangle can only be bisected by the contour curve.

If you orient the bisection lines naturally, they automatically form chains.
Assume you stand on the triangle face that has the vertices numbered counterclockwise.
Orient each bisection line so that the upper half stays on the right.

Closed bisection line chains form rings, so it does not matter where you start.
If you start an open bisection line chain from the middle, you will just split it in two.
If you want, you can re-merge such chains as the last step, with very little effort. I don't bother.

As long as you define all triangles in the mesh the same way, the algorithm works.
If you have holes, or your mesh does not reach the minimum at the edges, then you won't get closed curves.
To get closed curves, just add a border around the mesh at minimum level. (In 3D, the mesh is continuous, and does not have such "edges", so the contour curves will be closed.)

In practice, each element in the tri array contains three numbers. They specify the three vertices as indexes to the point array. The triangles must be defined counterclockwise. In 3D this means you must orient the triangles so that their normals point outwards. The samples (for example elevation) are always in the value array. The level parameter defines the contour sample value.

The points parameter is just for error checking; we could just trust the user that the points the refer to in the tri array exist.

Implementation-wise, the code is obviously quite ugly; especially the temporary edge list could be maintained much more efficiently. (Edges are completely defined by the point indices of the two endpoints -- two different nonnegative integers.)

There are many ways the above function can be used (and implemented, too!). If you have random samples, you can do e.g. Delaunay triangulation. With polygon meshes -- and for example ordinary rectangular grids -- you can just split each cell/loop into triangles.

Well, why not try that? Here is a small main program for the above two files. This one takes the number of samples horizontally and vertically as parameters, plus one or more contour specifications: value:HTMLColor or justvalue#RRGGBB, where value is the desired contour level. The output is almost valid SVG 1.1. Each contour is filled, so you will want to specify smaller values first. The program will read the samples from standard input, and assumes they form a regular rectangular grid. The rectangles are split along alternating diagonals, and the program will even add an outer border of minimum values to get closed contour curves. The code:

Code:

#include <stdio.h>
#include <string.h>
#include <errno.h>
#include "isocurves.h"
int main(int argc, char *argv[])
{
size_t points = 0;
point_t *point = NULL;
double *value = NULL;
size_t max_tris = 0;
size_t tris = 0;
tri_t *tri = NULL;
path_t *path, *curr;
size_t offset, stride, i, n;
double min, val;
int arg, o;
long nx, ny, x, y;
char dummy;
if (argc < 4) {
fprintf(stderr, "\n"
"Usage: %s xpoints ypoints level:color ...\n"
"\n"
"This program will read xpoints*ypoints scalar values from\n"
"standard input, and output an SVG image, with specified\n"
"levels drawn as contour plots, and filled with the desired color.\n"
"Example level and color definitions:\n"
"\t3.141:#ffffff\tWhite contour at 3.141\n"
"\t-2.542#ff0000\tRed contour at -2.542\n"
"\t0.0012#00ff00\tGreen contour at 0.0012\n"
"\n"
"The resulting contours contain a lot of redundant information.\n"
"This program does NOT try to optimize the contour paths.\n"
"\n"
, argv[0]);
return 1;
}
/* Parse grid size. */
if (sscanf(argv[1], "%ld %c", &nx, &dummy) != 1 || nx < 2L) {
fprintf(stderr, "%s: Invalid number of points along the x axis.\n", argv[1]);
return 1;
}
if (sscanf(argv[2], "%ld %c", &ny, &dummy) != 1 || ny < 2L) {
fprintf(stderr, "%s: Invalid number of points along the y axis.\n", argv[2]);
return 1;
}
/* Reserve one grid point as a border. */
stride = (size_t)(nx + 2L);
offset = stride + (size_t)1;
/* Allocate memory for the points, data, and triangles. */
points = stride * (size_t)(ny + 2L);
max_tris = (size_t)2 * points;
point = malloc(points * sizeof(point_t));
value = malloc(points * sizeof(double));
tri = malloc(max_tris * sizeof(tri_t));
if (!point || !value || !tri) {
fprintf(stderr, "Cannot allocate enough memory.\n");
return 1;
}
/* Initialize minimum to some ridiculous value. */
min = 1.0e307;
/* Read samples. */
fprintf(stderr, "Reading samples for a %ldx%ld grid: ", nx, ny);
fflush(stderr);
for (y = 0L; y < ny; y++)
for (x = 0L; x < nx; x++)
if (fscanf(stdin, "%lf", &val) > 0) {
if (min > val) min = val;
value[offset + x + stride * y] = val;
point[offset + x + stride * y].x = (double)x;
point[offset + x + stride * y].y = (double)y;
} else {
fprintf(stderr, "Error reading sample %ldx%ld.\n", x + 1L, y + 1L);
return 1;
}
/* Draw border with the minimum value, so we get full cycles. */
for (y = 0L; y < ny; y++) {
value[offset - 1 + stride * y] = min;
point[offset - 1 + stride * y].x = -1.0;
point[offset - 1 + stride * y].y = (double)y;
value[offset + nx + stride * y] = min;
point[offset + nx + stride * y].x = (double)nx;
point[offset + nx + stride * y].y = (double)y;
}
for (x = -1L; x <= nx; x++) {
value[offset + x - stride] = min;
point[offset + x - stride].x = (double)x;
point[offset + x - stride].y = -1.0;
value[offset + x + stride * ny] = min;
point[offset + x + stride * ny].x = (double)x;
point[offset + x + stride * ny].y = (double)ny;
}
fprintf(stderr, "Done.\nGenerating triangular mesh: ");
fflush(stderr);
for (y = -1L; y < ny; y++)
for (x = -1L; x < nx; x++) {
const size_t ul = offset + x + stride * y;
const size_t ur = ul + 1;
const size_t dl = ul + stride;
const size_t dr = dl + 1;
if ((x ^ y) & 1) {
tri[tris].p[0] = ul;
tri[tris].p[1] = ur;
tri[tris].p[2] = dl;
tris++;
tri[tris].p[0] = ur;
tri[tris].p[1] = dr;
tri[tris].p[2] = dl;
tris++;
} else {
tri[tris].p[0] = ul;
tri[tris].p[1] = dr;
tri[tris].p[2] = dl;
tris++;
tri[tris].p[0] = ul;
tri[tris].p[1] = ur;
tri[tris].p[2] = dr;
tris++;
}
}
fprintf(stderr, "%ld triangles.\n", (long)tris);
fflush(stderr);
printf("<?xml version=\"1.0\" standalone=\"no\" ?>\n"
"<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n"
"<svg width=\"20cm\" height=\"%.1fcm\" viewBox=\"0 0 %ld %ld\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n",
20.0*ny/nx, nx, ny);
for (arg = 3; arg < argc; arg++) {
if (sscanf(argv[arg], "%lf%n", &val, &o) < 1) {
fprintf(stderr, "%s: Cannot parse sample level. Skipping.\n", argv[arg]);
fflush(stderr);
continue;
}
while ((argv[arg][o] > 0 && argv[arg][o] < 33) ||
(argv[arg][o] == ':' || argv[arg][o] == '+' || argv[arg][o] == '/'))
o++;
if (!argv[arg][o]) {
fprintf(stderr, "%s: No fill color specified. Skipping.\n", argv[arg]);
fflush(stderr);
continue;
}
fprintf(stderr, "Generating contour %.3g (%s): ", val, (const char *)&argv[arg][o]);
fflush(stderr);
path = isocurves(points, point, value, tris, tri, val);
if (!path) {
fprintf(stderr, "%s.\n", strerror(errno));
fflush(stderr);
continue;
}
printf("<path d=\"");
while (path) {
curr = path;
path = path->next;
n = curr->points - (size_t)1;
printf("M %.5g,%.5g L", curr->point[0].x, curr->point[0].y);
for (i = 1; i < n ; i++)
printf(" %.5g,%.5g", curr->point[i].x, curr->point[i].y);
if (curr->point[0].x == curr->point[n].x && curr->point[0].y == curr->point[n].y)
printf(" Z");
else
printf(" %.5g,%.5g", curr->point[n].x, curr->point[n].y);
if (path)
putchar('\n');
free(curr);
}
printf("\" fill-rule=\"evenodd\" stroke=\"none\" fill=\"%s\" title=\"%.5g\" />\n", (const char *)&argv[arg][o], val);
fflush(stdout);
fprintf(stderr, "Done.\n");
fflush(stderr);
}
printf("</svg>\n");
return 0;
}

To compile, you can use e.g.

Code:

gcc -Wall -pedantic -std=c99 -O3 -o test main.c isocurves.c

and you can use e.g. awk to generate some test data:

The image is a map of the world (or even a limited region) and data are regularly spaced topography data.

grass is good
I do a lot of work using ISIS3 from NASA and USGS
one good tool for topography data is GDAL ( gdal_translate is the conversion tool ) http://www.gdal.org/

Hi nominal! A little update on this topic: I've launched your code more than two weeks ago and it's still running! I have 64 levels to compute and it's calculating the 41st now! The only inconvenience is that the output file size is more than 500M right now and I'm afraid I can't not manage such a huge object with standard graphic software and "only" 8G of RAM. Let's see at the end.

@John_VV: thanks for the suggestion. I will try a GIS software capable of exporting maps in SVG.

Thanks John. I'm using the ETOPO1 dataset from NOAA to generate my maps. It has an even too high resolution for my requirements. Anyway, good to know about the existence of another one, I will add it to my references/bookmarks.

Hi nominal! A little update on this topic: I've launched your code more than two weeks ago and it's still running!

Oh! If I'd known that, I'd written some hefty speed improvements. Specifically, the edge pairing is brute-force lookup, most likely quadratic (O(N^2)) right now. A list of edges per point would yield nearly O(1) behaviour. I suspect the edge search is the one that is eating all that CPU time for you right now.

I would like to enhance the program, but unfortunately I'm going almost-offline for two weeks in a couple of days. It certainly is an interesting problem!

Quote:

Originally Posted by colucix

I have 64 levels to compute and it's calculating the 41st now! The only inconvenience is that the output file size is more than 500M right now and I'm afraid I can't not manage such a huge object with standard graphic software and "only" 8G of RAM.

Each level is a single SVG path element. I recommend you first split each one into a separate SVG. They are trivial to merge back in afterwards. Each path can be split at "Z\nM" by restarting the path SVG element where the newline is; the coordinates are in the d attribute, and you can just duplicate all the other attributes. Each of those split paths will be a closed polygon, and can be manipulated independently of the others. For example, you might use inkscape to simplify the split paths.

I personally have a physics background, and I have an extremely deep-seated need to know (or be able to estimate) the error bounds, and errors introduced by the theories applied, and by my implementation. This is why I tend to do this kind of analysis from scratch. I do not know whether the various GIS explain in detail the methodology they use when computing the contour curves, and what kind of approximations and optimizations they do. Software that is designed to yield visually pleasing output is often not that accurate. Looking at standard maps, I do not have much confidence in those. (Personally, I'd also prefer to work on the original triangle meshes, not the regular rectangular grids extrapolated from that.)

I've considered using the slope at the contour curve points to determine bounds for simplifying the final polygons without compromising accuracy. Basically, the flatter the slope is, the more error there is in the approximation of the contour curve anyway -- although the contour curve will always be within the two original samples, one being above and the other below the contour curve. The nature of the samples (point samples, or averaged over a specific area around the point) affects this, and I'd need to brush up on various bits of theory first to get a real grip on that.

I wish I'd known your dataset takes that long to process. For one, it would have made much more sense to use a binary output format, with each polygon of each level in a separate file. Those are trivial to convert to SVG, but you could have experimented with the polygon optimization before converting the data to SVG. Right now, the rounding the text output uses may be a problem. You could also have started that immediately, when the first level (or even first polygon of first level) completed. (Perhaps you could have used Inkscape to simplify the polylines, and calculated the actual error that introduces.)

A complete data for each polygon is about eight (eleven) doubles per polyline: (x,y) for the next polyline endpoint, (x0,y0,z0) for the sample above the contour and (x1,y1,z1) for the sample below the contour. The two latter points are of course the original samples that yielded the point (x,y). If you have 3D points with a separate variable being evaluated, then it is eleven doubles per polyline (x,y,z,x0,y0,z0,v0,x1,y1,z1,v1). The contour level (elevation) is always the same for the entire polygon, so there is no need to save that for each point. Although the isocurves() function does not provide that in the result, it would be trivial to modify it to write the output directly to say basename-contourlevel.polygonnumber , as well as a helper program to convert that to SVG. (The inverse, converting back from SVG, is a bit more work, due to SVG having elements that transform and translate the path data.)

I'm quite interested in what you're doing, so I wouldn't mind exchanging ideas and suggestions about the way the data is processed; if you don't think it is suitable for LQ, you can always e-mail me directly.

LinuxQuestions.org is looking for people interested in writing
Editorials, Articles, Reviews, and more. If you'd like to contribute
content, let us know.