root/OpenSceneGraph/trunk/src/osgPlugins/gdal/DataSetLayer.cpp @ 13041

Revision 13041, 7.0 kB (checked in by robert, 3 years ago)

Ran script to remove trailing spaces and tabs

  • Property svn:eol-style set to native
Line 
1/* -*-c++-*- OpenSceneGraph - Copyright (C) 1998-2007 Robert Osfield
2 *
3 * This library is open source and may be redistributed and/or modified under
4 * the terms of the OpenSceneGraph Public License (OSGPL) version 0.0 or
5 * (at your option) any later version.  The full license is in LICENSE file
6 * included with this distribution, and on the openscenegraph.org website.
7 *
8 * This library is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11 * OpenSceneGraph Public License for more details.
12*/
13
14#include "DataSetLayer.h"
15
16#include <osg/Notify>
17#include <osg/io_utils>
18
19#include <cpl_string.h>
20#include <gdalwarper.h>
21#include <ogr_spatialref.h>
22
23#include <osgDB/ImageOptions>
24
25using namespace GDALPlugin;
26
27DataSetLayer::DataSetLayer():
28_dataset(0), _gdalReader(0)
29{
30}
31
32DataSetLayer::DataSetLayer(const std::string& fileName):
33_dataset(0), _gdalReader(0)
34{
35    setFileName(fileName);
36    open();
37}
38
39DataSetLayer::DataSetLayer(const DataSetLayer& dataSetLayer,const osg::CopyOp& copyop):
40    Layer(dataSetLayer), _gdalReader(dataSetLayer._gdalReader)
41{
42    if (dataSetLayer._dataset) open();
43}
44
45DataSetLayer::~DataSetLayer()
46{
47    close();
48}
49
50void DataSetLayer::open()
51{
52    if (_dataset) return;
53
54    if (getFileName().empty()) return;
55
56
57    OSG_NOTICE<<"DataSetLayer::open()"<<getFileName()<<std::endl;
58
59    _dataset = static_cast<GDALDataset*>(GDALOpen(getFileName().c_str(),GA_ReadOnly));
60
61    setUpLocator();
62}
63
64void DataSetLayer::close()
65{
66    OSG_NOTICE<<"DataSetLayer::close()"<<getFileName()<<std::endl;
67
68    if (_dataset)
69    {
70        GDALClose(static_cast<GDALDatasetH>(_dataset));
71
72        _dataset = 0;
73    }
74}
75
76unsigned int DataSetLayer::getNumColumns() const
77{
78    return _dataset!=0 ? _dataset->GetRasterXSize() : 0;
79}
80
81unsigned int DataSetLayer::getNumRows() const
82{
83    return _dataset!=0 ? _dataset->GetRasterYSize() : 0;
84}
85
86osgTerrain::ImageLayer* DataSetLayer::extractImageLayer(unsigned int sourceMinX, unsigned int sourceMinY, unsigned int sourceMaxX, unsigned int sourceMaxY, unsigned int targetWidth, unsigned int targetHeight)
87{
88    if (!_dataset || sourceMaxX<sourceMinX || sourceMaxY<sourceMinY) return 0;
89
90    if (!_gdalReader) return 0;
91
92    osg::ref_ptr<osgDB::ImageOptions> imageOptions = new osgDB::ImageOptions;
93    imageOptions->_sourceImageWindowMode = osgDB::ImageOptions::PIXEL_WINDOW;
94    imageOptions->_sourcePixelWindow.windowX = sourceMinX;
95    imageOptions->_sourcePixelWindow.windowY = sourceMinY;
96    imageOptions->_sourcePixelWindow.windowWidth = sourceMaxX-sourceMinX;
97    imageOptions->_sourcePixelWindow.windowHeight = sourceMaxY-sourceMinY;
98    imageOptions->_destinationPixelWindow.windowX = 0;
99    imageOptions->_destinationPixelWindow.windowY = 0;
100    imageOptions->_destinationPixelWindow.windowWidth = targetWidth;
101    imageOptions->_destinationPixelWindow.windowHeight = targetHeight;
102
103    osgDB::ReaderWriter::ReadResult result = _gdalReader->readImage(getFileName(),imageOptions.get());
104    osg::Image* image = result.getImage();
105    if (!image) return 0;
106
107    osgTerrain::ImageLayer* layer = new osgTerrain::ImageLayer;
108    layer->setFileName(getFileName());
109    layer->setImage(image);
110
111    return layer;
112}
113
114void DataSetLayer::setGdalReader(const osgDB::ReaderWriter* rw)
115{
116    _gdalReader = const_cast<osgDB::ReaderWriter*>(rw);
117}
118
119void DataSetLayer::setUpLocator()
120{
121    if (!isOpen()) return;
122
123    const char* pszSourceSRS = _dataset->GetProjectionRef();
124    if (!pszSourceSRS || strlen(pszSourceSRS)==0) pszSourceSRS = _dataset->GetGCPProjection();
125
126
127    osg::ref_ptr<osgTerrain::Locator> locator = new osgTerrain::Locator;
128
129    if (pszSourceSRS)
130    {
131        locator->setFormat("WKT");
132        locator->setCoordinateSystem(pszSourceSRS);
133    }
134
135    osg::Matrixd matrix;
136
137    double geoTransform[6];
138    if (_dataset->GetGeoTransform(geoTransform)==CE_None)
139    {
140#ifdef SHIFT_RASTER_BY_HALF_CELL
141        // shift the transform to the middle of the cell if a raster interpreted as vector
142        if (data->_dataType == VECTOR)
143        {
144            geoTransform[0] += 0.5 * geoTransform[1];
145            geoTransform[3] += 0.5 * geoTransform[5];
146        }
147#endif
148
149        matrix.set( geoTransform[1],    geoTransform[4],    0.0,    0.0,
150                    geoTransform[2],    geoTransform[5],    0.0,    0.0,
151                    0.0,                0.0,                1.0,    0.0,
152                    geoTransform[0],    geoTransform[3],    0.0,    1.0);
153
154
155        int nPixels = _dataset->GetRasterXSize();
156        int nLines = _dataset->GetRasterYSize();
157
158        locator->setTransform(
159            osg::Matrixd::scale(static_cast<double>(nPixels-1), static_cast<double>(nLines-1), 1.0) *
160            matrix);
161
162        locator->setDefinedInFile(true);
163
164        setLocator(locator.get());
165
166    }
167    else if (_dataset->GetGCPCount()>0 && _dataset->GetGCPProjection())
168    {
169        OSG_NOTICE << "    Using GCP's"<< std::endl;
170
171
172        /* -------------------------------------------------------------------- */
173        /*      Create a transformation object from the source to               */
174        /*      destination coordinate system.                                  */
175        /* -------------------------------------------------------------------- */
176        void *hTransformArg =
177            GDALCreateGenImgProjTransformer( _dataset, pszSourceSRS,
178                                             NULL, pszSourceSRS,
179                                             TRUE, 0.0, 1 );
180
181        if ( hTransformArg == NULL )
182        {
183            OSG_NOTICE<<" failed to create transformer"<<std::endl;
184            return;
185        }
186
187        /* -------------------------------------------------------------------- */
188        /*      Get approximate output definition.                              */
189        /* -------------------------------------------------------------------- */
190        double adfDstGeoTransform[6];
191        int nPixels=0, nLines=0;
192        if( GDALSuggestedWarpOutput( _dataset,
193                                     GDALGenImgProjTransform, hTransformArg,
194                                     adfDstGeoTransform, &nPixels, &nLines )
195            != CE_None )
196        {
197            OSG_NOTICE<<" failed to create warp"<<std::endl;
198            return;
199        }
200
201        GDALDestroyGenImgProjTransformer( hTransformArg );
202
203
204        matrix.set( adfDstGeoTransform[1],    adfDstGeoTransform[4],    0.0,    0.0,
205                    adfDstGeoTransform[2],    adfDstGeoTransform[5],    0.0,    0.0,
206                    0.0,                0.0,                1.0,    0.0,
207                    adfDstGeoTransform[0],    adfDstGeoTransform[3],    0.0,    1.0);
208
209        locator->setTransform(
210            osg::Matrixd::scale(static_cast<double>(nPixels-1), static_cast<double>(nLines-1), 1.0) *
211            matrix);
212
213        locator->setDefinedInFile(true);
214
215        setLocator(locator.get());
216    }
217    else
218    {
219        OSG_INFO << "DataSetLayer::setUpLocator(), No GeoTransform or GCP's - unable to compute position in space"<< std::endl;
220    }
221
222}
Note: See TracBrowser for help on using the browser.