0
  • I have a NetCDF with multiple climate data variables over an area for each day in a year, and I use it to generate a daily product for the area.

  • I have a set of GeoTIFFs covering a subset of that area, each GeoTIFF representing a day, containing different data.

  • I need to generate a new product that involves doing a calculation for each pixel in the GeoTIFF's area with the corresponding value from the product from the NetCDF.

  • Each pixel of the GeoTIFF has a spatial resolution of 3.125km x 3.125km and is projected in EPSG 3976 (NSIDC Sea Ice Polar Stereographic South).

  • The NetCDF is EPSG 4326 and comes with lat/lon grids within the structure. It has a native resolution of 0.28125 degrees (31km).

  • I have plotted them ontop of each other in QGIS (defining the NetCDF as EPSG 4326) and all seems in order.

I am struggling of where to start in developing a workflow so that I can make these GeoTIFFs and the NetCDF 'interact', so I can essentially multiply the output from my NetCDF product by the corresponding pixel in the GeoTIFF, at the GeoTIFF's resolution. I want to 'upsample' the NetCDF to match the GeoTIFF's resolution.

There are a lot of days of data so I want to script it, not do it manually using a tool. PyQGIS, Python, GDAL solutions would work.

Gdal info for the NetCDF:

Driver: netCDF/Network Common Data Format
Files: adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc
       adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc.aux.xml
Size is 512, 512
Metadata:
  NC_GLOBAL#Conventions=CF-1.6
  NC_GLOBAL#history=2021-04-15 19:39:59 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /cache/data8/adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc /cache/tmp/93fbbedb-7e77-4565-ba2b-05ef0e8b1c71-adaptor.mars.internal-1618515503.360069-11959-2-tmp.grib
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc":u10
  SUBDATASET_1_DESC=[720x13x85] u10 (16-bit integer)
  SUBDATASET_2_NAME=NETCDF:"adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc":v10
  SUBDATASET_2_DESC=[720x13x85] v10 (16-bit integer)
  SUBDATASET_3_NAME=NETCDF:"adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc":d2m
  SUBDATASET_3_DESC=[720x13x85] d2m (16-bit integer)
  SUBDATASET_4_NAME=NETCDF:"adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc":t2m
  SUBDATASET_4_DESC=[720x13x85] t2m (16-bit integer)
  SUBDATASET_5_NAME=NETCDF:"adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc":sp
  SUBDATASET_5_DESC=[720x13x85] surface_air_pressure (16-bit integer)
  SUBDATASET_6_NAME=NETCDF:"adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc":ssrd
  SUBDATASET_6_DESC=[720x13x85] surface_downwelling_shortwave_flux_in_air (16-bit integer)
  SUBDATASET_7_NAME=NETCDF:"adaptor.mars.internal-1618515503.3594973-11959-19-93fbbedb-7e77-4565-ba2b-05ef0e8b1c71.nc":strd
  SUBDATASET_7_DESC=[720x13x85] strd (16-bit integer)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)  

And for the geotiff:

