1

I was pulling some Himawari-8 images from RAMMB here and wanted to geo-reference the high res 5500 x 5500 images. I found a related question and answer using GDAL here which does work fine but I was hoping someone more familiar with GDAL and Rasterio could show me how to do the equivalent in Rasterio so I can write it fully in Python (I want to convert a lot of these images).

Here's an example of the working GDAL commands I wish to translate to Rasterio (slightly adapted from the linked question):

gdal_translate -a_srs "+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=140.7 +no_defs" -a_ullr -5500000 5500000 5500000 -5500000 full_disk_ahi_true_color_20210113014000.jpg temp.tif
gdalwarp -overwrite -t_srs "+proj=latlong +ellps=WGS84 +pm=140.7" -wo SOURCE_EXTRA=100 temp.tif Himawari8.tif

I'm hoping to get a Python snippet using Rasterio that generates the same geotif.

Owen Lamont
  • 13
  • 3
  • I added the GDAL commands to make this question more self-contained - I haven't attempted to implement this with Rasterio as I don't know where to begin with that. The GDAL commands work fine - I'd just prefer to be able to implement this in Python rather than calling out to a CLI application. – Owen Lamont Jan 16 '21 at 17:00
  • Thanks. Looks good, reopened. – Aaron Jan 16 '21 at 17:57

1 Answers1

1

Here's a gdal method (using python API, not CLI):

from osgeo import gdal

kwargs = {'format': 'VRT', 'outputSRS': '+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=140.7 +no_defs',
          'outputBounds': [-5500000, 5500000, 5500000, -5500000]}
src_ds = gdal.Translate('/vsimem/test.vrt', 'full_disk_ahi_true_color_20210116025000.jpg', **kwargs)


kwargs = {'format': 'GTiff', 'dstSRS': '+proj=latlong +ellps=WGS84 +pm=140.7', 'warpOptions': 'SOURCE_EXTRA=100'}
dst_ds = gdal.Warp('full_disk_ahi_true_color_20210116025000.tif', src_ds, **kwargs)
user2856
  • 65,736
  • 6
  • 115
  • 196
  • Thanks very much @user2856. I confirmed the GDAL API code produced identical results to the GDAL example code.

    I tested the Rasterio code too, and to my eye the transform it produced was somewhat similar but visually differently transformed tiff - albeit much smaller with a .aux.xml file. Looking at the result though in an image editor it was split the Earth image so there was a half Earth on each side and a lot of black space in the middle. When I dropped that geotiff as a layer into QGIS it didn't align with the original GDAL and GDAL API output.

     – Owen Lamont Jan 17 '21 at 04:47
  • I'll accept that answer since the GDAL API code sample solves my immediate problem.

    I had the misconception that being a higher level API Rasterio would have more succinct code but that doesn't appear to be the case. The Rasterio code look close but doesn't quite produce matching output to the raw GDAL calls.

     – Owen Lamont Jan 17 '21 at 04:50