Stephen, at first sight, I don't see anything wrong with your code.
You need to include "cpl_conv.h" for CPLMalloc(). And memory allocated with CPLMalloc() should be freed with CPLFree(). The doc is there : http://gdal.org/cpl__conv_8h.html (http://gdal.org/ -> API reference documentation -> Files -> cpl_conv.h ) But that doesn't likely account for the segmentation fault you have. This code should work just fine with a GTiff or VRT dataset. Do you use a recent version of GDAL ? If not, retry ;-) If still crashing, could you attach a debugger and show the stack trace at the place where it crashes and/or run valgrind on it ? Best regards, Even Le dimanche 22 août 2010 05:23:45, Stephen Woodbridge a écrit : > OK, well may to answer my own question by reading the GDAL API Tutorial > and taking a guess at some parts I came up with the following. This is > basically just a copy of the tutorial and at the end I added some code > that I'm hoping reads the correct pixel for my lat,lon and prints out > the elevation. > > The tutorial code for "Fetching a Raster Band" works on a GTiff file but > segv on a VRT file. Not sure why this happened, although there was a > compiler warning about: > > testit.c: In function 'main': > testit.c:80: warning: cast to pointer from integer of different size > > which would be: > > float *pafScanline; > pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize); > > I could not find docs on CPLMalloc, or any info on if I need to free > that and how! > > I would appreciate if someone could code review the code at the end and > comment on whether or not this is correct. > > Thanks, > -Steve > > #include <gdal.h> > > int main() { > GDALDatasetH hDataset; > char *pszFilename = "/u/srcdata/ned2.vrt"; > //char *pszFilename = "/u/srcdata/NED2/-90_28_-84_30.tif"; > > GDALDriverH hDriver; > double adfGeoTransform[6]; > > GDALRasterBandH hBand; > int nBlockXSize, nBlockYSize; > int bGotMin, bGotMax; > double adfMinMax[2]; > > float *pafScanline; > int nXSize; > > int nXOff, nYOff; > float elevation; > > float lon = -74.0; > float lat = 42.0; > > GDALAllRegister(); > > hDataset = GDALOpen( pszFilename, GA_ReadOnly ); > if( hDataset == NULL ) { > printf("Failed to open %s!\n", pszFilename); > exit(1); > } > > hDriver = GDALGetDatasetDriver( hDataset ); > printf( "Driver: %s/%s\n", > GDALGetDriverShortName( hDriver ), > GDALGetDriverLongName( hDriver ) ); > > printf( "Size is %dx%dx%d\n", > GDALGetRasterXSize( hDataset ), > GDALGetRasterYSize( hDataset ), > GDALGetRasterCount( hDataset ) ); > > if( GDALGetProjectionRef( hDataset ) != NULL ) > printf( "Projection is `%s'\n", GDALGetProjectionRef( hDataset ) > ); > > if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None ) > { > printf( "Origin = (%.6f,%.6f)\n", > adfGeoTransform[0], adfGeoTransform[3] ); > > printf( "Pixel Size = (%.6f,%.6f)\n", > adfGeoTransform[1], adfGeoTransform[5] ); > } > > hBand = GDALGetRasterBand( hDataset, 1 ); > GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize ); > printf( "Block=%dx%d Type=%s, ColorInterp=%s\n", > nBlockXSize, nBlockYSize, > GDALGetDataTypeName(GDALGetRasterDataType(hBand)), > GDALGetColorInterpretationName( > GDALGetRasterColorInterpretation(hBand)) ); > > adfMinMax[0] = GDALGetRasterMinimum( hBand, &bGotMin ); > adfMinMax[1] = GDALGetRasterMaximum( hBand, &bGotMax ); > if( ! (bGotMin && bGotMax) ) > GDALComputeRasterMinMax( hBand, TRUE, adfMinMax ); > > printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] ); > > if( GDALGetOverviewCount(hBand) > 0 ) > printf( "Band has %d overviews.\n", GDALGetOverviewCount(hBand)); > > if( GDALGetRasterColorTable( hBand ) != NULL ) > printf( "Band has a color table with %d entries.\n", > GDALGetColorEntryCount( > GDALGetRasterColorTable( hBand ) ) ); > > /* > * This section segv trying to do a memset in GDALRasterIO > * if I use the VRT driver, but works OK on the GTiff file. > */ > /* > nXSize = GDALGetRasterBandXSize( hBand ); > > pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize); > GDALRasterIO( hBand, GF_Read, 0, 0, nXSize, 1, > pafScanline, nXSize, 1, GDT_Float32, > 0, 0 ); > */ > > /* convert the lon, lat in a pixel location on the file and > * read that location into a float. > */ > nXOff = (lon - adfGeoTransform[0]) / adfGeoTransform[1] + 0.5; > nYOff = (lat - adfGeoTransform[3]) / adfGeoTransform[5] + 0.5; > GDALRasterIO( hBand, GF_Read, nXOff, nYOff, 1, 1, > &elevation, 1, 1, GDT_Float32, > 0, 0 ); > printf("Elevation at %.6f, %.6f = %.1f\n", lon, lat, elevation); > > GDALClose(hDataset); > return(0); > } > > On 8/21/2010 6:07 PM, Stephen Woodbridge wrote: > > So I have been able to download the NED data and I have US coverage. > > The data is in GeoTiff and all in a single band Int16 format. > > I have create a ned.vrt file. > > > > So I'm interested in writing a simple program in C that I can pass a > > lat, lon and get the elevation. I assume I need to open the .vrt file > > but looking at the GDAL API, I'm not sure how to map the lat, lon to a > > pixel and then how to fetch the pixel value. If it matter, I will > > probably extend this later to process a list of locations or maybe > > extend it to read a shapefile and process the shape through it. > > > > Does anyone have a sample program that does something like this. > > > > -Steve > > > > > > > > wood...@mappy:/u/srcdata$ gdalinfo -stats ned2.vrt > > Driver: VRT/Virtual Raster > > Files: ned2.vrt > > NED2/-100_25_-96_28.tif > > NED2/-102_28_-96_30.tif > > NED2/-102_30_-96_32.tif > > NED2/-102_32_-96_34.tif > > NED2/-102_34_-96_36.tif > > NED2/-102_36_-96_38.tif > > NED2/-102_38_-96_40.tif > > NED2/-102_40_-96_42.tif > > NED2/-102_42_-96_44.tif > > NED2/-102_44_-96_46.tif > > NED2/-102_46_-96_48.tif > > NED2/-102_48_-96_50.tif > > NED2/-108_28_-102_30.tif > > NED2/-108_30_-102_32.tif > > NED2/-108_32_-102_34.tif > > NED2/-108_34_-102_36.tif > > NED2/-108_36_-102_38.tif > > NED2/-108_38_-102_40.tif > > NED2/-108_40_-102_42.tif > > NED2/-108_42_-102_44.tif > > NED2/-108_44_-102_46.tif > > NED2/-108_46_-102_48.tif > > NED2/-108_48_-102_50.tif > > NED2/-114_30_-108_32.tif > > NED2/-114_32_-108_34.tif > > NED2/-114_34_-108_36.tif > > NED2/-114_36_-108_38.tif > > NED2/-114_38_-108_40.tif > > NED2/-114_40_-108_42.tif > > NED2/-114_42_-108_44.tif > > NED2/-114_44_-108_46.tif > > NED2/-114_46_-108_48.tif > > NED2/-114_48_-108_50.tif > > NED2/-120_32_-114_34.tif > > NED2/-120_34_-114_36.tif > > NED2/-120_36_-114_38.tif > > NED2/-120_38_-114_40.tif > > NED2/-120_40_-114_42.tif > > NED2/-120_42_-114_44.tif > > NED2/-120_44_-114_46.tif > > NED2/-120_46_-114_48.tif > > NED2/-120_48_-114_50.tif > > NED2/-123_33_-120_37.tif > > NED2/-125_37_-120_39.tif > > NED2/-125_39_-120_42.tif > > NED2/-126_42_-120_44.tif > > NED2/-126_44_-120_46.tif > > NED2/-126_46_-120_48.tif > > NED2/-126_48_-120_50.tif > > NED2/133_6.5_138.5_10.tif > > NED2/-139.0019842_51.9997277_-127.0018384_58.9998127.tif > > NED2/-139.0024054_58.9997472_-127.0018101_66.0000944.tif > > NED2/144_13_146_16.tif > > NED2/151_5_163.5_8.tif > > NED2/-157_18_-152_23.tif > > NED2/-157.5_18.75_-154.5_21.5.tif > > NED2/-160.5_21_-157.5_22.5.tif > > NED2/-162_18_-157_23.tif > > NED2/-171_-14.5_-169_-14.tif > > NED2/-68_17_-64_18.tif > > NED2/-72.0002778_40.9999994_-67.9999992_44.tif > > NED2/-72.0002778_44_-66_46.tif > > NED2/-72.0002778_46_-66_48.tif > > NED2/-78.0002778_33.0000047_-75.0000047_36.tif > > NED2/-78.0002778_36_-72_38.tif > > NED2/-78.0002778_38_-72_40.tif > > NED2/-78.0002778_40_-72_42.tif > > NED2/-78.0002778_42_-72_44.tif > > NED2/-78.0002778_44_-72_46.tif > > NED2/-84.0002778_24_-78_26.tif > > NED2/-84.0002778_26_-78_28.tif > > NED2/-84.0002778_28_-78_30.tif > > NED2/-84.0002778_30_-78_32.tif > > NED2/-84.0002778_32_-78_34.tif > > NED2/-84.0002778_34_-78_36.tif > > NED2/-84.0002778_36_-78_38.tif > > NED2/-84.0002778_38_-78_40.tif > > NED2/-84.0002778_40_-78_42.tif > > NED2/-84.0002778_42_-78_44.tif > > NED2/-85.0002778_44.0000047_-82.0000047_47.tif > > NED2/-90_28_-84_30.tif > > NED2/-90_30_-84_32.tif > > NED2/-90_32_-84_34.tif > > NED2/-90_34_-84_36.tif > > NED2/-90_36_-84_38.tif > > NED2/-90_38_-84_40.tif > > NED2/-90_40_-84_42.tif > > NED2/-90_42_-84_44.tif > > NED2/-90_44_-85_46.tif > > NED2/-90_46_-85_49.tif > > NED2/-96_28_-90_30.tif > > NED2/-96_30_-90_32.tif > > NED2/-96_32_-90_34.tif > > NED2/-96_34_-90_36.tif > > NED2/-96_36_-90_38.tif > > NED2/-96_38_-90_40.tif > > NED2/-96_40_-90_42.tif > > NED2/-96_42_-90_44.tif > > NED2/-96_44_-90_46.tif > > NED2/-96_46_-90_48.tif > > NED2/-96_48_-90_50.tif > > Size is 1180816, 284173 > > Coordinate System is: > > GEOGCS["WGS 84", > > DATUM["WGS_1984", > > SPHEROID["WGS 84",6378137,298.257223563, > > AUTHORITY["EPSG","7030"]], > > AUTHORITY["EPSG","6326"]], > > PRIMEM["Greenwich",0], > > UNIT["degree",0.0174532925199433], > > AUTHORITY["EPSG","4326"]] > > Origin = (-171.000000000113090,66.000094447240443) > > Pixel Size = (0.000283278659440,-0.000283278659440) > > Corner Coordinates: > > Upper Left (-171.0000000, 66.0000944) (171d 0'0.00"W, 66d 0'0.34"N) > > Lower Left (-171.0000000, -14.5000520) (171d 0'0.00"W, 14d30'0.19"S) > > Upper Right ( 163.4999735, 66.0000944) (163d29'59.90"E, 66d 0'0.34"N) > > Lower Right ( 163.4999735, -14.5000520) (163d29'59.90"E, 14d30'0.19"S) > > Center ( -3.7500132, 25.7500212) ( 3d45'0.05"W, 25d45'0.08"N) > > Band 1 Block=128x128 Type=Int16, ColorInterp=Gray > > Minimum=-83.000, Maximum=4665.000, Mean=25.333, StdDev=189.129 > > Metadata: > > STATISTICS_MINIMUM=-83 > > STATISTICS_MAXIMUM=4665 > > STATISTICS_MEAN=25.332848474484 > > STATISTICS_STDDEV=189.12886049417 > > _______________________________________________ > > gdal-dev mailing list > > gdal-dev@lists.osgeo.org > > http://lists.osgeo.org/mailman/listinfo/gdal-dev > > _______________________________________________ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
