source: subversion/applications/utils/srtm2shp/SRTMDataGrid.cpp @ 34385

Last change on this file since 34385 was 6182, checked in by nick, 12 years ago

new and internationally-capable srtm2shp added

File size: 9.9 KB
Line 
1/*
2    Copyright (C) 2005 Nick Whitelegg, Hogweed Software, nick@hogweed.org
3
4    This program is free software; you can redistribute it and/or modify
5    it under the terms of the GNU Lesser General Public License as published by
6    the Free Software Foundation; either version 2 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU Lesser General Public License for more details.
13
14    You should have received a copy of the GNU Lesser General Public License
15    along with this program; if not, write to the Free Software
16    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111 USA
17
18 */
19#include "SRTMDataGrid.h"
20//#include "functions.h"
21#include <cmath>
22#include "tomerc.h"
23
24
25// DataGrid constructor
26// Loads heights from one or more .hgt files.
27// Each rectangle in the input array "rects" defines the rectangle (in sampling
28// point indices) of a grid square. There's normally only one rectangle but
29// may be up to 4 if two latitude/longitude lines pass through the visible area.
30//
31// sampling_pts is a definition of the area in terms of nationally-indexed
32// sampling points; this function fills it in.
33//
34// Returns an array of all the heights indexed nationally.
35
36SRTMDataGrid::SRTMDataGrid(const std::string& srtmlocation,
37                                                        LATLON_TILE **rects,int w, int h, int f,
38                                                const std::string& outCoord, bool feet)
39{
40        // Do each input rectangle
41        int index_w=0, index_h=0,pts_w,pts_h;
42        samplewidth=0; 
43        sampleheight=0;
44
45        this->f=f;
46        this->outCoord=outCoord;
47        this->feet=feet;
48        this->srtmlocation=srtmlocation;
49
50        for(int hcount=0;hcount<h; hcount++)
51        {
52                sampleheight+=
53                                (rects[hcount][0].bottom-rects[hcount][0].top)+1;
54        }
55
56        for(int wcount=0; wcount<w; wcount++)
57        {
58                samplewidth+=
59                                (rects[0][wcount].right-rects[0][wcount].left)+1;
60        }
61
62        points=new SRTM_SAMPLE_POINT [sampleheight*samplewidth];
63
64        for(int hcount=0; hcount<h; hcount++)
65        {
66                index_w = index_h;
67                pts_h = (rects[hcount][0].bottom-rects[hcount][0].top)+1;
68                for(int wcount=0; wcount<w; wcount++)
69                {
70                        pts_w = (rects[hcount][wcount].right-rects[hcount][wcount].left)+1;
71                        doLoad(&rects[hcount][wcount],index_w);
72                        index_w += pts_w;
73                }
74                index_h += pts_h*samplewidth; 
75        }
76}
77
78
79
80void SRTMDataGrid::doLoad(LATLON_TILE *rect,int index)
81{
82        // Get the .hgt file for the current rectangle
83        char hgtfile[1024];
84        SRTMDataGrid::getHgtFilename(hgtfile,rect->origin);
85        FILE *fp=fopen(hgtfile,"rb");   
86
87        if(fp)
88        {
89
90                int width = (rect->right-rect->left)+1;
91                unsigned char *data = new unsigned char[width*2];
92
93                int datacount; 
94                double h;
95                EarthPoint curLatLon;
96
97                double frac = 0.0008333333333*f;
98                curLatLon.x=rect->origin.x+(((double)rect->left)/1200);
99                curLatLon.y=(rect->origin.y+1)-(((double)rect->top)/1200);
100                double origlong=curLatLon.x;
101                int i;
102
103                for(int row=rect->top; row<=rect->bottom; row++)
104                {
105                        curLatLon.x = origlong; 
106
107                        // 20/02/05 Only do every 'f' rows
108                        if((row-rect->top)%f == 0)
109                        {
110                                fseek(fp,(row*1201+rect->left)*2, SEEK_SET);
111                                fread (data,1,width*2,fp);
112                                datacount=0;
113                                i=0;
114                                for(int pt=row*1201+rect->left;pt<=row*1201+rect->right; pt++)
115                                {
116                                        // 20/02/05 Only do every 'f' columns
117                                        if( (pt-(row*1201+rect->left) )%f == 0) 
118                                        {
119                                                h=
120                                                ( ((double)data[datacount])*256+
121                                                  ((double)data[datacount+1]) ) * 
122                                                  (feet ? 3.28084: 1.0);
123                                                points[index+i].hgt = 
124                                                        (h>=1 && h<30000) ? h: 1;
125                                                points[index+i].earthPos=
126                                                                (outCoord=="Mercator" ? lltomerc(curLatLon):
127                                                                        curLatLon);
128
129                                                i++;
130                                        }
131                                        datacount+=2;
132                                        curLatLon.x +=  frac;
133                                }
134                                index+=samplewidth;
135                        }
136                        curLatLon.y -= frac;
137                }
138                delete[] data;
139                fclose(fp);
140        }
141}
142       
143// 211207 moved out of above method into its own as it is not always required,
144// e.g. for shapefiles
145void SRTMDataGrid::getScreenPoints(Map& map)
146{
147        for(int i=0; i<sampleheight*samplewidth; i++)
148                points[i].screenPos= map.getScreenPos(points[i].earthPos);
149}
150
151void SRTMDataGrid::getHgtFilename(char *hgtfile,EarthPoint& latlon)
152{
153        sprintf ( hgtfile,"%s/N%02d%s%03d.hgt",srtmlocation.c_str(),
154                                                int(latlon.y),(latlon.x<0 ? "W":"E"),
155                                                abs(int(latlon.x)));
156}
157
158
159SRTMDataGrid::~SRTMDataGrid()
160{
161        delete[] points;
162}
163
164
165void SRTMDataGrid::setPoint (int row,int col)
166{
167        pt = row*samplewidth+col;
168        edges[0][0] = pt;
169        edges[0][1] = pt+f;
170        edges[1][0] = pt;
171        edges[1][1] = pt+samplewidth*f;
172        edges[2][0] = pt+f;
173        edges[2][1] = pt+samplewidth*f+f;
174        edges[3][0] = pt+samplewidth*f;
175        edges[3][1] = pt+samplewidth*f+f;
176}
177
178double SRTMDataGrid::startHeight(int interval)
179{
180        double ht = min (
181                                        min(points[pt].hgt, points[pt+f].hgt),
182                                        min(points[pt+samplewidth*f].hgt, 
183                                                 points[pt+samplewidth*f+f].hgt )
184                                        );
185        return ceil(ht/interval) * interval;
186}
187
188double SRTMDataGrid::endHeight(int interval)
189{
190
191        double ht = max (
192                                        max(points[pt].hgt, points[pt+f].hgt),
193                                        max(points[pt+samplewidth*f].hgt, 
194                                                 points[pt+samplewidth*f+f].hgt )
195                                        );
196        return floor(ht/interval) * interval;
197}
198
199// checked : the correct edges are being considered each time.
200void SRTMDataGrid::getLine(LINE *lines,int *n_lines, int ht,bool screen)
201{
202        int go = 0;
203        int prevedges[2];
204        LINE line;
205        double eAh0, eAh1, eBh0, eBh1; 
206        double eAp0x, eAp0y, eAp1x, eAp1y, eBp0x, eBp0y, eBp1x, eBp1y;
207
208        // See addContour()
209        bool two_contours = (
210                        (points[pt].hgt<ht && points[pt+f].hgt>ht && 
211                         points[pt+samplewidth*f+f].hgt<ht && 
212                         points[pt+samplewidth*f].hgt > ht) 
213                         ||
214                        (points[pt].hgt>ht && points[pt+f].hgt<ht && 
215                         points[pt+samplewidth*f+f].hgt>ht && 
216                         points[pt+samplewidth*f].hgt < ht)
217                                ) ; 
218
219        for(int edge=0; edge<3; edge++)
220        {
221                if(between(ht,points[edges[edge][0]].hgt,points[edges[edge][1]].hgt))
222                {
223                        for(int edge2=edge+1; edge2<4; edge2++)
224                        {
225                                if(between
226                                        (ht,points[edges[edge2][0]].hgt,
227                                                 points[edges[edge2][1]].hgt))
228                                {
229                                        eAh0 = points[edges[edge][0]].hgt;
230                                        eAh1 = points[edges[edge][1]].hgt;
231                                        eBh0 = points[edges[edge2][0]].hgt;
232                                        eBh1 = points[edges[edge2][1]].hgt;
233
234                                        eAp0x = screen ? points[edges[edge][0]].screenPos.x:
235                                                                        points[edges[edge][0]].earthPos.x;
236                                        eAp0y = screen ? points[edges[edge][0]].screenPos.y:
237                                                                        points[edges[edge][0]].earthPos.y;
238                                        eAp1x = screen ? points[edges[edge][1]].screenPos.x:
239                                                                        points[edges[edge][1]].earthPos.x;
240                                        eAp1y = screen ? points[edges[edge][1]].screenPos.y:
241                                                                        points[edges[edge][1]].earthPos.y;
242                                        eBp0x = screen ? points[edges[edge2][0]].screenPos.x:
243                                                                        points[edges[edge2][0]].earthPos.x;
244                                        eBp0y = screen ? points[edges[edge2][0]].screenPos.y:
245                                                                        points[edges[edge2][0]].earthPos.y;
246                                        eBp1x = screen ? points[edges[edge2][1]].screenPos.x:
247                                                                        points[edges[edge2][1]].earthPos.x;
248                                        eBp1y = screen ? points[edges[edge2][1]].screenPos.y:
249                                                                        points[edges[edge2][1]].earthPos.y;
250
251
252                                        // We draw a line.
253                                        line.p[0].x =
254                                                        eAp0x + ( ((ht-eAh0) / (eAh1-eAh0))     
255                                                         * (eAp1x-eAp0x) );
256
257                                        line.p[0].y =
258                                                        eAp0y + ( ((ht-eAh0) / (eAh1-eAh0))     
259                                                         * (eAp1y-eAp0y) );
260
261                                        line.p[1].x =
262                                                        eBp0x + ( ((ht-eBh0) / (eBh1-eBh0))     
263                                                         * (eBp1x-eBp0x) );
264
265                                        line.p[1].y =
266                                                        eBp0y + ( ((ht-eBh0) / (eBh1-eBh0))     
267                                                         * (eBp1y-eBp0y) );
268
269                                        if(two_contours==false) 
270                                        {
271                                                lines[0] = line;
272                                                *n_lines=1;
273                                        }
274                                        else
275                                        {
276                                                SRTMDataGrid::addContour
277                                                                (lines,line,edge,edge2,prevedges,&go,
278                                                                 n_lines);
279                                        }
280                                }
281                        }
282                }
283        }
284}
285
286// This function will be called when the special case of two opposite corners
287// having a height above a contour, and the other two below, occurs. In this
288// case - and this case only - two contours of a given height will be drawn
289// through a quadrangle. This special case confuses the standard algorithm
290// no end :-) Thus, we need to ensure that the two contours are drawn
291// on the opposite side of the quadrangle (any other combination wouldn't
292// make sense), and this function does that. Which two opposite sides
293// doesn't actually matter; in fact, it's impossible to tell which two would
294// be correct.
295void SRTMDataGrid::addContour
296        (LINE *lines,LINE line,int edge1,int edge2,int *prevedges,int *go,
297        int *n_lines)
298{
299        /* Edge numbering assumed to be:
300
301             0
302           +---+
303          1|   |2
304           +---+
305             3
306
307        Don't draw contours through opposite edges. These don't apply in this
308        special case.
309
310        */     
311
312        if(!((edge1==0 && edge2==3) || (edge1==1 && edge2==2)))
313        {
314                if(*go==0)
315                {
316                        prevedges[0] = edge1;
317                        prevedges[1] = edge2;
318                        *go = 1;
319                        *n_lines = 1;
320                        lines[0] = line;
321                }
322                else if( (edge1==2 && edge2==3 && prevedges[0]==0 && prevedges[1]==1)
323                         || (edge1==1 && edge2==3 && prevedges[0]==0 && prevedges[1]==2)
324                          && *go==1)
325                {
326                        lines[1] = line;
327                        *n_lines = 2;
328                }
329        }
330}       
331
332// finds the minimum distance between a point and the list of previous points
333int SRTMDataGrid::hgtptDistance(vector<int>& prevs)
334{
335        int dist=sqrt(samplewidth*sampleheight), result; 
336        int dx, dy;
337        for(int count=0; count<prevs.size(); count++)
338        {
339                dx = prevs[count]%samplewidth - pt%samplewidth;
340                dy = prevs[count]/sampleheight - pt/sampleheight;
341                result= sqrt(dx*dx + dy*dy);
342                dist = (result<dist) ? result:dist;
343        }
344        return dist;
345}
346
347Colour SRTMDataGrid::getHeightShading (double shadingres)
348{
349        double est_ht = floor (((points[pt].hgt+points[pt+samplewidth*f+f].hgt)/2+
350                         (points[pt+f].hgt+points[pt+samplewidth*f].hgt)/2 
351                   ) / 2);
352        return doGetHeightShading(est_ht,shadingres);
353}
354
355Colour SRTMDataGrid::doGetHeightShading(double ht,double shadingres)
356{
357        Colour colour;
358        ht=round(ht/shadingres);
359        double a=1250/shadingres,
360        b=250/shadingres,
361        c=2000/shadingres,
362        d=500/shadingres,
363        e=1750/shadingres;
364       
365        colour.r = (ht<a) ? 191+(ht*(64/a)) : 
366        ((ht<c) ? 255-((ht-(1250/shadingres))*(16/b)) : 
367                                  255-((ht-d)*(16/d)) );
368
369        colour.g=(ht<c) ? 255-(ht*(16/b)) : 255-((ht+c)*(16/d));
370
371        colour.b = (ht<e) ? 191-(ht*(16/b)) :
372        ((ht<c) ? 95+((ht-e)*(16/b)) :
373                                  95+((ht-(1500/shadingres))*(16/d)) );
374        return colour;
375}
376
Note: See TracBrowser for help on using the repository browser.