Revision: 70123
http://sourceforge.net/p/brlcad/code/70123
Author: starseeker
Date: 2017-08-24 16:11:49 +0000 (Thu, 24 Aug 2017)
Log Message:
-----------
We'll need to be able to project data... start working on figuring out how to
do that.
Modified Paths:
--------------
brlcad/trunk/src/libgcv/plugins/gdal/gdal.cpp
Modified: brlcad/trunk/src/libgcv/plugins/gdal/gdal.cpp
===================================================================
--- brlcad/trunk/src/libgcv/plugins/gdal/gdal.cpp 2017-08-24 14:18:25 UTC
(rev 70122)
+++ brlcad/trunk/src/libgcv/plugins/gdal/gdal.cpp 2017-08-24 16:11:49 UTC
(rev 70123)
@@ -43,6 +43,7 @@
/* GDAL headers */
#include <gdal.h>
+#include <gdalwarper.h>
#include <gdal_utils.h>
#include <cpl_conv.h>
#include <cpl_string.h>
@@ -131,7 +132,111 @@
return 1;
}
+/* Pull the latitude and longitude out of the dataset. Based on gdalinfo
code, so
+ * this function (gdal_ll) is:
+ *
+ * Copyright (c) 1998, Frank Warmerdam
+ * Copyright (c) 2007-2015, Even Rouault <even.rouault at spatialys.com>
+ * Copyright (c) 2015, Faza Mahamood
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ */
HIDDEN int
+gdal_ll(GDALDatasetH hDataset, double *lat_center, double *long_center)
+{
+ double longitude, latitude;
+ double adfGeoTransform[6];
+ const char *pzp = NULL;
+ OGRCoordinateTransformationH htfm = NULL;
+ double gx = GDALGetRasterXSize(hDataset)/2.0;
+ double gy = GDALGetRasterYSize(hDataset)/2.0;
+ if(GDALGetGeoTransform(hDataset, adfGeoTransform) == CE_None ) {
+ longitude = adfGeoTransform[0] + adfGeoTransform[1] * gx +
adfGeoTransform[2] * gy;
+ latitude = adfGeoTransform[3] + adfGeoTransform[4] * gx +
adfGeoTransform[5] * gy;
+ pzp = GDALGetProjectionRef(hDataset);
+ } else {
+ return 0;
+ }
+
+ if(pzp != NULL && strlen(pzp) > 0) {
+ OGRSpatialReferenceH hpj, hll = NULL;
+ hpj = OSRNewSpatialReference( pzp );
+ if(hpj != NULL) hll = OSRCloneGeogCS(hpj);
+ if(hll != NULL) {
+ CPLPushErrorHandler(CPLQuietErrorHandler);
+ /* The bit we need... */
+ htfm = OCTNewCoordinateTransformation(hpj, hll);
+ CPLPopErrorHandler();
+ OSRDestroySpatialReference(hll);
+ }
+ if(hpj != NULL) OSRDestroySpatialReference(hpj);
+ }
+
+ if(htfm != NULL) {
+ if (OCTTransform(htfm,1,&longitude,&latitude,NULL)) {
+ if (lat_center) (*lat_center) = latitude;
+ if (long_center) (*long_center) = longitude;
+ } else {
+ return 0;
+ }
+ } else {
+ if (lat_center) (*lat_center) = latitude;
+ if (long_center) (*long_center) = longitude;
+ }
+ bu_log("lat: %f, long: %f\n", latitude, longitude);
+
+ return 1;
+}
+
+
+/* Get the UTM zone and corresponding EPSG number - see
+ *
https://gis.stackexchange.com/questions/241696/how-to-convert-from-lat-lon-to-utm-with-gdaltransform
+ * and the linked posts in the answers for more info. */
+HIDDEN int
+gdal_utm_zone(struct conversion_state *state)
+{
+ int zone = INT_MAX;
+ double longitude = 0.0;
+ if (gdal_ll(state->hDataset, NULL, &longitude)) {
+ zone = floor((longitude + 180) / 6) + 1;
+ }
+ bu_log("zone: %d\n", zone);
+ return zone;
+}
+
+
+/* Get corresponding EPSG number of the UTM zone - see
+ *
https://gis.stackexchange.com/questions/241696/how-to-convert-from-lat-lon-to-utm-with-gdaltransform
+ * and the linked posts in the answers for more info. */
+HIDDEN int
+gdal_utm_epsg(struct conversion_state *state, int zone)
+{
+ int epsg = INT_MAX;
+ double latitude = 0.0;
+ if (gdal_ll(state->hDataset, &latitude, NULL)) {
+ epsg = (latitude > 0) ? 32600 + zone : 32700 + zone;
+ }
+ bu_log("epsg: %d\n", epsg);
+ return epsg;
+}
+
+HIDDEN int
gdal_read(struct gcv_context *context, const struct gcv_opts *gcv_options,
const void *options_data, const char *source_path)
{
@@ -154,7 +259,51 @@
(void)get_dataset_info(state);
- /* Read in the data */
+ int zone = gdal_utm_zone(state);
+ (void)gdal_utm_epsg(state, zone);
+
+#if 0
+ GDALDataType eDT =
GDALGetRasterDataType(GDALGetRasterBand(state->hDataset,1));
+ GDALDriverH hD = GDALGetDriverByName("MEM");
+ const char *src_Wkt = GDALGetProjectionRef(state->HDataset);
+ OGRSpatialReference oSRS;
+ oSRS.SetUTM(gdal_utm_zone(state));
+ oSRS.SetWellKnownGeogCS();
+ const char *dst_Wkt;
+ oSRS.exportToWkt(&dst_Wkt);
+ void *hTsfA = GDALCreateGenImgProjTransformer(state->hDataset, src_Wkt,
NULL dst_Wkt, FALSE, 0, 1);
+ double adfDstGeoTransform[6];
+ int nPixels=0, nLines=0;
+ CPLErr eErr = GDALSuggestedWarpOutput(state->hDataset,
GDALGenImgProjTransform, hTsfA, adfDstGeoTransform, &nPixels, &nLines);
+ GDALDestroyGenImgProjTransformer(hTsfA);
+ if (eErr != CE_None) return 0;
+
+ GDALDatasetH pjdata = GDALCreate(hD, NULL, nPixels, nLines,
GDALGetRasterCount(state->hDataset), eDT, NULL);
+ if (pjdata == NULL) return 0;
+
+ GDALSetProjection(pjdata, dst_Wkt);
+ GDALSetGeoTransform(pjdata, adfDstGeoTransform);
+ GDALWarpOptions *w_opts = GDALCreateWarpOptions();
+ w_opts->hSrcDS = state->hDataset;
+ w_opts->hDstDS = pjdata;
+ w_opts->nBandCount = 1;
+ w_opts->panSrcBands = (int *)CPLMalloc(sizeof(int) * w_opts->nBandCount);
+ w_opts->panSrcBands[0] = 1;
+ w_opts->panDstBands = (int *)CPLMalloc(sizeof(int) * w_opts->nBandCount);
+ w_opts->panDstBands[0] = 1;
+ w_opts->pfnProgress = GDALTermProgress;
+ // Establish reprojection transformer.
+ w_opts->pTransformerArg = GDALCreateGenImgProjTransformer(state->hDataset,
GDALGetProjectionRef(state->hDataset), pjdata, GDALGetProjectionRef(pjdata),
FALSE, 0.0, 1 );
+ w_opts->pfnTransformer = GDALGenImgProjTransform;
+ // Initialize and execute the warp operation.
+ GDALWarpOperation oOperation;
+ oOperation.Initialize(w_opts);
+ oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(pjdata),
GDALGetRasterYSize(pjdata));
+ GDALDestroyGenImgProjTransformer(w_opts->pTransformerArg);
+ GDALDestroyWarpOptions(w_opts);
+#endif
+
+ /* Read the data into something a DSP can process */
unsigned int xsize = GDALGetRasterXSize(state->hDataset);
unsigned int ysize = GDALGetRasterYSize(state->hDataset);
unsigned short *uint16_array = (unsigned short *)bu_calloc(xsize*ysize,
sizeof(unsigned short), "unsigned short array");
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
BRL-CAD Source Commits mailing list
brlcad-commits@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/brlcad-commits
