LinuxQuestions.org
Share your knowledge at the LQ Wiki.
Go Back   LinuxQuestions.org > Forums > Linux Forums > Linux - Software
User Name
Password
Linux - Software This 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


Reply
  Search this Thread
Old 09-09-2011, 03:30 PM   #1
colucix
LQ Guru
 
Registered: Sep 2003
Location: Bologna
Distribution: CentOS 6.5 OpenSuSE 12.3
Posts: 10,509

Rep: Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981
Creating vector image from gridded data


Hi all,

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.
 
Old 09-10-2011, 02:30 AM   #2
Nominal Animal
Senior Member
 
Registered: Dec 2010
Location: Finland
Distribution: Xubuntu, CentOS, LFS
Posts: 1,723
Blog Entries: 3

Rep: Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948
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.
 
Old 09-10-2011, 03:02 AM   #3
terry-duell
Member
 
Registered: Jan 2007
Location: Melbourne, Australia
Distribution: Fedora 32 x86_64
Posts: 514

Rep: Reputation: 47
Quote:
Originally Posted by colucix View Post
Hi all,

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.

Cheers,
Terry
 
Old 09-10-2011, 03:50 AM   #4
colucix
LQ Guru
 
Registered: Sep 2003
Location: Bologna
Distribution: CentOS 6.5 OpenSuSE 12.3
Posts: 10,509

Original Poster
Rep: Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981
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 View Post
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.
Attached Thumbnails
Click image for larger version

Name:	mediterraneo.jpg
Views:	14
Size:	64.6 KB
ID:	7965  
 
Old 09-10-2011, 03:51 AM   #5
colucix
LQ Guru
 
Registered: Sep 2003
Location: Bologna
Distribution: CentOS 6.5 OpenSuSE 12.3
Posts: 10,509

Original Poster
Rep: Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981
Quote:
Originally Posted by terry-duell View Post
Have a look at GRASS. This is a very capable GIS which should be able to handle your data.
Thanks, Terry! If GIS software is capable of handling SVG, maybe that's a solution. I will give it a try.
 
Old 09-11-2011, 07:00 AM   #6
Nominal Animal
Senior Member
 
Registered: Dec 2010
Location: Finland
Distribution: Xubuntu, CentOS, LFS
Posts: 1,723
Blog Entries: 3

Rep: Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948
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.
Code:
#ifndef   ISOCURVES_H
#define   ISOCURVES_H
#include <stdlib.h>	/* For size_t */

typedef struct point	 point_t;
struct point {
	double		 x;
	double		 y;
#ifdef USE_3D
	double		 z;
#endif
};

typedef struct tri	 tri_t;
struct tri {
	size_t		 p[3];
};

typedef struct path	 path_t;
struct path {
	struct path	*next;
	size_t		 points;
	struct point	 point[];
};

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);

static inline point_t between(const point_t p0, const point_t p1, const double t)
{
	if (t <= 0.0)
		return p0;

	if (t >= 1.0)
		return p1;

#ifdef USE_3D
	return (struct point){ .x = p0.x + t * (p1.x - p0.x),
	                       .y = p0.y + t * (p1.y - p0.y),
	                       .z = p0.z + t * (p1.z - p0.z) };
#else
	return (struct point){ .x = p0.x + t * (p1.x - p0.x),
	                       .y = p0.y + t * (p1.y - p0.y) };
#endif
}

#endif /* ISOCURVES_H */
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.
  1. 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.
  2. 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.
  3. 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.
  4. 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.
  5. 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:
Code:
N=50

awk -v n=$N 'BEGIN {
	u=8.0/n
	for (y=0;y<n;y++)
		for (x=0;x<n;x++)
			printf("%.5g\n", cos((x*u-4.0)**2)*cos((y*u-4.0)**2))
}' | ./test $N $N '-0.5#900' '0.0#f00' '0.5#ff0' '0.75#fff' > test.svg

Last edited by Nominal Animal; 09-11-2011 at 07:04 AM. Reason: Increased N=50 in the final example.
 
1 members found this post helpful.
Old 09-11-2011, 04:55 PM   #7
John VV
LQ Muse
 
Registered: Aug 2005
Location: A2 area Mi.
Posts: 17,502

Rep: Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617
Quote:
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/

GRASS uses gdal and proj4
 
Old 09-28-2011, 03:19 AM   #8
colucix
LQ Guru
 
Registered: Sep 2003
Location: Bologna
Distribution: CentOS 6.5 OpenSuSE 12.3
Posts: 10,509

Original Poster
Rep: Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981
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.
 
Old 09-28-2011, 03:49 AM   #9
John VV
LQ Muse
 
Registered: Aug 2005
Location: A2 area Mi.
Posts: 17,502

Rep: Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617Reputation: 2617
for the whole Earth the best topo data set is the nasa's "srtm_ramp2"
http://mirrors.arsc.edu/nasa/topogra...0x43200.bin.gz

but you have been at it for a week ... so
 
Old 09-28-2011, 03:55 AM   #10
colucix
LQ Guru
 
Registered: Sep 2003
Location: Bologna
Distribution: CentOS 6.5 OpenSuSE 12.3
Posts: 10,509

Original Poster
Rep: Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981Reputation: 1981
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.
 
Old 09-28-2011, 01:51 PM   #11
Nominal Animal
Senior Member
 
Registered: Dec 2010
Location: Finland
Distribution: Xubuntu, CentOS, LFS
Posts: 1,723
Blog Entries: 3

Rep: Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948Reputation: 948
Quote:
Originally Posted by colucix View Post
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 View Post
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.
 
  


Reply


Thread Tools Search this Thread
Search this Thread:

Advanced Search

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off



Similar Threads
Thread Thread Starter Forum Replies Last Post
creating live usb image using moblin image creator 2 mobquest Linux - Newbie 0 07-09-2009 12:35 AM
convert image into a vector Chris106 Linux - Newbie 9 11-30-2008 07:35 AM
Image manager with vector support Hellmark Linux - Software 4 03-19-2008 02:01 AM
LXer: Creating a dd/dcfldd Image Using Automated Image & Restore (AIR) LXer Syndicated Linux News 0 03-07-2007 06:16 PM
creating data cd image problems with dd command gobabushka Linux - General 1 12-22-2006 12:33 PM

LinuxQuestions.org > Forums > Linux Forums > Linux - Software

All times are GMT -5. The time now is 12:24 PM.

Main Menu
Advertisement
My LQ
Write for LQ
LinuxQuestions.org is looking for people interested in writing Editorials, Articles, Reviews, and more. If you'd like to contribute content, let us know.
Main Menu
Syndicate
RSS1  Latest Threads
RSS1  LQ News
Twitter: @linuxquestions
Open Source Consulting | Domain Registration