SHAPE file for draw polyf
Graziano Giuliani
giuliani at LAMMA.RETE.TOSCANA.IT
Wed Dec 14 10:21:53 EST 2005
Dear all,
I have made a little hack to grads code (version 1.9b4 codebase, with some
patches applied of whom I did speak about in previous message of 4 November,
just some cosmetic to build against recent OpenDAP and some clean-up for
Debian install with gcc 4.0) to directly use a shapefile for the draw polyf
function.
The Shapefile format is a working and interchange format promulagated by ESRI
for simple vector data with attributes. The file format actually consists of
three files.
XXX.shp - holds the actual vertices.
XXX.shx - hold index data pointing to the structures in the .shp file.
XXX.dbf - holds the attributes in xBase (dBase) format.
There exist a simple C library to access data available from
http://shapelib.maptools.org/ maintained by Frank Warmerdam which is LGPL (or
MIT style, You can choose).
As a lot of "shapefiles" are available on the web, (freegis for example, and
also gshhs) I wanted to bypass the ASCII lpolyxxxx files and the basemap gs
script.
Sooo...
43a44,47
> #if USESHAPE == 1
> #include <shapefil.h>
> #endif
>
1081a1086,1166
> #if USESHAPE == 1
> if (i==3) {
> SHPHandle hSHP;
> int nShapeType, nshp, nobj, npart, pstart, nvert, ivert;
> double adfMinBound[4], adfMaxBound[4];
> float shplon, shplat;
> SHPObject *psShape;
> char shpname[1024];
> char message[2048];
>
> sscanf(nxtwrd(cmd), "%s", shpname);
> sprintf(message, "Using shapefile %s\n", gxgnam(shpname));
> gaprnt (0,message);
> hSHP = SHPOpen(gxgnam(shpname), "rb");
> if (hSHP == NULL)
> {
> gaprnt (0,"SHAPEFILE not found.\n");
> gaprnt (0,"DRAW error: Syntax is DRAW POLYF shapefile object
part\n");
> return (1);
> }
> SHPGetInfo(hSHP, &nshp, &nShapeType, adfMinBound, adfMaxBound);
> cmd = nxtwrd(cmd);
> cmd = nxtwrd(cmd);
> if (intprs(cmd, &nobj) == NULL)
> {
> gaprnt (0,"DRAW error: Invalid object number in shapefile\n");
> return (1);
> }
> cmd = nxtwrd(cmd);
> if (intprs(cmd, &npart) == NULL)
> {
> gaprnt (0,"DRAW error: Invalid part of object number\n");
> return (1);
> }
> if (nobj < 0 || nobj >= nshp)
> {
> gaprnt (0,"DRAW error: Invalid object number in shapefile\n");
> return (1);
> }
> psShape = SHPReadObject(hSHP, nobj);
> if (npart < 0 || npart >= psShape->nParts)
> {
> gaprnt (0,"DRAW error: Invalid part of object number\n");
> return (1);
> }
> pstart = psShape->panPartStart[npart];
> if (npart == psShape->nParts-1)
> {
> nvert = psShape->nVertices-pstart;
> }
> else
> {
> nvert = psShape->panPartStart[npart+1]-pstart;
> }
> xy = (float *)malloc(2*sizeof(float)*(nvert+2));
> if (xy==NULL) {
> gaprnt (0,"DRAW error: Memory allocation error\n");
> return (1);
> }
> i = 0;
> for (ivert = 0; ivert < nvert; ivert ++)
> {
> shplon = psShape->padfX[ivert+pstart];
> shplat = psShape->padfY[ivert+pstart];
> if ((shplon >= pcm->dmin[0] && shplon <= pcm->dmax[0]) &&
> (shplat >= pcm->dmin[1] && shplat <= pcm->dmax[1]))
> {
> gxconv (shplon,shplat,&x,&y,2);
> xy[i] = x;
> xy[i+1] = y;
> i += 2;
> }
> }
> if (i == 0)
> {
> gaprnt (0,"DRAW warning: Outside window.\n");
> return (0);
> }
> }
> else if (i<6) {
> #else
1082a1168
> #endif
1086,1096c1172,1185
< xy = (float *)malloc(sizeof(float)*(i+2));
< if (xy==NULL) {
< gaprnt (0,"DRAW error: Memory allocation error\n");
< return (1);
< }
< i = 0;
< while ( (cmd = nxtwrd(cmd)) != NULL) {
< if ( valprs(cmd,xy+i) == NULL ) {
< gaprnt (0,"DRAW error: Invalid polyf coordinate\n");
< free (xy);
< return(1);
---
> else {
> xy = (float *) malloc(sizeof(float)*(i+2));
> if (xy==NULL) {
> gaprnt (0,"DRAW error: Memory allocation error\n");
> return (1);
> }
> i = 0;
> while ( (cmd = nxtwrd(cmd)) != NULL) {
> if ( valprs(cmd,xy+i) == NULL ) {
> gaprnt (0,"DRAW error: Invalid polyf coordinate\n");
> free (xy);
> return(1);
> }
> i++;
1098d1186
< i++;
this way a new syntax is added to draw polyf:
draw polyf shapefile object partition
as in
draw polyf shapes/gshhs/gshhs_land 128 0
using shapefiles from
http://www.ngdc.noaa.gov/mgg/shorelines/data/gshhs/version1.3/shapefiles
Just in case I attach another simple program to be linked with same shapelib
to convert any shapefile or directly the gshhs format or any simple ascii
format into a map data for the GrADS program.
Cheers,
Graziano.
--
\ | /
(@ @)
-------------------------o00-(_)-00o -----------------------------
LaMMA - Laboratorio per la Meteorologia e la Modellistica Ambientale
Laboratory for Meteorology and Environmental Modelling
Via Madonna del Piano, 50019 Sesto Fiorentino (FI)
tel: + 39 055 4483049
fax: + 39 055 444083
web: www.lamma.rete.toscana.it
e-mail: giuliani at lamma.rete.toscana.it
-------------- next part --------------
A non-text attachment was scrubbed...
Name: create_map.c
Type: text/x-csrc
Size: 14238 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20051214/124fa573/attachment.bin
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Makefile
Type: text/x-makefile
Size: 463 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20051214/124fa573/attachment-0001.bin
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20051214/124fa573/attachment-0002.bin
More information about the gradsusr
mailing list