Python Forum
Python Forum

Problem with masking precipitation data into a shapefile
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Thread Modes
Problem with masking precipitation data into a shapefile
atrwmb109
Unladen Swallow
Posts: 1
Threads: 1
Joined: Jul 2019
Reputation: 0
#1
Jul-22-2019, 05:10 PM (This post was last modified: Jul-22-2019, 09:14 PM by Yoriz.)
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! 
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] = 1
This is the full error I have been getting:
Click here to reveal error

(By the way, line 97 points to if Ngl.gc_inout(data_lat[k,l], data_lon[k,l], lats, lons) != 0:)
Find
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  numpy masking/filtering nilamo 3 6,881 Dec-04-2021, 10:28 PM
Last Post: nilamo
  Problem with saving data and loading data to mysql kirito85 4 6,260 Feb-08-2019, 10:53 AM
Last Post: kirito85

Forum Jump:

User Panel Messages

Log Out
My Profile
Pay your profile a visit
User Control Panel
Do some changes on your profile
My Messages
View private messages unread
Avatar
Change avatar
Signature
Change signature
Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020