3

I am trying to save a 2D numpy array as a raster, and then split the raster using shapefiles. When I create the raster, I either have the problem of the latitude being flipped upside down, or that the shapefiles won't plot on the basemap.

Code - Shapefiles plot correctly on basemap, though latitude is flipped upside down.

# Create raster
xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]
nrows, ncols = np.shape(array)
geotransform = ([xmin, 0.1, 0, ymax, 0, -0.1])
output_raster = gdal.GetDriverByName('GTiff').Create('raster.tif', ncols, nrows, 1, gdal.GDT_Float32)  
output_raster.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
output_raster.SetProjection( srs.ExportToWkt() )
output_raster.GetRasterBand(1).WriteArray(array)
output_raster.FlushCache()
output_raster = None

# Plot
ds = gdal.Open('raster.tif')
zz = ds.ReadAsArray()
gt = ds.GetGeoTransform()
x = np.linspace(gt[0], gt[0] + gt[1] * zz.shape[1], zz.shape[1])
y = np.linspace(gt[3], gt[3] + gt[5] * zz.shape[0], zz.shape[0])
xx, yy = np.meshgrid(x, y)
plt.figure(1)
plt.subplot(111)
plt.tight_layout()
ax1 = plt.gca()
m = Basemap(llcrnrlon=xx.min(), llcrnrlat=yy.min(), urcrnrlon=xx.max(), urcrnrlat=yy.max(), projection = 'cyl')
m.readshapefile('shapefile', 'states')
cs = m.contourf(xx, yy, zz)

Code change to make latitude correct, though shapefiles do not plot at all on basemap.

geotransform = ([xmin, 0.1, 0, ymax, 0, 0.1])

I am after both the correct latitude and the shapefiles to plot.

Can anyone help?

EDIT: The next step in the script is clip the raster data with a shapefile, creating new rasters for each polygon (Code below). When I use the first method above these new rasters contain the data (though upside down), and when using the second method the data in them is all 0.

# Creates tif files for each state polygon
inraster = 'raster.tif'
inshape = 'shapefile'

ds = ogr.Open(inshape)
lyr = ds.GetLayer(0)

lyr.ResetReading()
ft = lyr.GetNextFeature()

while ft:

    country_name = ft.GetFieldAsString('NAME_1')

    outraster = inraster.replace('.tif', '_%s.tif' % country_name.replace(' ', '_'))    
    subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, 
                 '-crop_to_cutline', '-cwhere', "'NAME_1'='%s'" % country_name])

    ft = lyr.GetNextFeature()

ds = None
  • Wondering if reversing the sign of gt[5] in y = np.linspace(gt[3], gt[3] + gt[5] * zz.shape[0], zz.shape[0]) could help. – user30184 Aug 7 '16 at 12:56
  • Still the same issue. Also, this change wouldn't affect the raster data itself, only the basemap. My next step in the script is to clip the raster data using a shapefile, so the solution I am after is within the raster creation. – luke.conibear Aug 8 '16 at 4:38
  • Is the raster.tif file flipped also if you open it with QGIS? – user30184 Aug 8 '16 at 4:49
  • When I open both rasters in QGIS, they are both the same way up as when displayed in basemap. – luke.conibear Aug 8 '16 at 5:29
  • 1
    I am not sure but I think you have two alternatives: a) Use your array as it is, set the origin in the GeoTransform to lower left corner, and use positive 0.1 as the north step size or b) Flip the array so that the first row is the north-most one, keep origin at top-left and use negative step size -0.1. I am remembering that is some parts of the code GDAL may automatically change the positive y-increment to negative (does not believe that user really can have a south up image which may make alternative a) to fail. – user30184 Aug 9 '16 at 8:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.