I'm trying to work out the details of using pyproj to translate the lat lon coordinates from the USGS Earthquake feed to xy coordinates corresponding to those positions on a basemap of a determined width and height and origin at the upper left. My code so far naively uses PIL/pillow to translate lon lat to xy positions and draw circles without taking into account the EPSG:4326 projection of the basemap:
from PIL import Image from PIL import ImageDraw # open basemap image file basemap = Image.open(basemap_path).convert('RGBA) # resize to desired map size basemap.thumbnail(width, height, Image.LANCZOS) # get proportional height width_bmp, height_bmp = basemap.size # create background frame and paste basemap on it img=Image.new('RGB',(width, height), color = '#000000') img.paste(basemap, (0,0), basemap) draw = ImageDraw.Draw(img, 'RGBA') width_scale = width/360 height_scale = height_bmp/180 # usgs data has been parsed into a list for quake in earthquake_list: lon = float(quake["longitude"]) lat = float(quake["latitude"]) mag = float(quake["mag"]) # want to use pyproj to translate coordinates here instead of the following cx = (lon + 180) * width_scale cy = (90-lat) * height_scale r = scaleRadius(mag) # draw earthquake circles draw.ellipse((cx-r, cy-r, cx+r, cy+r), fill = colormap(mag)) draw = ImageDraw.Draw(img) img.save(filepath, quality=100)
I've gone through the documentation and have a general sense that I would use something like:
cx, cy = pyproj.transform("EPSG:4326", "xy", lon, lat)
with some width and height dimensions. Running the line above gives the error:
TypeError: p1 must be a Proj class
Update I do get values for:
p = Proj(proj = 'longlat', ellps='WGS84') cx,cy = p(lon, lat)
In this case a (lon, lat) input of
(-110, 39) gives the output
1.923, 0.688 which isn't right. What I need is a projection function that takes the lon, lat and the width and height of the image and outputs the x, y coordinates. I would write this on my own, but I can't find a suitable equation for that conversion and was hoping to find this capability in proj.
UPDATE Here is the basemap I am using.
And here is the resulting drawing. I'm assuming the centerpoints for the circles marking earthquakes are off, because I am using Cartesian coordinates for lon, lat, but the basemap was created with a projection and thus must have a distortion from x, y Cartesian coordinates.
Also, I realize I need to put my legend on a logarithmic scale. Haven't gotten to that yet.