Driver: GTiff/GeoTIFF
Files: clip_asi-AMSR2-s3125-20200301-v5.4.tif
Size is 117, 163
Coordinate System is:
PROJCRS["IDL GeoTIFF Suport
Projection = Polar Stereographic
True scale:      -70.0000deg
Gunnar Spreen, Apr 2004",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Polar Stereographic (variant B)",
        METHOD["Polar Stereographic (variant B)",
            ID["EPSG",9829]],
        PARAMETER["Latitude of standard parallel",-70,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8832]],
        PARAMETER["Longitude of origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8833]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",north,
            MERIDIAN[90,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            MERIDIAN[0,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            ORDER[2],
            LENGTHUNIT["metre",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (-1810937.500000000000000,-554687.500000000000000)
Pixel Size = (3125.000000000000000,-3125.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2020:03:02 05:15:07
  TIFFTAG_DOCUMENTNAME=/ssmi/www/htdocs/data/amsr2/asi_daygrid_swath/s3125/2020/mar/Amundsen/asi-AMSR2-s3125-20200301-v5.4.tif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=IDL 8.1, ITT Visual Information Solutions
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-1810937.500, -554687.500) (107d 1'47.10"W, 72d38'40.82"S)
Lower Left  (-1810937.500,-1064062.500) (120d26'14.82"W, 70d47' 7.34"S)
Upper Right (-1445312.500, -554687.500) (110d59'45.57"W, 75d46'46.47"S)
Lower Right (-1445312.500,-1064062.500) (126d21'39.97"W, 73d32'29.83"S)
Center      (-1628125.000, -809375.000) (116d25'58.56"W, 73d19'46.09"S)
4
  • gdalwarp all datasets to the same origin and spacing and then use raster math to multiply?
    – wingnut
    Apr 17 at 3:49
  • 1
    Does this answer your question? Raster Subtraction using rasterio or gdal
    – snowman2
    Apr 22 at 0:29
  • It keeps getting marked as a duplicate, but the question it supposedly duplicates 1) was posted after this one 2) just asks about raster subtraction, when the key issue here is them being different resolutions, projections and file formats.
    – Beardsley
    Apr 22 at 1:12
  • I updated the answer in this question to address the differences. Hopefully it is helpful. I will remove my close vote due to the differences you mentioned.
    – snowman2
    Apr 22 at 2:21
0

Here is my workflow of how I did this in case it is useful to anyone else:

  • Open relevant NETCDF layers and do calculations to give me the product I need to use with the geotiff. [In python]
  • Save this product as a new netcdf with the lat/lon data from the original netcdf.[In python]
  • Use gdal_translate to convert the new netcdf to multiband geotiff, and add EPSG 4326 projection info to it.
  • Use gdalwarp to convert the new geotiff's projection to EPSG 3976.
  • Use gdalwarp to convert the new geotiff's spatial resolution to the same as the other geotiff.
  • Use gdalwarp to set the corner coordinates of the new geotiff to the same as those of the other geotiff, so that the pixels are aligned.
  • Use QGIS's 'Clip raster by extent' to clip the new geotiff to the same dimensions as the other one, so that they have the same dimensions/number of pixels. (This could be done with gdal too but there was a problem for some reason.
  • Use python/numpy to do my calculations I needed - between the netcdf output (i.e. the new geotiff) and the other geotiff.
0

Step 1: Use rio.reproject_match to ensure they have the same resolution/projection/boundary.

import rioxarray
raster = rioxarray.open_rasterio("path_to_raster.tif", masked=True)
netcdf = rioxarray.open_rasterio("path_to_netcdf.nc", masked=True)
netcdf.rio.write_crs("EPSG:4326", inplace=True)
netcdf_match = netcdf.rio.reproject_match(raster)

Note: It is recommended to use assign_coords after rio.reproject_match to make the coordinates the exact same due to tiny differences in the coordinate values due to floating precision (issue 298).

netcdf_match = netcdf_match.assign_coords({
    "x": raster.x,
    "y": raster.y,
})

Step 2: Do calculations

calculated = netcdf_match["u10"] - raster
calculated.rio.to_raster(output_path)

Here is another example of raster calculations using reproject_match to do raster calculations: https://carpentries-incubator.github.io/geospatial-python/07-raster-calculations/index.html

10
  • You have applied this same answer to two questions. Have you considered whether the questions might be duplicates?
    – PolyGeo
    Apr 22 at 0:06
  • They are definitely very similar. Could be considered a duplicate.
    – snowman2
    Apr 22 at 0:18
  • I don't know the subject matter well enough to propose that they are duplicates, and if I vote that way my vote is binding. You have the vote to close and re-open privilege so if you voted for a duplicate then you would be able to get up to four second opinions before it was closed as one.
    – PolyGeo
    Apr 22 at 0:23
  • Since neither have been selected as answered, are there any issues marking as duplicate?
    – snowman2
    Apr 22 at 0:27
  • What do you mean by "selected as answered"? This particular question has been both answered and that answer has been accepted as answered. There's an FAQ on duplicates at meta.stackexchange.com/questions/10841/…
    – PolyGeo
    Apr 22 at 0:34

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.