LL2LAEA Subroutine

public subroutine LL2LAEA(x, y, lon_in, lat_in, projection_attr)

Arguments

Type IntentOptional Attributes Name
real, intent(out) :: x
real, intent(out) :: y
real, intent(in) :: lon_in
real, intent(in) :: lat_in
double precision, intent(in) :: projection_attr(10)

Called by

proc~~ll2laea~~CalledByGraph proc~ll2laea LL2LAEA proc~ll2proj LL2PROJ proc~ll2proj->proc~ll2laea proc~read_country_bounding_box_data read_country_bounding_box_data proc~read_country_bounding_box_data->proc~ll2laea proc~uemep_preaggregate_shipping_asi_data uEMEP_preaggregate_shipping_asi_data proc~uemep_preaggregate_shipping_asi_data->proc~ll2laea proc~uemep_read_industry_data uEMEP_read_industry_data proc~uemep_read_industry_data->proc~ll2laea proc~uemep_read_netcdf_population uEMEP_read_netcdf_population proc~uemep_read_netcdf_population->proc~ll2laea proc~uemep_read_receptor_data uEMEP_read_receptor_data proc~uemep_read_receptor_data->proc~ll2laea proc~uemep_read_roadlink_data_ascii uEMEP_read_roadlink_data_ascii proc~uemep_read_roadlink_data_ascii->proc~ll2laea proc~uemep_read_shipping_asi_data uEMEP_read_shipping_asi_data proc~uemep_read_shipping_asi_data->proc~ll2laea proc~uemep_set_subgrid_select_latlon_centre uEMEP_set_subgrid_select_latlon_centre proc~uemep_set_subgrid_select_latlon_centre->proc~ll2laea proc~uemep_calculate_emissions_for_emep uEMEP_calculate_emissions_for_EMEP proc~uemep_calculate_emissions_for_emep->proc~uemep_read_industry_data proc~uemep_calculate_emissions_for_emep->proc~uemep_read_roadlink_data_ascii proc~uemep_calculate_emissions_for_emep->proc~uemep_read_shipping_asi_data proc~uemep_region_mask_new uEMEP_region_mask_new proc~uemep_region_mask_new->proc~ll2proj proc~uemep_subgrid_emep_from_in_region uEMEP_subgrid_EMEP_from_in_region proc~uemep_subgrid_emep_from_in_region->proc~ll2proj program~uemep uEMEP program~uemep->proc~read_country_bounding_box_data program~uemep->proc~uemep_preaggregate_shipping_asi_data program~uemep->proc~uemep_read_industry_data program~uemep->proc~uemep_read_netcdf_population program~uemep->proc~uemep_read_receptor_data program~uemep->proc~uemep_read_roadlink_data_ascii program~uemep->proc~uemep_read_shipping_asi_data program~uemep->proc~uemep_set_subgrid_select_latlon_centre program~uemep->proc~uemep_calculate_emissions_for_emep program~uemep->proc~uemep_region_mask_new program~uemep->proc~uemep_subgrid_emep_from_in_region

Source Code

    subroutine LL2LAEA(x, y, lon_in, lat_in, projection_attr)
        ! https://epsg.io/3035
        ! grid_mapping_name = lambert_azimuthal_equal_area
        ! Map parameters:
        ! longitude_of_projection_origin
        ! latitude_of_projection_origin
        ! false_easting - This parameter is optional (default is 0)
        ! false_northing - This parameter is optional (default is 0)lat_stand1_lambert=projection_attributes(1)
        double precision, intent(in) :: projection_attr(10)
        real, intent(in) :: lon_in, lat_in
        real, intent(out) :: x, y

        ! Local variables
        real :: semi_major_axis
        real :: deg2rad,rad2deg
        real :: lat0, lat0_in, lon0, lon0_in
        real :: false_easting, false_northing
        real :: lon, lat
        real :: inv_f, f, a, e, q_p, q_0, q, beta, beta_0, R_q, D, B
        
        lon0_in = projection_attr(1)
        lat0_in = projection_attr(2)
        false_easting = projection_attr(3)
        false_northing = projection_attr(4)
        semi_major_axis = projection_attr(5)
        inv_f = projection_attr(6) !flattening

        deg2rad = PI/180.0
        rad2deg = 180.0/PI
        a = semi_major_axis
        f = 1.0/inv_f
        e = sqrt(2.0*f - f*f)

        lat0 = lat0_in*deg2rad
        lon0 = lon0_in*deg2rad
        lon = lon_in*deg2rad
        lat = lat_in*deg2rad

        q_p = (1.0 - e*e)*(1.0/(1.0-e*e) - 1.0/(2.0*e)*log((1.0 - e)/(1.0 + e)))
        q_0 = (1.0 - e*e)*(sin(lat0)/(1.0 - e*e*sin(lat0)**2) - 1.0/(2.0*e)*log((1.0 - e*sin(lat0))/(1.0 + e*sin(lat0))))
        q = (1.0 - e*e)*(sin(lat)/(1.0 - e*e*sin(lat)**2) - 1.0/(2.0*e)*log((1.0 - e*sin(lat))/(1.0 + e*sin(lat))))
        beta_0 = asin(q_0/q_p)
        beta = asin(q/q_p)
        R_q = a*sqrt(q_p/2.0)
        D = a*(cos(lat0)/sqrt(1.0 - e*e*sin(lat0)**2)/(R_q*cos(beta_0)))
        B = R_q*sqrt(2.0/(1.0 + sin(beta_0)*sin(beta) + cos(beta_0)*cos(beta)*cos(lon - lon0)))
        x = false_easting + B*D*cos(beta)*sin(lon - lon0)
        y = false_northing + B/D*(cos(beta_0)*sin(beta) - sin(beta_0)*cos(beta)*cos(lon - lon0))
    end subroutine LL2LAEA