Dave Seems to me this is a point not a raster process you are looking for
e.g. How do I make a surface from a bunch of discrete points ? So instead of having 3 bands of raster you want a stream of XYZ triplets you can then submit these to gdal_grid or other tools that work with discrete data pseudo python assuming you have read your files into numpy arrays XYZ = array( zip( bandX.flatten(), bandY.flatten(), bandZ.flatten() ) ) convert XYZ to any OGR format call gdal_grid HTH Norman On Jan 22, 2013, at 10:12 PM, David Hoese <dho...@gmail.com> wrote: > 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 <dho...@gmail.com >> <mailto:dho...@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 <dho...@gmail.com >>> <mailto:dho...@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@lists.osgeo.org <mailto: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
