GDAL – How to Convert NetCDF to GeoTIFF Upside Down

climate-data-operatorsgdalgdal-translatenetcdf

I want to convert monthly precipitation netCDF file, downloaded from TerraClimate data to GeoTIFF and available as 1 file per year. Example: http://thredds.northwestknowledge.net:8080/thredds/fileServer/TERRACLIMATE_ALL/data/TerraClimate_ppt_2020.nc (+-103MB).

To do the conversion for each month, I use CDO to select the date. Below is the complete script:

for t in `cdo showdate TerraClimate_ppt_2020.nc`; do
   cdo seldate,$t TerraClimate_ppt_2020.nc dummy.nc
   gdal_translate -of GTiff -a_ullr -180 90 180 -90 -a_srs EPSG:4326 dummy.nc ../tif/TerraClimate_ppt_$t.tif
done

Unfortunately the result when I open to QGIS is upside down. Is there any configuration options from CDO or GDAL that I should use in above script to avoid the upside down result?

Below is the metadata from ncdump and gdalinfo, also result in qgis.

Best Answer

I think gdalmdimtranslate is a perfect use case for your problem (New in version 3.1).

Reposting what worked for you

gdalmdimtranslate -of GTiff -co COMPRESS=LZW -co PREDICTOR=1 dummy.nc ../TerraClimate_ppt_$t.tif

Related Solutions

[GIS] Create CF-compliant NetCDF from ESRI Grid

In ArcGIS 10.1, I just tried converting two rasters (geographic and albers) using ArcToolbox=>Multidimensional Tools=>Raster to NetCDF and the resulting files seem to be CF-Compliant -- they work just fine with the Unidata THREDDS Data Server.

Try drilling down to the OPeNDAP data links or Godiva2 viewer links here: http://geoport.whoi.edu/thredds/catalog/usgs/data2/rsignell/data/bathy/arc_netcdf/catalog.html to see these two netCDF files.

And of course you can just download them using the HTTPServer links.

Here's what the geographic file looks like in Unidata's Tools-UI: enter image description here

But oh no, here's what the Albers file looks like: enter image description here

Things don't look good. The topography is not plotting in the right location -- it's shifted by quite a bit.

Taking a look at the ArcGIS 10.1-produced NetCDF file, we see that there is no ellipsoid information in the grid_mapping variable, in this case, the variable albers_conical_equal_area:

netcdf dods://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/data/bathy/arc_netcdf/traster_albers.nc {
  dimensions:
    x = 1272;
    y = 1361;
  variables:
    double x(x=1272);
      :long_name = "x coordinate of projection";
      :standard_name = "projection_x_coordinate";
      :units = "Meter";

    double y(y=1361);
      :long_name = "y coordinate of projection";
      :standard_name = "projection_y_coordinate";
      :units = "Meter";

    int albers_conical_equal_area;
      :grid_mapping_name = "albers_conical_equal_area";
      :longitude_of_central_meridian = -96.0; // double
      :latitude_of_projection_origin = 23.0; // double
      :false_easting = 0.0; // double
      :false_northing = 0.0; // double
      :standard_parallel = 29.5, 45.5; // double

    float traster_View_ProjectRaster(y=1361, x=1272);
      :_CoordinateAxes = "x y y x ";
      :long_name = "traster_View_ProjectRaster";
      :esri_pe_string = "PROJCS[\"USA_Contiguous_Albers_Equal_Area_Conic_USGS_version\",GEOGCS[\"GCS_North_American_1983\",DATUM[\"D_North_American_1983\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Albers\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-96.0],PARAMETER[\"Standard_Parallel_1\",29.5],PARAMETER[\"Standard_Parallel_2\",45.5],PARAMETER[\"Latitude_Of_Origin\",23.0],UNIT[\"Meter\",1.0]]";
      :coordinates = "x y";
      :grid_mapping = "albers_conical_equal_area";
      :units = "Meter";
      :missing_value = -3.4028235E38f; // float

  // global attributes:
  :Conventions = "CF-1.0";
  :Source_Software = "Esri ArcGIS";
}

The CF Conventions specify that the ellipsoid representation should be encoded using semi_major_axis and (semi_minor_axis or inverse_flattening). Since semi_major_axis and inverse_flattening values are right there in the esri_pe_string attribute, we can fix this file up using NcML thusly:

<netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2"
  location="traster_albers.nc">
<variable name="albers_conical_equal_area">
    <attribute name="semi_major_axis" type="double" value="6378137.0"/>
    <attribute name="inverse_flattening" type="double" value="298.257222101"/>
</variable>
</netcdf>

If we try plotting again, it now looks correct: enter image description here

This should be easy for ESRI to fix because all they have to do is pull the ellipsoid values out of the esri_pe_string. The resulting CF-Compliant file would then look like this:

netcdf dods://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/data/bathy/arc_netcdf/traster_albers_fixed2.ncml {
  dimensions:
    x = 1272;
    y = 1361;
  variables:
    double x(x=1272);
      :long_name = "x coordinate of projection";
      :standard_name = "projection_x_coordinate";
      :units = "Meter";

    double y(y=1361);
      :long_name = "y coordinate of projection";
      :standard_name = "projection_y_coordinate";
      :units = "Meter";

    int albers_conical_equal_area;
      :grid_mapping_name = "albers_conical_equal_area";
      :longitude_of_central_meridian = -96.0; // double
      :latitude_of_projection_origin = 23.0; // double
      :false_easting = 0.0; // double
      :false_northing = 0.0; // double
      :standard_parallel = 29.5, 45.5; // double
      :semi_major_axis = 6378137.0; // double
      :inverse_flattening = 298.257222101; // double

    float traster_View_ProjectRaster(y=1361, x=1272);
      :_CoordinateAxes = "x y y x ";
      :long_name = "traster_View_ProjectRaster";
      :esri_pe_string = "PROJCS[\"USA_Contiguous_Albers_Equal_Area_Conic_USGS_version\",GEOGCS[\"GCS_North_American_1983\",DATUM[\"D_North_American_1983\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Albers\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-96.0],PARAMETER[\"Standard_Parallel_1\",29.5],PARAMETER[\"Standard_Parallel_2\",45.5],PARAMETER[\"Latitude_Of_Origin\",23.0],UNIT[\"Meter\",1.0]]";
      :coordinates = "x y";
      :grid_mapping = "albers_conical_equal_area";
      :units = "Meter";
      :missing_value = -3.4028235E38f; // float

  // global attributes:
  :Conventions = "CF-1.0";
  :Source_Software = "Esri ArcGIS";
}

gdal – How to Add Alpha Channel to MBTiles Using gdal_translate

For the issue "turned pngs", you may be able to use https://github.com/roblabs/gdal2tilesp

which has support for xyz, from the notes: Generates Leaflet layer for use with the XYZ default format. Where the tms parameter is set to false

I've also tried gdal2tiles.py and mb-util utilities. they made good mbtile from point of transparency but all the pngs are turned upside down see figure2 (gdalinfo shows that origin has wrong Y coordinate with "-")

python gdal2tilesp.py -h

-o OUTPUT_CACHE, --output=OUTPUT_CACHE
                    Format for output cache. Values allowed are tms and
                    xyz, being xyz the default value

Best Answer

Related Solutions

[GIS] Create CF-compliant NetCDF from ESRI Grid

gdal – How to Add Alpha Channel to MBTiles Using gdal_translate

Related Question