[gdal-dev] gridding binary data

[David Hoese]

Jan, please CC the mailing list in your replies.

I have asked questions about not regularly spaced data on this mailing 
list before, I just can't find the thread atm. In your first point, do 
you mean my input data is "effectively a projected surface"?  If so, 
then no it is not projected, there are lat/lon values for each 
individual pixel and they aren't evenly spaced.  If you mean the output 
is projected, then yes, my software does the projecting from lat/lon 
space into a projection's metered (or degree) space. The data can then 
be interpolated and treated like gridding a scatter plot (see the 
griddata function in scipy and matlab).

The problem with GDAL and data that is not equally spaced geographically 
is that it describes the data with a geotransform.  The geotransform 
(http://www.gdal.org/gdal_datamodel.html) only knows the origin of the 
data and the pixel size (ignoring the 2 rotation elements). So for my 
satellite data (VIIRS and MODIS), the data can not be properly 
represented. My data starts out at one resolution near the center of 
each row and loses quality as you go outward due to...the earth's shape :).

So if I understood you correctly, yes it is a hassle because you have to 
look at each lat/lon point. If you don't look at each point (as I have 
done with gdalwarp and its geotranforms) the produced image isn't 
"correct".  The quality is wrong and you have problems with "the bow-tie 
effect".  As far as I know, there are only 2 pieces of software 
currently in development that handle satellite data in this way. One is 
my software which isn't fully developed and uses a modified version of 
some old software to do the remapping (which is what I'm trying to 
replace, but it works). The second is pytroll, but doesn't meet the 
needs of my bosses and is slower than the old software I'm using.

I hope at least some of this makes sense. Thanks for sticking with my 
attempt at explanation. My final goal would be to add to GDAL the 
ability or a utility to remap satellite data. Sorry this got so long, I 
wasn't sure how much you knew about GDAL or satellite observations.

-Dave

