[gdal-dev] Working with NED elevation data - part 2

Brent Fraser bfraser at geoanalytic.com
Mon Aug 23 10:17:56 EDT 2010


Steve,

   Excellent info!  I had been contemplating using GDAL's VRT format to create 
an elevation web service so I could display X,Y and Z on mouse moves in web 
mapping, but was unsure about the performance.  Looks like fcgi could work.

   FYI, I think I'll use Google's <ElevationResponse> XML format 
(http://code.google.com/apis/maps/documentation/elevation/), but I may extend to 
have some kind of quality/accuracy and vertical datum values too.  Here's a 
stock Google example:

<ElevationResponse>
  <status>OK</status>
  <result>
   <location>
    <lat>39.7391536</lat>
    <lng>-104.9847034</lng>
   </location>
   <elevation>1608.8402100</elevation>
  </result>
  <result>
   <location>
    <lat>36.4555560</lat>
    <lng>-116.8666670</lng>
   </location>
   <elevation>-50.7890358</elevation>
  </result>
</ElevationResponse>


Or perhaps the less verbose USGS format 
(http://gisdata.usgs.gov/XMLWebServices2/Elevation_Service_Methods.php):

<Data_ID>[Source Layer ID]</Data_ID>
<Elevation>[Decimal Elevation Value]</Elevation>
<Units>[FEET or METERS]</Units>

Thanks again...

Best Regards,
Brent Fraser

Stephen Woodbridge wrote:
> Even,
> 
> It seems the adding:
> 
> #include <cpl_conv.h>
> 
> resolves the segv issue.
> 
> I would recommend changing the example in:
> 
> http://www.gdal.org/gdal_tutorial.html
> 
> To include that in the section "Reading Raster Data" because the example 
> C code uses CPLMalloc(). It would also be very helpful to add:
> 
> ...
> CPLFree(pafScanline);
> 
> to the example code. And a link to http://gdal.org/cpl__conv_8h.html 
> which is not otherwise referenced.
> 
> Anyway, code is working great now. I am generating an elevation profile 
> from a number of points along a path.
> 
> 1. tried to get the points via USGS NED webservice: 132 point took 4-6 
> min and 132 points is a short path.
> 
> So I downloaded the USGS NED data and created a .vrt file for the layer:
> 
> 2. C code to open vrt file and handle requests via local CGI took 4 sec
> 3. C code to open vrt file and handle requests via local FCGI took < 1 sec
> 
> And each profile time above has some overhead because it has to compute 
> a route in pgrouting to get the points needed to generate the profile.
> 
> http://imaptools.com:8080/routeloops/route_loops6.php?x=-71.387764&y=42.620002&length=4828.032&heading=-1&shape=-1&rand=0.06593320066404917&allow=0%2C1&zoom=11&lat=42.62&lon=-71.38776&layers=B00TT&ll=-71.387764%2042.620002&len=3&rs=0.4646901273751023&al=0,1&u=E&cw=1&f=profile 
> 
> 
> Thanks for the help and support everyone.
> 
> -Steve
> 
> On 8/22/2010 12:43 PM, Even Rouault wrote:
>> 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
>>>>
>>>>
>>>>
>>>> woodbri at 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 at lists.osgeo.org
>>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>>
>>> _______________________________________________
>>> gdal-dev mailing list
>>> gdal-dev at lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
> 



More information about the gdal-dev mailing list