[gdal-dev] reproject python numpy binary swath/lat/lon

[Etienne Tourigny]

On Thu, Jun 14, 2012 at 4:15 AM, Rutger <kassies at gmail.com> wrote:
> Hey David,
>
>
> David Hoese wrote
>>
>> On 6/13/12 10:00 AM, gdal-dev-request at .osgeo wrote
>> I'm not sure at the moment, it's VIIRS data.  I'm playing around with
>> different existing remapping tools, I have one that works called
>> ms2gt(ll2cr/fornav) which you might be able to google for the old
>> versions, but it's slightly outdated and is single purpose.  I would
>> like to use GDAL if it can handle it, because of the various formats
>> that it supports as output.
>>
>
> I actually already did some test with GDAL geolocation arrays and VIIRS
> data. It work flawless for my target area; utm31 for the Netherlands. My
> VRT's looked like this:
>
> lat.vrt
> <VRTDataset rasterXSize="6400" rasterYSize="1536">
>  <VRTRasterBand dataType="Float32" band="1">
>    <Metadata />
>    <SimpleSource>
>      <SourceFilename
> relativeToVRT="0">HDF5:GITCO_npp_d20120315_t1208597_e1210238_b01974_c20120321152146856862_noaa_ops.h5://All_Data/VIIRS-IMG-GEO-TC_All/Latitude</SourceFilename>
>      <SourceBand>1</SourceBand>
>      <SourceProperties RasterXSize="6400" RasterYSize="1536"
> DataType="Float32" BlockXSize="6400" BlockYSize="1536" />
>      <SrcRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
>      <DstRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
>    </SimpleSource>
>  </VRTRasterBand>
> </VRTDataset>
>
> lon.vrt
> <VRTDataset rasterXSize="6400" rasterYSize="1536">
>  <VRTRasterBand dataType="Float32" band="1">
>    <Metadata />
>    <SimpleSource>
>      <SourceFilename
> relativeToVRT="0">HDF5:GITCO_npp_d20120315_t1208597_e1210238_b01974_c20120321152146856862_noaa_ops.h5://All_Data/VIIRS-IMG-GEO-TC_All/Longitude</SourceFilename>
>      <SourceBand>1</SourceBand>
>      <SourceProperties RasterXSize="6400" RasterYSize="1536"
> DataType="Float32" BlockXSize="6400" BlockYSize="1536" />
>      <SrcRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
>      <DstRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
>    </SimpleSource>
>  </VRTRasterBand>
> </VRTDataset>
>
> band5.vrt
> <VRTDataset rasterXSize="6400" rasterYSize="1536">
>   <Metadata domain="GEOLOCATION">
>     <MDI key="X_DATASET">lon.vrt</MDI>
>     <MDI key="X_BAND">1</MDI>
>     <MDI key="Y_DATASET">lat.vrt</MDI>
>     <MDI key="Y_BAND">1</MDI>
>     <MDI key="PIXEL_OFFSET">0</MDI>
>     <MDI key="LINE_OFFSET">0</MDI>
>     <MDI key="PIXEL_STEP">1</MDI>
>     <MDI key="LINE_STEP">1</MDI>
>   </Metadata>
>  <VRTRasterBand dataType="UInt16" band="1">
>    <Metadata />
>    <SimpleSource>
>      <SourceFilename
> relativeToVRT="0">HDF5:SVI05_npp_d20120315_t1208597_e1210238_b01974_c20120316203244033373_noaa_ops.h5://All_Data/VIIRS-I5-SDR_All/BrightnessTemperature</SourceFilename>
>      <SourceBand>1</SourceBand>
>      <SourceProperties RasterXSize="6400" RasterYSize="1536"
> DataType="UInt16" BlockXSize="6400" BlockYSize="1536" />
>      <SrcRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
>      <DstRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
>    </SimpleSource>
>  </VRTRasterBand>
> </VRTDataset>
>
> And then reprojection with gdalwarp like:
> gdalwarp -s_srs "epsg:4326" -t_srs "+proj=utm +zone=31 +datum=WGS84
> +units=m" \
> -tr 250 250 -te 392875 5574125 808875 5940125 band5.vrt band5.tif
>
> This works without any errors for me, and the geolocation seems accurate.
> But as said, i didnt like the resulting image quality. I have attached a
> simple comparison between a gdalwarp nearest neighbor, bilinear and
> Pyresample gaussian reprojection. The gdalwarp reprojection seem very
> patchy, and the bilinear flag (same with cubic) doesnt seem to have much
> affect. Im not sure whether this is to be expected with the current
> implementation or if i made any mistakes while specifying the VRT files.

Did you try the other resampling  methods - cubicspline and lanczos ?

It seems that gaussian resampling is only available for gdaladdo.

>
>
> David Hoese wrote
>>
>> So I've played around with this for a while now and have gotten
>> somewhere.  I was able to get data into the geotiff, but with a lot of
>> guess work and small tweaking.  My t_srs is now:
>> "+a=6378137 +lon_1=-95.000 +b=6356752.314 +lon_0=-123.044 +proj=lcc
>> +lat_1=25.000 +lat_0=59.844 +ellps=WGS84 +no_defs +wktext +units=m"
>>
>> All of which are probably not necessary, but at this point it doesn't
>> hurt.  The origin in that projection is the upper left corner of my
>> grid.  So I now specify -te as something like:
>> "-te 0.0 -5120000.0 5120000.0 0.0"
>>
>> Which gets me something in the geotiff.  I've played around with
>> different extents and my best result comes from me manually changing the
>> extents to match a picture I have from another system.  The grid I'm
>> trying to match to is the first one on this page:
>> http://www.nws.noaa.gov/noaaport/html/icdtb48_2.html
>>
>> I have tried using pyproj and plugging in different corner points to get
>> the extent points in meters, but when I use those I get skewed images
>> (the grid should be 5120 x 5120).
>>
>> So at this point, I just need to find a way to get the extents figured
>> out in a programmatic way (like 1km * 5120 grid pixels = 5120000
>> meters).  If you have any more advice that'd be great, otherwise I'll
>> just have to keep playing around and asking people around my work.
>> Thanks a lot for the help so far.
>>
>> -Dave
>>
>
> Since it worked for my case, i think its likely that the problems you are
> having are related to the projection and area definition. I have no
> experience with the AWIPS grid you are mentioning. What i would do in such a
> case is to get a sample file with the correct projection and boundingbox,
> and just read the metadata with:
> gdalinfo -proj4 samplefile.tif
>
> Optionally add some of the -norat -noct etc flags to hide the large tables.
> The proj4 string can be directly used by other gdal tools.
>
>
> Regards,
> Rutger
>
> --
> View this message in context: http://osgeo-org.1560.n6.nabble.com/gdal-dev-reproject-python-numpy-binary-swath-lat-lon-tp4978609p4981318.html
> Sent from the GDAL - Dev mailing list archive at Nabble.com.
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev