Unable to warp HDF5 files

Question

I have numerous HDF5 files(.nc) which I need to batch process warp using one of the gdal utilities gdalwarp. When I tried to warp the files an error occurred:

INPUT:

gdalwarp -geoloc -te 109.975 3.475 135.025 25.025 HDF5:"@file"://geophysical_data/chlor_a %out_path%\@fname.tif"

RESULT:

ERROR 1: Unable to compute a GEOLOC_ARRAY based transformation between pixel/lin e and georeferenced coordinates for HDF5:A2015045060500.L2_LAC_OC.nc://geophysical _data/chlor_a.

Update1:

Just to make it clear, do you mean in lat.vrt, lon.vrt and chlor.vrt I should remove the GCP Id's and MDI key and insert this section:

   <metadata domain="GEOLOCATION">  
     <mdi key="X_DATASET">oc-long.vrt</mdi>  
     <mdi key="X_BAND">1</mdi>  
     <mdi key="Y_DATASET">oc-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>

in between this section?

<VRTDataset rasterXSize="1116" rasterYSize="1610">

  ###### metadata section here #######

 <VRTRasterBand dataType="Float32" band="1">
    <SimpleSource>
      <SourceFilename relativeToVRT="0">HDF5:A2015194044000.L2_LAC.SeAHABS.nc://geophysical_data/chlor_a</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1116" RasterYSize="1610" DataType="Float32" BlockXSize="373" BlockYSize="17" />
      <SrcRect xOff="0" yOff="0" xSize="1116" ySize="1610" />
      <DstRect xOff="0" yOff="0" xSize="1116" ySize="1610" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

Does it work if you separate -geoloc and -te in two commands? — AndreJ, Jul 14 '15 at 09:24
Do you mean run gdalwarp with -geoloc first then run with -te again? Haven't tried it yet. — lovelyvm, Jul 14 '15 at 09:41
From my answer here http://gis.stackexchange.com/questions/128040/warping-images-with-unknown-projection/128068#128068 I guess that gdalwarp -geoloc -of GTIFF -t_srs EPSG:4326 -te ... should run. But you have to use the right values for te, taken from the metadata as reported by gdalinfo. — AndreJ, Jul 14 '15 at 10:23
I have tried this with HDF4 files and it worked but when the format became HDF5 it showed that error above. — lovelyvm, Jul 14 '15 at 11:35
The <metadata> should be removed too, or closed before the new section. The lon.vrt and lat.vrt get no GEOLOCATION, only the chlor.vrt. — AndreJ, Jul 20 '15 at 05:19

AndreJ · Accepted Answer · 2019-04-23T06:09:31.107

After some testing, I think that the geoloc is not working properly. So I used the alternative method using manually created vrt files:

Create a file named lon.vrt:

<VRTDataset rasterXSize="1354" rasterYSize="2030">
  <SRS>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]</SRS>
  <VRTRasterBand dataType="Float32" band="1">
    <SimpleSource>
      <SourceFilename relativeToVRT="1">HDF5:A2015045060000.L2_LAC_OC.nc://navigation_data/longitude</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1354" RasterYSize="2030" DataType="Float32" BlockXSize="452" BlockYSize="21" />
      <SrcRect xOff="0" yOff="0" xSize="1354" ySize="2030" />
      <DstRect xOff="0" yOff="0" xSize="1354" ySize="2030" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

Same for the latitudes in lat.vrt:

<VRTDataset rasterXSize="1354" rasterYSize="2030">
  <SRS>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]</SRS>
  <VRTRasterBand dataType="Float32" band="1">
    <SimpleSource>
      <SourceFilename relativeToVRT="1">HDF5:A2015045060000.L2_LAC_OC.nc://navigation_data/latitude</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1354" RasterYSize="2030" DataType="Float32" BlockXSize="452" BlockYSize="21" />
      <SrcRect xOff="0" yOff="0" xSize="1354" ySize="2030" />
      <DstRect xOff="0" yOff="0" xSize="1354" ySize="2030" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

and for the data chlor.vrt:

<VRTDataset rasterXSize="1354" rasterYSize="2030">
       <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 band="1" datatype="Float32">  
    <SimpleSource>
      <SourceFilename relativeToVRT="1">HDF5:A2015045060000.L2_LAC_OC.nc://geophysical_data/chlor_a</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1354" RasterYSize="2030" DataType="Float32" BlockXSize="452" BlockYSize="21" />
      <SrcRect xOff="0" yOff="0" xSize="1354" ySize="2030" />
      <DstRect xOff="0" yOff="0" xSize="1354" ySize="2030" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

Do the warping with:

gdalwarp -geoloc -t_srs EPSG:4326 chlor.vrt chlor-out.tif

and the result fits to the shorelines around Borneo:

enter image description here

Alternatively to creating the vrts manually, you can create them with GDAL:

gdal_translate -of VRT HDF5:A2015045060000.L2_LAC_OC.nc://navigation_data/longitude lon.vrt
gdal_translate -of VRT HDF5:A2015045060000.L2_LAC_OC.nc://navigation_data/latitude lat.vrt
gdal_translate -of VRT HDF5:A2015045060000.L2_LAC_OC.nc://geophysical_data/chlor_a chlor.vrt

With a good text editor, remove the GCP lists from all of them, and insert only into the chlor.vrt this section instead:

   <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>

Then run

gdalwarp -geoloc -t_srs EPSG:4326 -overwrite chlor.vrt chlor-vrt.tif

to get the same picture as above.

Another solution, working with manually edited GCP points, can be found in my answer for Using GDALwarp for reprojecting netCDF file?

I want the output pixel size to be 2488 2141 because later on I'm gonna be using the gdal_calc which requires the image to have the same dimension. How will I do it? Where should I insert the -ts in gdalwarp? — lovelyvm, Jul 14 '15 at 23:03
You can automate the process by building the vrt files with python code, but I can't help you out on that. Since the swath images are rotated, the size will be different for every flight and day. — AndreJ, Jul 15 '15 at 04:17
Does that mean it isn't possible with batch process? Can I reduce the output pixel size by just adding -ts in gdalwarp? — lovelyvm, Jul 16 '15 at 00:49
I don't think you can write the GCP substitution into a batch file. For the -ts, just try it. You might need to add -srcnodata as well. — AndreJ, Jul 16 '15 at 03:54

Unable to warp HDF5 files

1 Answers1

Linked