3

I have created a plot in Python 3.7 over Alaska. Firstly, the pixels covered a lot more than shown in the image (below). Now I have set the lat lon extend to fit the plot 'better'.

  1. Is there a way to optimize this, so the plot actually fits within the border of Alaska?
  2. Can I do it prior to applying it in Python, hence in QGIS?

plot

Here is the Python script

import netCDF4 as nc
import xarray as xr
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

ds = nc.Dataset(path + 'test.nc')

lons = ds.variables["lon"][:]
lats = ds.variables["lat"][:]
prob = ds.variables["test"][:]

lon_0 = -140
lat_0 = 65

m = Basemap(width = 5000000, height = 3500000, resolution = 'l', projection = 'stere', lat_ts = 40, lat_0 = lat_0, lon_0 = lon_0)

lon, lat = np.meshgrid(lons, lats)
xi, yi = m(lon, lat)

# Plot Data
cs = m.pcolor(xi,yi,np.squeeze(prob))

# Add Grid Lines
m.drawparallels(np.arange(-80., 81., 10.), labels=[1,0,0,0], fontsize=10)
m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10)

# Add Coastlines, States, and Country Boundaries
m.drawcoastlines()
m.drawcountries()
2
  • do you want to a) crop away the sea, or b) make the pixels align better to the map (hard). also please include which libraries you have imported, (rasterio, Gdal, numpy) so we can give a more directed answer. Dec 2, 2020 at 14:18
  • Crop away the sea @GevaertJoep, please. I'll add the libraries now, thanks.
    – Thomas
    Dec 2, 2020 at 15:40

1 Answer 1

1

You will need some more libraries - either GDAL or Rasterio. Here is an example from the documentation:

import geopandas
import rasterio
import rasterio.mask

shapefile = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
shapes = [feature["geometry"] for feature in shapefile]
del shapefile

with rasterio.open("tests/data/RGB.byte.tif") as src:
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta
out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

with rasterio.open("RGB.byte.masked.tif", "w", **out_meta) as dest:
    dest.write(out_image)

Where I used geopandas in stead of fiona because it has an earth map in the standard database. I have not tested it, but good luck!

1
  • 1
    It does not completely work. I get the error of __enter__ when using the with parameter in geopandas. And, you are working with a tif file. I need it to be NetCDF, thanks.
    – Thomas
    Dec 3, 2020 at 14:10

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.