pro project_lambert, lon=lon_, lat=lat_, x=x,y=y,kmperpix=kmperpix, minmaplat=minmaplat,center_lon=center_lon_,inverse=inverse, south=south

    if n_elements(center_lon_) eq 0 then center_lon_=0
    screen_x=map_size_lambert(kmperpix=kmperpix,minmaplat=minmaplat,maxrho=maxrho)
    screen_y=screen_x
    
    ;  print,'screen_x:',screen_x,' kmperpix:',kmperpix
    ;  print,' minmaplat:',minmaplat, ' maxrho:',maxrho
    
    ;If the south flag is undefined then set it automatically if possible
    ; TODO Remove the following line and make calling procedure set south if needed
    if n_elements(south) eq 0 && n_elements(lat_) gt 0 && max(lat_,/nan) le 0 then south = 1
    
    if keyword_set(south) then BEGIN
;       center_lon=-center_lon_ -180 
       center_lon = -center_lon_
    ENDIF ELSE BEGIN
       center_lon = center_lon_
    ENDELSE
    
    if ~keyword_set(inverse) then begin
        if keyword_set(south) then begin
            lat=-lat_
            lon=-lon_
        end else begin
            lat=lat_
            lon=lon_
        end
        rho=2d*sin((90d -lat)/2d*!const.dtor)/maxrho*screen_x/2d
        theta=(lon-center_lon-90d)*!const.dtor
        x=rho*cos(theta)+screen_x/2d
        y=rho*sin(theta)+screen_y/2d
    end else begin
        x0=x-screen_x/2d
        y0=y-screen_y/2d
        theta=atan(y0,x0)
        rho=sqrt(x0^2+y0^2)
        lon=theta/!const.dtor+center_lon+90d
        ;Constrain longitude to +/-180deg
        w=where(lon gt 180d,nw)
        if nw gt 0 then lon[w]-=360
        w=where(lon lt -180d,nw)
        if nw gt 0 then lon[w]+=360
        lat=90d -2d*asin(rho*maxrho/screen_x)/!const.dtor
        if keyword_set(south) then begin
            lat_=-lat
            lon_=-lon
        end else begin
            lat_=lat
            lon_=lon
        end
    end
    
end