On 1/22/2013 4:34 PM, Jan Heckman wrote:
> Hi Dave,
> It sounds as you have effectively a projected surface. The projected
> (satellite) image consists of regularly spaced pixels? So that the
> pixels would ideally be warped back, perhaps not image-wise but
> pixel-wise, which would be a bit of a hassle...  Don't see yet that the
> warp can't be done in one go, but that shows I'm still missing stuff, or
> that the transformation has to be reformulated to get it to work.
> So I have some hope that fiddling with gdalwarp might produce something
> useful. Then, the sheer amount of data may ultimately be best served by
> a dedicated routine using parallel code. You seem to have figured out
> the geometry pretty much already.
> Maybe I'm not getting it yet?
> Jan
>
> On Tue, Jan 22, 2013 at 10:52 PM, David Hoese <dhoese at gmail.com
> <mailto:dhoese at gmail.com>> wrote:
>
>     Hey Jan,
>
>     Thanks for the response. My process will need to be automated so
>     excel won't help.  I also will be working with potentially large
>     amounts of data (200MB - ~4GB).  I would like the output to be in a
>     flat binary format (just the data), but if I can pull the data out
>     of another format that will do for now.  The output data will be
>     used in other software that I have written to place the data in any
>     format I request (geotiff, AWIPS compatible netCDF files, etc).  I
>     was hoping to use gdal_grid for its algorithms to do the
>     interpolation from input swath pixels to output grid pixels.  I am
>     working with "raw" satellite observation data that is not
>     equidistant (so it can't be put through gdalwarp which expects an
>     equidistant grid as input).
>
>     The x and y arrays are created by software I have written.  The
>     values in these arrays state the x and y offset from the grid origin
>     of where the image pixel is.  For example, a single pixel might be
>     at column 1.14 and row 1.6 of the grid I am mapping to. I would like
>     gdal_grid (or some other code) to use an algorithm (nearest
>     neighbor, averaging, etc) to interpolate this image pixel value onto
>     the surrounding grid points.  Does that make sense? It's difficult
>     for me to describe since I've been working with it so much lately.
>
>     The start-to-finish goal is essentially what gdalwarp does when
>     remapping one grid to another, but this doesn't make the assumption
>     that the data points are equidistant.  The gridding utility I'm
>     currently using is old and only offers a single algorithm for
>     interpolation.  That algorithm does not produce "pretty" images for
>     a specific instrument's data so I'm looking for alternatives without
>     having to hack up the old utility.  Maybe looking at it this way
>     will make more sense:
>
>     (lat,lon,projection) -> my software -> (y,x)
>     (x,y,image_in) -> gdal_grid -> (gridded_image_out)
>
>     Maybe I'll just write the python wrapper for GDALGridCreate since
>     the rest of my software is python.  As far as I know GDAL doesn't
>     have anything else that does gridding and doesn't have anything that
>     can do non-equidistant data points. Thanks for the help so far. Man
>     this email got long fast.
>
>     -Dave
>
>
>     On 1/22/13 2:58 PM, Jan Heckman wrote:
>>     Hi Dave,
>>     So you have effectively three arrays, one for pixel values, one
>>     for x, one for y (probably offsets), which are coordinated and you
>>     want to generate a raster or grid which you can display / use in a
>>     GIS?
>>     Sort of sparse array to be converted to a full array, effectively?
>>     It sounds crazy, maybe, but consider converting these to separated
>>     text (numbers), bring them into excel (columns), tell (arc)gis to
>>     convert to a (point) shapefile and finally, if necessary, convert
>>     that to a raster.
>>     If this is the job, it can alternatively be done quite simply by
>>     writing a small procedure (C++ or python (or something else, I'm
>>     sure)), without using excel etc, as long as you know the output
>>     format (and maybe have a geoprocessor to call upon). A simple
>>     output format is .bil. This is roughly simlar to shapefiles, with
>>     some supporting files to use it in a GIS. Possibly flat 32 bit
>>     float (.flt) is understood by gdal as well.
>>     Sorry that I can't help you with the gdal/ogr commands.
>>     If you don't get a helpful gdal command, give me a yell with an
>>     example.
>>     Jan
>>
>>     On Sun, Jan 20, 2013 at 9:51 PM, David Hoese <dhoese at gmail.com
>>     <mailto:dhoese at gmail.com>> wrote:
>>
>>         I have 3 flat binary files that I would like to use to grid
>>         data. The 3 binary files are the "scattered" image data pixels
>>         and the other two specify the fractional column and row for
>>         each input pixel in the output grid.  So the image pixel in
>>         the first file at offset 0 is mapped to the grid column at
>>         offset 0 (from the columns file) and the grid row at offset 0
>>         (from the rows file).
>>
>>         One way I think I can get the data gridded is to use the
>>         "gdal_grid" utility, but I'm not really sure how to get the
>>         flat binary files to be accepted by "gdal_grid".  I thought I
>>         could use VRT files, but I'm not exactly sure how I would do
>>         that.  I wasn't sure since the documentation here
>>         http://www.gdal.org/ogr/drv_vrt.html says that I need to
>>         specify a "SrcLayer", but binary files don't have layers. So
>>         then would I need an intermediate VRT file to define layers
>>         for the binaries?  If I need an intermediate VRT maybe I
>>         should just push the data into a geotiff (with no projection
>>         info) and try using that?  Could someone clear this up for me?
>>
>>         My other question is has there been any work done to make a
>>         python wrapper for GDALGridCreate().  I found email threads
>>         and a ticket (http://trac.osgeo.org/gdal/ticket/2023) talking
>>         about creating one, but then conversation just stops and those
>>         are from 5 years ago.  I didn't notice anything that screamed
>>         "this is impossible" when skimming through those threads and
>>         the ticket.  I'm also sure a lot has changed for GDAL and
>>         numpy in the last 5 years.  If it's a time thing I'm willing
>>         to add it if someone can point me in the right direction to
>>         start (how-to on adding to GDAL).
>>
>>         Thanks.
>>
>>         -Dave
>>         _______________________________________________
>>         gdal-dev mailing list
>>         gdal-dev at lists.osgeo.org <mailto:gdal-dev at lists.osgeo.org>
>>         http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>>
>
>