0

I am trying to map global greenhouse gas emissions data from the European Commission's Emissions Database for Global Atmospheric Research (EDGAR) and Northern Arizona University's Fossil Fuel Data Assimilation System (FFDAS) in QGIS. I am also trying to perform a spatial intersection between a rasterized/polygonized version of these layers (individually) with census designated boundary (CDB) shapefiles, which use the Lambert Conformal Conic projection.

However, the netCDF (.nc) data files from EDGAR and FFDAS data repositories linked are not georeferenced, and the FFDAS emissions map is flipped upside down. To address this, I've manually set the layer CRS of both the netCDF layers and the shape shapefile (.shp) layer to EPSG:4326 - WGS 84 (default setting). But, this has the unintended effect of making the global emissions maps and the county census shapefile map appear in wildly different locations, with no overlap, and also appear in unrealistic sizes. Not only is a single county many times larger than the global map, the county appears many "worlds" away. The problem is visualized below.

The discrepancy with the EDGAR file: EDGAR

The discrepancy with the FFDAS file: FFDAS


Reprojection

When I reproject the shapefile layer (Vector > Data Management Tools > Reproject Layer...), from Lambert Conformal Conic to WGS 84 (the target CRS), the county becomes a more reasonable size but is still in a different plane.

FFDAS2

When I attempt to assign the EDGAR/FFDAS layer as input to project to WGS 84 (by Raster > Projections > Assign Projections...), I get the following error:

ERROR 6: The NETCDF driver does not support update access to existing datasets.
Process returned error code 255

When I try to do the same thing with the Warp method, I get the following error:

ERROR 1: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for /Users/toshiro/Desktop/Projects/CCAR/Data/FFDAS/ffdas_flux_2010.nc. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.
Process returned error code 1

I also tried to reproject the EDGAR emissions using the below Python code (as provided by @Keagan Allan). I used a variation that replaced 'lon' with 'longitude' to reproject the FFDAS emissions file.

import xarray as xr
import rioxarray
    
with xr.open_dataset(r'Filepath') as src:
    ds = src
    ds = ds.rio.write_crs("epsg:4326", inplace=True)
    print (ds.coords['lon'])
    ds.coords['lon'] = (ds.coords['lon'] + 180) % 360 - 180
    ds = ds.sortby(ds.lon)
    print (ds.coords['lon'])

    ds.to_netcdf(r'Filepath')

While this makes the EDGAR file larger, and moves it in a slightly different location, it still does not overlap with the reprojected county shapefile (shown below).

EDGAR reprojected

For the FFDAS file, no visible change is apparent: the map is still upside down and there is no overlap with the county shapefile (shown below).

FFDAS reprojected


Metadata

I found in Panoply the following relevant metadata for FFDAS (non-reprojected, as the original file seems to be of the correct projection, per @Keagan):

dimensions:
    latitude = 1800;
    latitude_edge = 1801;
    longitude = 3600;
    longitude_edge = 3601;
  variables:
    float latitude(latitude=1800);
      :units = "degrees_n";

    float latitude_edge(latitude_edge=1801);
      :units = "degrees_n";

    float longitude(longitude=3600);
      :units = "degrees_e";

    float longitude_edge(longitude_edge=3601);
      :units = "degrees_e";

Likewise, here is the relevant metadata for the (reprojected) EDGAR file:

dimensions:
    lat = 1800;
    lon = 3600;
  variables:
    float lon(lon=3600);
      :_FillValue = NaNf; // float

    float emi_co2(lat=1800, lon=3600);
      :standard_name = "tendency_of_atmosphere_mass_content_of_carbon_dioxide_due_to_emission";
      :long_name = "Emissions of CO2 - ";
      :units = "kg m-2 s-1";
      :cell_method = "time: mean (interval: 1 year, 365 days)";
      :total_emi_co2 = "   3.41281e+013 kg/year";
      :comment = " (see http://edgar.jrc.ec.europa.eu/methodology.php#12sou for the definitions of the single sources)";
      :_FillValue = NaNf; // float
      :grid_mapping = "spatial_ref";

    float lat(lat=1800);
      :standard_name = "latitude";
      :long_name = "latitude";
      :units = "degrees_north";
      :comment = "center_of_cell";
      :_FillValue = NaNf; // float

    int spatial_ref;
      :crs_wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]";
      :semi_major_axis = 6378137.0; // double
      :semi_minor_axis = 6356752.314245179; // double
      :inverse_flattening = 298.257223563; // double
      :reference_ellipsoid_name = "WGS 84";
      :longitude_of_prime_meridian = 0.0; // double
      :prime_meridian_name = "Greenwich";
      :geographic_crs_name = "WGS 84";
      :grid_mapping_name = "latitude_longitude";
      :spatial_ref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]";

