Draw PNG with points from a shapefile using PIL

point_shp_2_png.py

import sys
from osgeo import gdal, ogr
from PIL import Image, ImageDraw

# Filename of input OGR file
vector_fn = 'RG2017_rgi_centroids.shp'

# Open the data source and read in the extent
datasource = ogr.Open(vector_fn)

if datasource is None:
  print ('Could not open file')
  sys.exit(1)

vector_layer = datasource.GetLayer()

mapXmin, mapXmax, mapYmin, mapYmax = vector_layer.GetExtent()
# http://gdal.org/python/osgeo.ogr.Layer-class.html#GetExtent

print ("mapXmin:", mapXmin)
print ("mapXmax:", mapXmax)
print ("mapYmin:", mapYmin)
print ("mapYmax:", mapYmax)

print(mapXmin, mapYmin, mapXmax, mapYmax)

print ("W:", mapXmax - mapXmin)
print ("H:", mapYmax - mapYmin)

# Define pixel_size
pixel_size = 0.01  # meters are one pixel
# Create the target data source
target_Width = int(abs(mapXmax - mapXmin) / pixel_size)
target_Height = int(abs(mapYmax - mapYmin) / pixel_size)


print ("target_Width:", target_Width)
print ("target_Height:", target_Height)


image = Image.new('RGBA', (target_Width, target_Height))
draw = ImageDraw.Draw(image)
r = 6  # radio
# Loop through the features in the layer
feature = vector_layer.GetNextFeature()
while feature:

    # get the x,y coordinates for the point
    geom = feature.GetGeometryRef()

    xx = geom.GetX()
    yy = geom.GetY()
    x = int((xx - mapXmin) / pixel_size)
    y = int((mapYmax - yy) / pixel_size)

    color = (18, 173, 42)  # dark green

    draw.ellipse((x - r, y - r, x + r, y + r), fill=color, outline=None)

    feature.Destroy()
    feature = vector_layer.GetNextFeature()

# close the data source and text file
datasource.Destroy()

image = image.resize((int(target_Width*0.5), int(target_Height*0.5)), Image.ANTIALIAS)
image.save('circles.png')