Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active April 11, 2025 05:36
Show Gist options
  • Save mdsumner/12ee320594bbf613d4f939c523b61c43 to your computer and use it in GitHub Desktop.
Save mdsumner/12ee320594bbf613d4f939c523b61c43 to your computer and use it in GitHub Desktop.
path = 'https://eoimages.gsfc.nasa.gov/images/imagerecords/73000/73909/world.topo.bathy.200412.3x5400x2700.png'

import rioxarray
import rasterio
## assign by bbox (respected by open, when we open via rioxarray)
dsn_bbox = f'vrt:///vsicurl/{path}?a_ullr=-180,90,180,-90&a_srs=EPSG:4326'
ds = rasterio.open(dsn_bbox)
gt = ds.get_transform()
## ul corner x, x pixel width, shear, ul corner y, shear, y pixel height (from top)
#-180.00000000    0.06666667    0.00000000   90.00000000    0.00000000   -0.06666667

## assign by gt (respected by open, when we open via rioxarray)
a_gt = ",".join([str(g) for g in gt])
dsn_gt =   f'vrt:///vsicurl/{path}?a_gt={a_gt}'

## assign GCPs (not respected by open, we'll have to "warp")
## it does seem a bit funny to index 0,shape[0] 0,shape[1] but this seems to do the job and that is strictly where the outer margins are
dsn_gcp = f'vrt:///vsicurl/{path}?gcp=0,0,-180,90&gcp=5400,2700,180,-90'

## return the first and last xy values on xarray
def reportcoords(dsn):
  ds = rioxarray.open_rasterio(dsn)
  shape = ds.shape
  return([ds.x.values[[0, shape[2] - 1]], ds.y.values[[0, shape[1] - 1]]])


## these are the first and last centre-points as per xarray
reportcoords(dsn_bbox)
##[array([-179.96666667,  179.96666667]), array([ 89.96666667, -89.96666667])]

reportcoords(dsn_gt)
##[array([-179.96666667,  179.96666667]), array([ 89.96666667, -89.96666667])]

## nothing has happened, we have to "warp", it's really a no-op in terms of reprojection but resolves the georeferencing
reportcoords(dsn_gcp)
##[array([   0, 5399]), array([   0, 2699])]


from rasterio.vrt import WarpedVRT

vrt = WarpedVRT(rasterio.open(dsn_gcp), crs='EPSG:4326')
vrt.get_transform()
#[-180.0, 0.06666666666666667, 0.0, 90.0, 0.0, -0.06666666666666667]

## (not sure how to get that VRT text so we can open with rioxarray, but it's possible I usually use 
## the api more directly for that, I don't understand rasterio much)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment