2

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 9:24
  • Do you mean run gdalwarp with -geoloc first then run with -te again? Haven't tried it yet. – user Jul 14 '15 at 9:41
  • 1
    From my answer here gis.stackexchange.com/questions/128040/… 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. – user 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 5:19
3

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

  1. 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>
  1. 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>
  1. 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>
  1. 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? – user Jul 14 '15 at 23:03
  • Will this be possible in batch process? – user Jul 14 '15 at 23:21
  • 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 4:17
  • Does that mean it isn't possible with batch process? Can I reduce the output pixel size by just adding -ts in gdalwarp? – user Jul 16 '15 at 0: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 3:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.