load "./shapefile_utils.ncl"

undef("plot_shapefile")
function plot_shapefile(wks,shapefile_name)
local f, shp_type, sres, geometry_type
begin
;---Open file
  f = addfile(shapefile_name,"r")
;---Check validity of shapefile
  if(.not.isatt(f,"geometry_type")) then
    print("plot_shapefile: Error: This doesn't appear to be a valid shapefile")
    print("No 'geometry_type' attribute is present.")
    exit
  else if(.not.any(f@geometry_type.eq.(/"point","polyline","polygon"/))) then
    print("plot_shapefile: Warning: the 'geometry_type' attribute on '" + shapefile_name + "'") 
    print("             is not set to 'point', 'polyline', or 'polygon'.")
    print("             Will try to plot this data using 'point', but you")
    print("             may want to examine the file further to make sure")
    print("             it is correct.")
    geometry_type = "point"
  else
    geometry_type = f@geometry_type
  end if
  end if

;---Read lat/lon and create map
  lat = f->y
  lon = f->x
  map = create_shapefile_map(wks,shapefile_name,geometry_type,lat,lon)
  
  sres = True
  if(geometry_type.eq."point") then
    sres@gsMarkerIndex = 16
    map@id = gsn_add_shapefile_polymarkers(wks,map,shapefile_name,sres)
  else if(geometry_type.eq."polyline") then
    map@id = gsn_add_shapefile_polylines(wks,map,shapefile_name,sres)
  else if(geometry_type.eq."polygon") then
    sres@gsFillColor = "cyan"
    map@id1 = gsn_add_shapefile_polygons(wks,map,shapefile_name,sres)
    map@id2 = gsn_add_shapefile_polylines(wks,map,shapefile_name,sres)
  end if
  end if
  end if

  return(map)
end

begin
    data1d=asciiread("jma.csv",-1,"float")
    shp_filename="Bicol_region.shp" 
    print_shapefile_info(shp_filename)
    lat = data1d(0::2)
    lon = data1d(1::2)
    printMinMax(lat,0)
    printMinMax(lon,0)

;---Create dummy data array so we can attach lat/lon attributes to it.
    npts = dimsizes(lat)
    data = new(npts,float)
    data@lat1d = lat
    data@lon1d = lon

;---Return a mask array (0s and 1s) indicating which lat/lon pairs fall inside the shapefile region
    opt = True
    opt@return_mask = True
    data_mask = shapefile_mask_data(data,shp_filename,opt)

;---Open device to draw graphics.
    wks = gsn_open_wks("x11","Bicol")

    map = plot_shapefile(wks,shp_filename)

    mkres = True
    mkres@gsMarkerIndex = 16
    mkres@gsMarkerColor = "Red"
    mkres@gsMarkerSizeF = 10.
    id_all = gsn_add_polymarker(wks,map,lon,lat,mkres)

    ii = ind(data_mask.eq.1)
    if(.not.ismissing(ii(0))) then
      nii = dimsizes(ii)
      print("Found " + nii + " lat,lon pairs inside the shapefile regions.")
      print("    " + lat(ii) + " / " + lon(ii))
      mkres@gsMarkerColor = "DarkSlateBlue"
      id_mask = gsn_add_polymarker(wks,map,lon(ii),lat(ii),mkres)
    else
      print("Didn't find any lat/lon values inside the shapefile regions")
    end if

    draw(map)
    frame(wks)
end