Hi there! I am new to this forum but I did read the rules beforehand so if I make any mistakes please let me know! I am currently an undergrad researcher going into my junior year of college and new to Python programming. I have been trying to load in a file with years of precip data, and average over the time dimension to get a 2D field of precipitation with latitudes and longitudes. I plan to mask the precip data into a shapefile, that I found of Long Island online. I have been working with this code for some time, and have had many issues along the way. One issue I had was that once the lats and lons were put into a list, the code was not taking a true min/max of the lats and lons. I attempted to fix this by setting the max/min lons and lats manually to a desired number, but I have been getting this ambiguous error that I do not know how to fix. Any help would be greatly appreciated!
(By the way, line 97 points to
import numpy as np import Ngl import xarray as xr import shapefile import fiona import pyproj f = "/Users/austinreed/Desktop/URECA2019/LongIslandAnnualtotalprecipitationncfilesnew/*.nc" DS1 = xr.open_mfdataset(f,concat_dim="time") #this command opens the file using xarray data_lat = DS1.lat.values[0,:,:] #you must read in the latitude and longitude variables for plotting, but remember they might not be named "lat" and "lon" in your file data_lon = DS1.lon.values[0,:,:] print(data_lat) print(data_lon) #data_lat = y #data_lon = x nlats = len(data_lat[:,0]) nlons = len(data_lon[0,:]) #lon_grid, lat_grid = np.meshgrid(x,y) wgs84 = pyproj.Proj("+init=EPSG:4326") # LatLon with WGS84 datum used by GPS units and Google Earth daymet = pyproj.Proj("+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +ellps=WGS84 +units=m +no_defs ") #xx, yy = pyproj.transform(daymet, wgs84, lon_grid, lat_grid) data = DS1.prcp.values #this is the variable that you want to read in, for example, temperature or precipitation so change "var_name" to the name of your variable in the netcdf file annual_sums=np.add.reduceat(data, np.arange(0, len(data), 365)) time_avg= (np.mean(annual_sums,axis=0)) fp = "/Users/austinreed/Desktop/URECA2019/longisland_counties/longisland_counties.shp" #path to your shapefile mask = np.zeros((nlats,nlons)) for shapes in fiona.open(fp): geom_type = shapes['geometry']['type'] geom_coords = shapes['geometry']['coordinates'] #coordinates should be the latitude and longitudes of Long Island print(np.shape(geom_coords)) if geom_type == 'Polygon': test_array = np.full_like(time_avg,np.nan) #test_array = np.full_like(precip, np.nan) part_squeeze = np.squeeze(geom_coords) part_shape = np.shape(part_squeeze) lons = [] lats = [] for i in range(part_shape[0]): lons.append(part_squeeze[i][0]) lats.append(part_squeeze[i][1]) max_lat = max(lats) min_lat = min(lats) max_lon = max(lons) min_lon = min(lons) for k in range(nlats): for l in range(nlons): if data_lat[k,l] <= max_lat and data_lat[k,l] >= min_lat and data_lon[k,l] <= max_lon and data_lon[k,l] >= min_lon: if Ngl.gc_inout(data_lat[k,l], data_lon[k,l], lats, lons) != 0: test_array[k,l] = time_avg[k,l] mask[k,l] = 1 elif geom_type == 'MultiPolygon': #test_array = np.full_like(precip, np.nan) test_array = np.full_like(time_avg,np.nan) for part in geom_coords: part_squeeze = np.squeeze(part) part_shape = np.shape(part_squeeze) print(part_shape) lons = [] lats = [] for i in range(part_shape[0]): lons.append(part_squeeze[i][0]) lats.append(part_squeeze[i][1]) min_lon = -74.041867 #Left Longitude Limit max_lon = -71.85615 #Right Longitude Limit min_lat = 40.541722 #Lower Latitude Limit max_lat = 41.292377 for k in range(nlats): for l in range(nlons): if data_lat[k,l] <= max_lat and data_lat[k,l] >= min_lat and data_lon[k,l] <= max_lon and data_lon[k,l] >= min_lon: if Ngl.gc_inout(data_lat[k,l], data_lon[k,l], lats, lons) != 0: test_array[k,l] = time_avg[k,l] mask[k,l] = 1This is the full error I have been getting:
(By the way, line 97 points to
if Ngl.gc_inout(data_lat[k,l], data_lon[k,l], lats, lons) != 0:
)