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>
Improve this question
5
  • Does it work if you separate -geoloc and -te in two commands?
    – AndreJ
    Jul 14, 2015 at 9:24
  • Do you mean run gdalwarp with -geoloc first then run with -te again? Haven't tried it yet.
    – lovelyvm
    Jul 14, 2015 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, 2015 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, 2015 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, 2015 at 5:19

1 Answer 1

Reset to default
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?

Improve this answer
5
  • 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, 2015 at 23:03
  • Will this be possible in batch process?
    – lovelyvm
    Jul 14, 2015 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, 2015 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?
    – lovelyvm
    Jul 16, 2015 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, 2015 at 3:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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