1

I am trying to re-project a netCDF file using the gdalwarp command. However, the gdalwarp is giving me a black output.

Below is my code:

import gdal
import osr
import netCDF4
from osgeo.gdalconst import *
import os
import subprocess

gdal.AllRegister()
file_nc = 'input.nc'
ds = gdal.Open(file_nc, GA_ReadOnly)
Dataset = ds.GetSubDatasets()

substr = 'Rrs_655'
band = ''
names = [i[0] for i in Dataset]
for word in names[:]:
if word.endswith(substr):
band = word
b = 'gdal_translate -a_scale 0.000002 -a_offset 0.05 -a_nodata -32767 '+ band + ' output13.tif'
print(b)
subprocess.run(b, shell = True)
cmd = 'gdalwarp output13.tif outfile2.tif -t_srs "+proj=longlat +ellps=WGS84"'
subprocess.run(cmd, shell = True)

The gdalinfo of output13.tif band is:

Files: output13.tif
Size is 7821, 7931
Coordinate System is `'
GCP Projection =
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 7931.0)
Upper Right ( 7821.0, 0.0)
Lower Right ( 7821.0, 7931.0)
Center ( 3910.5, 3965.5)
Band 1 Block=7821x1 Type=Int16, ColorInterp=Gray
NoData Value=-32767
Offset: 0.05, Scale:2e-06
Metadata:
geophysical_data_Rrs_655_add_offset=0.050000001
geophysical_data_Rrs_655_long_name=Remote sensing reflectance at 655 nm
geophysical_data_Rrs_655_scale_factor=2e-06
geophysical_data_Rrs_655_solar_irradiance=1550.3799
geophysical_data_Rrs_655_standard_name=surface_ratio_of_upwelling_radiance_emerging_from_sea_water_to_downwelling_radiative_flux_in_air
geophysical_data_Rrs_655_units=sr^-1
geophysical_data_Rrs_655_valid_max=25000
geophysical_data_Rrs_655_valid_min=-30000
geophysical_data_Rrs_655__FillValue=-32767

The gdalinfo if the file outfile2.tif is:

Files: outfile2.tif
Size is 12259, 12307
Coordinate System is:
GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (-20696.575935146218399,3354.169516545325678)
Pixel Size = (1.950162663734099,-1.950162663734099)
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -20696.576, 3354.170) (Invalid angle,Invalid angle)
Lower Left ( -20696.576, -20646.482) (Invalid angle,Invalid angle)
Upper Right ( 3210.468, 3354.170) (Invalid angle,Invalid angle)
Lower Right ( 3210.468, -20646.482) (Invalid angle,Invalid angle)
Center ( -8743.054, -8646.156) (Invalid angle,Invalid angle)
Band 1 Block=12259x1 Type=Int16, ColorInterp=Gray
NoData Value=-32767
Metadata:
geophysical_data_Rrs_655_add_offset=0.050000001
geophysical_data_Rrs_655_long_name=Remote sensing reflectance at 655 nm
geophysical_data_Rrs_655_scale_factor=2e-06
geophysical_data_Rrs_655_solar_irradiance=1550.3799
geophysical_data_Rrs_655_standard_name=surface_ratio_of_upwelling_radiance_emerging_from_sea_water_to_downwelling_radiative_flux_in_air
geophysical_data_Rrs_655_units=sr^-1
geophysical_data_Rrs_655_valid_max=25000
geophysical_data_Rrs_655_valid_min=-30000
geophysical_data_Rrs_655__FillValue=-32767
Improve this question
8
  • By the corner coordinates ( Upper Left ( 0.0, 0.0) and so on) output13.tif it is not really georeferenced. That should be fixed first.
    – user30184
    Apr 17, 2019 at 18:45
  • Hello user30184, can you suggest how that can be done?
    – npal
    Apr 17, 2019 at 18:52
  • I am not so familiar with netCDF and without test data I can't do anything. Try with gdalinfo and see if even it can read the georeferencing from your subdatasets.
    – user30184
    Apr 17, 2019 at 19:13
  • Here is the link to the input.nc file: drive.google.com/file/d/1xw6OB0oZhOde0Lb4x3YJLe4hyyDfO2vu/… The commands that I am using are: 1. gdal_translate -a_scale 0.000002 -a_offset 0.05 -a_nodata -32767 HDF5:"input.nc"://geophysical_data/Rrs_655 output13.tif 2. gdalwarp output13.tif outfile2.tif -t_srs "+proj=longlat +ellps=WGS84" Thank you!
    – npal
    Apr 17, 2019 at 19:25
  • You seem to have HDF5 data and driver page gdal.org/frmt_hdf5.html says There is no universal way of storing georeferencing in HDF files. However, some product types have mechanisms for saving georeferencing, and some of these are supported by GDAL. Currently supported are (subdataset_type shown in parenthesis): HDF5 OMI/Aura Ozone (O3) Total Column 1-Orbit L2 Swath 13x24km (Level-2 OMTO3). I guess that you have other sort of data and GDAL can't georeference it and also gdalwarp will fail. The not-georeferenced tiff is probably the best that you can get but I may be wrong.
    – user30184
    Apr 17, 2019 at 19:59

1 Answer 1

Reset to default
1

Gdalwarp stumbles over the nodata values in the latitude and longitude bands of the netcdf file. The related bug issue is fixed in GDAL 2.4.2: https://github.com/OSGeo/gdal/issues/1451

As a workaround, I did this:

Extract the desired band to a vrt:

 gdal_translate -of VRT HDF5:"input.nc"://geophysical_data/Rrs_655 -a_nodata -32767 input.vrt

Open the file in a text editor and remove all GCP that have coordinates of

X="-3.276700000000E+04" Y="-3.276700000000E+04"

Use gdalwarp WITHOUT the -geoloc parameter:

 gdalwarp -t_srs EPSG:4326 input.vrt input.tif

to get this result: enter image description here

Improve this answer
13
  • I was wondering if the above method valid for batch re-projection of netCDF files?
    – npal
    Apr 21, 2019 at 18:26
  • Basically, yes. The vrt file is XML, and it should be possible to use an XML parser to eliminate the blocking GCP lines. But I can not help you on coding that in Python or C.
    – AndreJ
    Apr 21, 2019 at 18:45
  • Thank you for the solution. I wanted to ask how did you choose those blocking GCP lines? I am trying to automate the process and was wondering if you could tell me how to certain which GCP lines to exclude from the code.
    – npal
    Apr 22, 2019 at 13:25
  • You have to remove all GCP that have the nodata value -32767 for X and Y instead of useful degrees. For automation, you could read every GCP line, check for a valid range of +/-180, discard the GCPs outside of that, and write the others to a new vrt file.
    – AndreJ
    Apr 22, 2019 at 16:46
  • Thank you AndreJ it works like a charm. The GDAL library has a bug which does not read the GCP and nodata values properly for HDF5 files. I have reported the issue and it has been fixed in the version 2.4.2.
    – npal
    Apr 22, 2019 at 17:17

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.