Are there any options to solve this issue? How can I properly georeference these .nc files such that the county shapefile overlaps with its correct location on the global map, and with accurate proportions?

6
  • 1
    This is a variant of the most common question here in GIS SE. If two datasets which should line up don't, it's because one or both have had their coordinate reference system metadata mangled. The cause is nearly always someone clobbering the correct metadata, though it's not always the fault of the OP. Your first task is to undo whatever changes you have made, then to Edit the Question to report the spatial reference and envelope metadata of the datasets.
    – Vince
    Aug 15 at 0:41
  • 1
    When you say you manually set the coordinate system of the Shapefile. Did you assign it, or did you "Reproject" it? You need to reun the shapefile through the Reproject tool to convert from one system to another. Aug 15 at 0:41
  • @KeaganAllan I've tried both -- neither quite do the trick. I'll update the post with details.
    – ttoshiro
    Aug 15 at 0:49
  • @Vince Thanks for the response, but where/how can I find spatial reference and envelope metadata -- is it contained in the file? I'm not very familiar, sorry.
    – ttoshiro
    Aug 17 at 22:22
  • @ttoshiro - what version of QGIS are you using? Aug 18 at 0:14

1 Answer 1

2

If you look at the EDGAR data (NetCDF) - you will see that the longitude value range is from 0 to 360 degrees. Looking at the metadata of this dataset using an xarray, I saw the units are degrees east.

My guess why your data is not aligning, is the WGS84 system uses negative values for West and South, and positive values for East and North.

If you reproject your county boundary to WGS84, it will convert the coordinates of your boundary into positive or negative values based on their location.

I downloaded one of the tiles and loaded it into QGIS and found the following (I havent styled the data, but you can see the extent of the data):

enter image description here

You can see the western most extent of the data is 0, and is placed on the Prime Meridian. We need to change the NetCDF data from 0-360 longitude to -180 to 180 longitudes. There are a number of ways to do this. I decided to use Python and the xarray library - because I am comfortable with Python and I found some good posts showing examples:

Converting a NetCDF from 0 to 360 to -180 to 180 via xarray

Using rioxarray to assign spatial reference (EPSG:4326) to NetCDF built from CSV with lat long columns

The code I used:

import xarray as xr
import rioxarray
    
with xr.open_dataset(r'Path to NetCDF.nc') as src:
    ds = src
    ds = ds.rio.write_crs("epsg:4326", inplace=True)
    print (ds.coords['lon'])
    ds.coords['lon'] = (ds.coords['lon'] + 180) % 360 - 180
    ds = ds.sortby(ds.lon)
    print (ds.coords['lon'])

    ds.to_netcdf(r'Path to NetCDF_Test.nc')

The output of this script:

enter image description here

and checking the extents of the layer:

enter image description here

The way this is written, it could be used in a loop and run over the whole dataset / time series. I havent checked the second dataset, but this should get the EDGAR data in a format you can use.

** UPDATE **

I have downloaded the "Counties" shapefile, as well as the FFDAS data. Using the method described in this answer, I converted the 0 to 360 degree EDGAR data and added it to a QGIS map, set to WGS84:

enter image description here

I haven't styled it, but you can see your counties shapefile in the correct location in the dataset generated using the methodology described above.

I downloaded the FFDAS files and loaded one - I did not change the projection, I used it as is.

enter image description here

Again, I haven't styled it correctly, but you can see the data is in the correct location.

enter image description here

  1. The FFDAS looks to already be in WGS84 (or something equivalent).
  2. The EDGAR data needs to be converted from 0 - 360 to -180 to 180.
  3. Your counties data can be used as is, the "on the fly" projection should correctly sort it out. If you want to have it all the same, then you need to set the projetion first to "NAD_1983_Lambert_Conformal_Conic" and then Reproject to WGS84
3
  • Thanks for the response. While this makes the EDGAR emissions map larger and in a different location, it still doesn't overlap with the county shapefile, and neither does the FFDAS emissions map (which is still upside down). I'll update the post with this information.
    – ttoshiro
    Aug 17 at 22:04
  • 1
    @ttoshiro - I have updated my answer - You may need to follow Vince's advice in the very first comment of your question - The data lines up for me. Aug 18 at 1:42
  • Thank you -- this is really interesting. What QGIS version are you using? Are you adding the .nc files as raster layers? Also, I believe I added the spatial reference/envelope metadata to the original post, if that helps.
    – ttoshiro
    Aug 18 at 19:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.