Traducir EPSG: 4326 coordenadas lon lat a xy usando pyproj
Estoy tratando de resolver los detalles del uso de pyproj para traducir las coordenadas lat lon de la alimentación del terremoto de USGS a las coordenadas xy correspondientes a esas posiciones en un mapa base de un ancho, alto y origen determinados en la parte superior izquierda. Mi código hasta ahora usa ingenuamente PIL / pillow para traducir long lat a posiciones xy y dibujar círculos sin tener en cuenta la proyección EPSG: 4326 del mapa base:
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)
Revisé la documentación y tengo la sensación general de que usaría algo como:
cx, cy = pyproj.transform("EPSG:4326", "xy", lon, lat)
con algunas dimensiones de ancho y alto. Ejecutar la línea de arriba da el error:
TypeError: p1 must be a Proj class
Actualización obtengo valores para:
p = Proj(proj = 'longlat', ellps='WGS84')
cx,cy = p(lon, lat)
En este caso, una entrada (lon, lat) de (-110, 39)da la salida 1.923, 0.688que no es correcta. Lo que necesito es una función de proyección que tome la longitud, la latitud y la anchura y la altura de la imagen y genere las coordenadas x, y. Escribiría esto por mi cuenta, pero no puedo encontrar una ecuación adecuada para esa conversión y esperaba encontrar esta capacidad en proj.
ACTUALIZAR Aquí está el mapa base que estoy usando.
Y aquí está el dibujo resultante. Supongo que los puntos centrales de los círculos que marcan terremotos están apagados, porque estoy usando coordenadas cartesianas para lon, lat, pero el mapa base se creó con una proyección y, por lo tanto, debe tener una distorsión de las coordenadas cartesianas x, y.
Además, me doy cuenta de que necesito poner mi leyenda en una escala logarítmica. Aún no he llegado a eso.
Respuestas
Cuando superpongo su mapa base en una proyección equirectangular, encaja perfectamente. Entonces, esto significa que su mapa base realmente solo usa la simple "proyección" de trazar longitudes como xy latitudes como y. Lo que significa que lo único de lo que debe preocuparse es verificar el tamaño de la imagen, lo que ya hace en sus fórmulas.
cx = (lon + 180) * width_scale
cy = (90-lat) * height_scale
Estos funcionarán correctamente para sus necesidades, al píxel más cercano. De hecho, verifiqué los terremotos en USGS en esa fecha, ¡y sus puntos están todos en la ubicación correcta!