0

I have been utilizing code that looks like the following to output GeoTIFFs for some time (rudimentary, but gets the job done):

def output_geotiff(data, outfile):

    # Get our dimensions
    cols = data.shape[1]
    rows = data.shape[0]

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(outfile, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((minlon, 0.01, 0, minlats, 0, 0.01))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(data)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

Where the function would output a GeoTIFF with EPSG:4326 at an equal lat/lon spacing of 0.01 degrees.

However, I now have a 2D grid that has dimensions of 30-arcseconds per point. I change the EPSG from 4326 to 1043, and then the values of 0.01 to 30. However, when I run this I get an error such that:

ERROR 1: PROJ: proj_create_from_database: crs not found

Is it not possible to output a GeoTIFF in arcseconds - does a conversion need to be done of sorts?

0

1 Answer 1

2

1043 is not a coordinate reference system (CRS) - EPSG.org (authoritative), EPSG.io.

Not all EPSG codes are CRSs, they can also identify units, datums, transformations and various other types.

There are 3600 arc seconds in a degree, so as your data is 30 arcsec resolution, just specify the resolution as 0.083333 (30/3600) and leave CRS as 4326.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.