source: subversion/applications/utils/srtm2shp/SRTMConGen.cpp @ 29022

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

new and internationally-capable srtm2shp added

File size: 9.8 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 "SRTMConGen.h"
20#include "tomerc.h"
21
22
23SRTMConGen::SRTMConGen(const std::string& srtmlocation,
24                                                EarthPoint& bottomleft,EarthPoint& topright,
25                                                bool feet,int f)
26{
27        makeGrid(srtmlocation,bottomleft,topright,feet,f);
28        inCoord=outCoord="latlon";
29}
30
31void SRTMConGen::makeGrid(const std::string& srtmlocation,
32                                                        EarthPoint& bottomleft,
33                                                        EarthPoint& topright,bool feet,int f)
34{
35        int w, h;
36
37        // Get the bounding rectangles for all lat/long squares
38        LATLON_TILE **tiles = get_latlon_tiles (bottomleft,topright,&w,&h);
39
40        // Get the sampled heights from the .hgt file, and the screen coordinates
41        sampledata = new SRTMDataGrid(srtmlocation,tiles,w,h,f,outCoord,feet);
42
43        for(int hcount=0; hcount<h; hcount++)
44                delete[] tiles[hcount];
45        delete[] tiles;
46}
47
48// Given the latitide and longitude of the bottom left and top right of the
49// visible area, this function returns an array of rectangles definining the
50// SRTM point indices for all latitude/longitude squares in the visible area.
51// This will normally be just one, but if, e.g., 51N and 1W both crossed the
52// visible area, it could be up to 4.
53LATLON_TILE ** SRTMConGen::get_latlon_tiles(EarthPoint& bottomleft,
54                                EarthPoint& topright,int *w,int *h)
55{
56
57        if (inCoord=="Mercator")
58        {
59                bottomleft = merctoll(bottomleft);
60                topright = merctoll(topright);
61        }
62
63        // Get the latitude/longitude square of each rectangle
64        LATLON_TILE **rect=getrects(bottomleft,topright,w,h);
65        EarthPoint llsq;
66
67        // Fill in the actual bounds
68        for(int hcount=0; hcount<*h; hcount++)
69        {
70                for(int wcount=0; wcount<*w; wcount++)
71                {
72                        llsq = rect[hcount][wcount].origin;
73                        rect[hcount][wcount].left =  bottomleft.x> llsq.?
74                                        floor((bottomleft.x - llsq.x)*1200): 0;
75                        rect[hcount][wcount].right = (topright.x < llsq.x+1) ?
76                                        1+floor((topright.x - llsq.x)*1200) : 
77                                                        1200;
78       
79                        rect[hcount][wcount].top =(topright.y < llsq.y+1 ) ?
80                                        floor(((llsq.y+1)-topright.y)*1200) :  0;
81                        rect[hcount][wcount].bottom = bottomleft.y > llsq.?
82                                        1+floor(((llsq.y+1)-bottomleft.y)*1200) :
83                                                1200;
84                }
85        }
86
87        return rect;
88}
89
90// Given the latitude and longitude of the bottom left and top right of the
91// visible map area, this function returns the appropriate number of
92// rectangles specific to a given grid square. For example, if both 51N and 1W
93// passed through the visible area, four rectangles would be generated, one
94// for the 51N/1W square, one for the 50N/1W square etc. Only the base latitude
95// and longitude are filled in at this stage. The dimensions of each rectangle
96// within the visible area are filled in by the calling function.
97LATLON_TILE** SRTMConGen::getrects
98        (const EarthPoint& bottomleft,const EarthPoint& topright, int *w,int *h)
99{
100        LATLON_TILE** rects;
101        *h=0;
102        rects=new LATLON_TILE* [int(floor(topright.y)-floor(bottomleft.y))+1];
103
104        for(int lat=floor(topright.y); lat>=floor(bottomleft.y); lat--)
105        {
106                *w=0;
107                rects[*h]=new LATLON_TILE
108                                [int(floor(topright.x)-floor(bottomleft.x))+1];
109                for(int lon=floor(bottomleft.x); lon<=floor(topright.x); lon++)
110                {
111                        rects[*h][*w].origin.y=lat;
112                        rects[*h][*w].origin.x=lon;
113                        (*w)++;
114                }
115                (*h)++;
116        }
117        return rects;
118}
119
120
121void SRTMConGen::generate(DrawSurface *ds, Map& map)
122{
123        sampledata->getScreenPoints(map);
124        std::map<int,vector<int> > last_pt;
125
126        for(int row=0; row<sampledata->getHeight()-1; row++)
127        {
128                // Do each point of the current row
129                for(int col=0; col<sampledata->getWidth()-1; col++)
130                {
131                        do_contours(ds,row,col,50, last_pt);
132                }
133        }
134}
135
136void SRTMConGen::do_contours (DrawSurface *ds,int row,int col, 
137                                int interval, std::map<int,vector<int> >&last_pt )
138{
139        sampledata->setPoint (row,col);
140        int start_ht = sampledata->startHeight(interval),
141                end_ht = sampledata->endHeight(interval);
142                       
143
144        LINE lines[2];
145        int n_line_pts;
146        char htstr[1024];
147
148        Colour colour, contour_colour(192,192,0), mint(0,192,64);
149
150        for(int ht=start_ht; ht<=end_ht; ht+=interval)
151        {
152                n_line_pts=0;
153                sampledata->getLine(lines,&n_line_pts,ht);
154
155
156                // draw line
157                if(n_line_pts!=0)
158                {
159                        for(int count=0; count<n_line_pts; count++)
160                        {
161                                colour = (ht%(interval*5)) ?  contour_colour : mint;
162                                if( (last_pt[ht].size()==0) || 
163                                                (sampledata->hgtptDistance(last_pt[ht])>20)) 
164                                {
165                                        // 08/02/05 changed parameters for slope_angle()
166                                        // 12/02/05 put all the text drawing code in angle_text()
167                                        double angle=slope_angle(lines[count].p[0].x, 
168                                                                        lines[count].p[0].y,
169                                                                        lines[count].p[1].x,
170                                                                        lines[count].p[1].y);
171                                        int i = (lines[count].p[1].x > lines[count].p[0].x) ? 0:1;
172                                        sprintf(htstr,"%d",ht);
173                                        ds->drawAngleText(8, angle, 
174                                                lines[count].p[i].x, lines[count].p[i].y, colour.r,
175                                                colour.g, colour.b, htstr);
176                                        last_pt[ht].push_back(sampledata->getPoint());
177                                }
178                                       
179                                ds->drawContour(lines[count].p[0].x,
180                                                lines[count].p[0].y,
181                                                lines[count].p[1].x,
182                                                lines[count].p[1].y, 
183                                                colour.r, colour.g, colour.b); 
184                        }       
185                }
186        }
187}
188
189
190// This method by Artem
191
192void SRTMConGen::merge_segments
193        (Contour & contour, LineString & line, unsigned maxlength)
194{
195   bool first = true;
196   unsigned max = contour.size()*contour.size();
197   unsigned i = 0;
198   double start_x;
199   double start_y;
200   double end_x;
201   double end_y;
202   
203   while (contour.size() && i++ < max)
204   {
205      LINE * seg = *(contour.end() - 1);
206      double x0 = seg->p[0].x;
207      double y0 = seg->p[0].y;
208      double x1 = seg->p[1].x;
209      double y1 = seg->p[1].y;
210      contour.pop_back();
211     
212      if (first)
213      {
214         first = false;
215         start_x = x0;
216         start_y = y0;
217         end_x = x1;
218         end_y = y1;
219         line.add_last(x0,y0);
220         line.add_last(x1,y1);
221      }
222      else if (start_x == x0 && start_y == y0)
223      {
224         start_x = x1;
225         start_y = y1;
226         line.add_first(x1,y1);
227      }
228      else if (start_x == x1 && start_y == y1)
229      {
230         start_x = x0;
231         start_y = y0;
232         line.add_first(x0,y0);
233      }
234      else if (end_x == x0 && end_y == y0)
235      {
236         end_x = x1;
237         end_y = y1;
238         line.add_last(x1,y1);
239      }
240      else if (end_x == x1 && end_y == y1)
241      {
242         end_x = x0;
243         end_y = y0;
244         line.add_last(x0,y0);
245      } 
246      else
247      {
248         contour.push_front(seg);
249      }
250      if (line.length() > maxlength) break;
251   }
252} 
253
254void SRTMConGen::generateShp (const char* shpname,int interval)
255{
256   SHPHandle shp = SHPCreate(shpname,SHPT_ARC);
257   DBFHandle dbf = DBFCreate(shpname);
258   int htidx = DBFAddField(dbf,"height",FTInteger,255,0);
259   int mjridx = DBFAddField(dbf,"major",FTInteger,255,0);
260   appendShp(shp,dbf,interval,htidx,mjridx);
261   DBFClose(dbf);
262   SHPClose(shp);
263}
264
265void SRTMConGen::appendShp(SHPHandle shp,DBFHandle dbf,int interval,
266                                                                int htidx, int mjridx)
267{
268   double xs[2], ys[2];
269   
270   SegmentCache segments;
271   Range heights;
272   for(int row=0; row<sampledata->getHeight()-1; row++)
273   {
274      // Do each point of the current row
275      for(int col=0; col<sampledata->getWidth()-1; col++)
276      {
277         sampledata->setPoint (row,col);
278         int start_ht = sampledata->startHeight(interval),
279            end_ht = sampledata->endHeight(interval);
280
281
282         LINE lines[2];
283         int n_line_pts;
284         char htstr[1024];
285
286         for(int ht=start_ht; ht<=end_ht; ht+=interval)
287         {
288            n_line_pts=0;
289            sampledata->getLine(lines,&n_line_pts,ht,false);
290            // draw line
291            if(n_line_pts!=0)
292            {
293               for(int count=0; count<n_line_pts; count++)
294               {
295                  segments.insert(std::make_pair(ht,lines[count]));
296                  heights.insert(ht);
297               }
298            }
299         }
300      }
301   }
302   
303   shape_writer writer(shp,dbf,htidx,mjridx,interval);
304   //std::for_each(segments.begin(),segments.end(),writer);
305   
306   Range::iterator itr = heights.begin();
307   Range::iterator end = heights.end();
308   
309   while (itr != end)
310   {
311      int height = *itr;
312      Contour contour;
313      SegmentCache::iterator pos;
314     
315//      std::cout << "height = " <<  height << "\n";
316      for ( pos = segments.lower_bound(height);
317            pos != segments.upper_bound(height);
318            ++pos)
319      {
320         contour.push_back(&(pos->second));
321      }
322     
323      while (contour.size() > 0)
324      {
325         LineString line(height);
326         merge_segments(contour,line,1000);
327         writer(line);
328      }
329      ++itr;
330   }
331   
332}
333
334void SRTMConGen::generateShading(DrawSurface *ds,double shadingres,Map& map)
335{
336        sampledata->getScreenPoints(map);
337        Colour colour;
338        for(int row=0; row<sampledata->getHeight()-1; row++)
339        {
340                // Do each point of the current row
341                for(int col=0; col<sampledata->getWidth()-1; col++)
342                {
343                        sampledata->setPoint(row,col); 
344                        colour = sampledata->getHeightShading(shadingres);
345                        ds->heightShading(sampledata->getTopLeft().x,
346                                                          sampledata->getTopLeft().y,
347                                                          sampledata->getTopRight().x,
348                                                          sampledata->getTopRight().y,
349                                                          sampledata->getBottomRight().x,
350                                                          sampledata->getBottomRight().y,
351                                                          sampledata->getBottomLeft().x,
352                                                          sampledata->getBottomLeft().y, 
353                                                          colour.r,colour.g,colour.b);
354                }
355        }
356}
357
Note: See TracBrowser for help on using the repository browser.