texts: Contained Query does not work in spatial coordinates in GeoDjango?

texts:
These are my Models. I imported the shapefiles to this Models

class Parcel(models.Model):
    objectid = models.BigIntegerField()
    apn = models.CharField(max_length = 80)
    lot_size = models.BigIntegerField()
    neibrhc = models.CharField(max_length = 80)
    street_nbr = models.CharField(max_length = 80)
    street_nam = models.CharField(max_length = 80)
    city = models.CharField(max_length = 80)
    zip = models.BigIntegerField()
    tra = models.CharField(max_length = 80)
    subdivisio = models.CharField(max_length = 80)
    lot = models.CharField(max_length = 80)
    landuse = models.CharField(max_length = 80)
    lu_general = models.CharField(max_length = 80)
    lu_specif = models.CharField(max_length = 80)
    lu_detail = models.CharField(max_length = 80)
    lu_use = models.CharField(max_length = 80)
    lu_sec_use = models.CharField(max_length = 80)
    area = models.FloatField()
    perimeter = models.FloatField()
    shape_are = models.FloatField()
    shape_len = models.FloatField()
    geom = models.MultiPolygonField(srid = 3492)

    def __str__(self):
        return 'City: %s' % self.city

class Primary(models.Model):
    id = models.CharField(max_length = 12, primary_key = True)
    address = models.CharField(max_length = 30)
    city = models.CharField(max_length = 30)
    zip = models.CharField(max_length = 10)
    state = models.CharField(max_length = 2)
    unit = models.CharField(max_length = 13)
    primpgon = models.BigIntegerField()
    numpgons = models.BigIntegerField()
    x = models.FloatField()
    y = models.FloatField()
    censusbloc = models.CharField(max_length = 15)
    objectid = models.BigIntegerField()
    geom = models.MultiPolygonField(srid = 4326)

from django.contrib.gis.geos
import Polygon
parcel = Parcel.objects.filter(objectid = "22520").values()
geojson = parcel[0]['geom'].geojson
format_l = ast.literal_eval(geojson)
coordinates_data = format_l['coordinates'][0][0]

poly = Polygon(coordinates_data, srid = 3492)
poly_parcel = poly.transform(4326, clone = True)
polygon = poly_parcel.coords
print(polygon)

I define the Polygon as below

(((-121.49440799279022, 38.55725848782723),
    (-121.49438844924772, 38.557303481514126),
    (-121.4943760310021, 38.5573320694682),
    (-121.49436263531841, 38.557362909896675),
    (-121.49402385986245, 38.557269114460084),
    (-121.49406987333441, 38.55716268909225),
    (-121.49440799279022, 38.55725848782723)), )

My Aim is I want to locate the Primary data within the Parcel Data. So I tried this.

    Primary.objects.filter(geom__contains=polygon).values()
    #and also tried with contained lookup
    Primary.objects.filter(geom__contained=polygon).values()

When I try this querysets it throws the "ValueError: Cannot use object with type tuple for a spatial lookup parameter."

I tried the with List also, but I am getting the same error.


texts:
Sorry if this is basic.  I have an excel sheet with lat/long in it and a field with a radius, in meters, included.  I would like to create a circle or ellipse feature class around these points.  Example:

ObjID      Lat     Long       Radius(m)
0         26.514  -115.253      6
1         33.516  -111.235      54

I can create a point feature class out of this but i do not know how to create a circle or ellipse bases on the "Radius (m)" field. I am using Arc10.0

Thanks!!


texts: ArcGIS API for JavaScript - Create simple buffer or circle graphic against one point

texts:
I am stuck, trying to figure out how to create a simple graphic for adding a circle based on one point in my web map. I need the circle to be centered on my single point, with a radius of 5 miles. I saw examples for the geodesic circle and that didn't work for me, as well a creating a buffer, which I don't think I need. I just need a simple circle graphic. Is there a simple way to do this? I have a map with one point and need this graphic for a boundary. I apologize as I am new to this and trying to learn this API quickly.

thanks


texts: Can layers with different extents in ArcMap 10 be made to line up?

texts:
I am currently trying to get some insect sampling coordinates that are represented by points to project on top of a polygon that shows all the different plots of land in a ranch. The coordinates for the points are in decimal degrees and were collected using iPhones. Both layers have the same exact Projection and Geographic coordinate systems, but do have different extents. They also have the same false easting, false northing, central meridian and datum. Is there a method of getting the points to line up on the polygon by changing extents, or is there another method that I can use?

Here are the extents of both layers if it helps.

For the polygon:

top: 3050776.094191 m
left: 631422.082049 m
bottom: 3031869.963465 m
right: 644361.972775 m

The points:

top: 27.573800 m
left: 97.000000 m
bottom: 27.000000 m
right: 97.640800 m


texts: Can ArcGIS Online Wep Maps support symbols with leader lines for overlapping features?

texts: I'm creating a web map with ArcGIS Online and have some points with the same address. I've seen some web maps where the overlapping points will separate with leader lines when you zoom closely to them (see image). This makes it easier to see that there are multiple features present and also makes it easier to navigate their pop-outs (rather than scrolling through them individually). Can this type of visualization be done with ArcGIS Online?


texts: Cut polygon layer into a new subset by the boundaries of another polygon layer

texts:
ArcGIS 10.2.1

I have a feature class of wetland mapping. It is a series of various polygons, all irregular shapes and sizes.

I also have a parcel boundary feature class that is basically even sized rectangular polygons.

I am trying to divide the wetland polygons into a "subset", where the polygons are "cut" by the boundaries of the parcels, that way I can get area estimates per parcel of how much wetlands each parcel has.


texts: Change decimal commas to points in QGIS created dbf-file

texts:
There is a bug in the QGIS desktop application where two merged shapefiles produce a dbf-file with commas as decimal separators.

I need to have this merged file imported to PostGIS database, but the import using shp2pgsql and psql fails to error ERROR:  invalid input syntax for type numeric: "6817749,945313"

The QGIS version used is the one in Ubuntu 16.04-repositories (2.8.6-Wien)

The bug is a bit old https://issues.qgis.org/issues/8332

How can I update the dbf-file to have the correct decimal separators?


texts: How to use “left” and “strpos” to get the left part of a string with an expression?

texts: I have a table in QGIS with full names and I want to parse into first name(s) and last name. I’ve created the columns and am using Expression Builder to do this. Here's what I'm trying and it's not working. left( "name",strpos( "name", " ") ) with name being the full name. I've tried just putting in strpos( "name", " ") to see if things are working and nothing shows up. I click on Update Selected and nothing changes. I've tried putting the whole expression in last_name = "name",strpos( "name", " ") ) and no luck. QGIS on Mac OSX. PostGIS


texts: Determining active layer regardless of its name using ArcPy?

texts:
ArcPy has a method of defining an active map regardless of its filename or path:

mxd = arcpy.mapping.MapDocument("CURRENT")

I wonder if there is a similar way to define an active layer. I need it so I can use it in a labeling script. I specifically want to avoid hard-coding the layer’s name.


texts: Are large companies successfully using SQL Server GeoSpatial functionality?

texts:
I am working on a project where I have to allow users to find things that are near by.  In the past I have just manually stored lat/longs and done my own calculations but I'm curious if it's worth looking into using the GeoSpatial features in SQL Server 2008.

Are people out there having success using this for their spatial queries?  Are there examples of large companies out there that are using this?  Are there any common problems that I should be aware of?


texts: Error while ploting multiple geopandas (ArcMap and KML) in the same projection?

texts:
I appreciate your help on:

Finding possible errors on projections management
Guide me on which should be the best way to know which "labels - colors" are within the path line  (not exposed here)

--
Source 1: Geopanda imported from ArcMap in .shp format - It contains all shapes (4k) and is used as the basic map. It's very populated in terms of ArcMap columns and attributes
Source 2: Geopanda imported from KML coordinates in .csv format - It's a ROV path in UTM 31N zone.
Steps done:
Import and plot ArcMap representation in GeoPandas. Done! "The same" as ArcMap

Import and plot KML path in GeoPandas. Done!

Make them work in the same CRS
# print(gdf.crs) == > epsg:32631

# Source 2 is directly readed from CSV. Need to assign geocode.
# For a non GIS guy, makes sense to apply || epsg: 32631 or 25831

# print(track.crs) == > 32633   Value determined automatically. As the input for source 2 is just KML coordinates
track = track.to_crs(gdf.crs)

print(gdf.crs)     # After conversion, both are the same: epsg:32631
print(track.crs)

# Follow the convention of plotting both geoPandas by sharing an axes
fig, axes = plt.subplots(figsize=(15,15))
gdf.plot(ax=axes, color=gdf['Color'])
track.plot(ax=axes)

Incorrect output:
Can you see the miniature of SOURCE 1?

I assume this point is my KML path:

What I expected
(picture from ArcMap):


texts: Creating point in shapefile from latitude/longitude using GeoTools?

texts:
I m new to geotools and i find bit difficult to plot a point using the latitude and longitude values. can any one heple me out to solve my problem.

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory():
SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);
Point point = geometryFactory.createPoint(new Coordinate(longitude, latitude));
featureBuilder.add(point);

but this code is not working properly. could any one explain why..?


texts: Postprocessing.py error when running TauDEm on new machine

texts:
When running any operation using the TauDEm provider I get the following errors and no output raster file is produced.

2018-07-06T10:19:55 0   TauDEM execution console output
2018-07-06T10:19:55 2   Error loading result layer:
            Traceback (most recent call last):
              File "C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python/plugins\processing\gui\Postprocessing.py", line 77, in handleAlgorithmResults
                isRaster)
              File "C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python/plugins\processing\tools\dataobjects.py", line 186, in load
                + '\nCheck the processing framework log to look for errors')
            RuntimeError: Could not load layer: C:\Users\DAVID~1.LAV\AppData\Local\Temp\processing43e2c87d54c34dc48591433ce19aeb4e\879e56b7c22445d6881832ccfde7ae20\fel.tif
            Check the processing framework log to look for errors

I have uninstalled qgis and re-installed a fresh install after deleting my qgis2 folder and also everything I could find in the registry. I have installed 2.18.14 as this was the version I had before my computer went up in smoke yesterday.

This new computer is also running Windows 10 as my previous one.


texts: Setting different context menu actions for each layer type with PyQGIS

texts:
I want to set different actions(context menu) for vectors,rasters,etc .

class MyMenuProvider(QgsLayerTreeViewMenuProvider):
def __init__(self, view):
    QgsLayerTreeViewMenuProvider.__init__(self)
    self.view = view

  def createContextMenu(self):
    if not self.view.currentLayer():
      return None
    m = QMenu()
    m.addAction("Show Attribute Table", self.openAttributeTable)
    return m

  def openAttributeTable(self):
    #iface1.openFeatureForm(iface1.activeLayer(), feature, False)
    iface1.showAttributeTable(iface1.activeLayer())

then using it as

 view = iface1.layerTreeView()
 provider = MyMenuProvider(view)
 view.setMenuProvider(provider)

But this sets menu provider for all layers , I want to have different actions for different layers


texts: Determining when ArcGIS feature was created?

texts:
In ArcCatalog I'm looking at a feature's table fields.

The DATECREATED column is empty for all features.

How can one determine when a specific feature was created?


texts: Iterating over several layers and running process using PyQGIS?

texts:
I am trying to find road intersections within each of my memory Polyline layers (which are named road_layer_i for i in (0-99)) so this is what I tried:

import processing

self=qgis.utils
layers = self.iface.mapCanvas().layers()

for layer in layers:
   for i in range(0,100):
      if layer.name()=='road_layer_%s' % (i) :
        processing.runalg('qgis:lineintersections', layer, layer, 'id', 'id', 'C:/Phoenix/intersection_%s.shp' % (i))
        IntersectionLayer=QgsVectorLayer('C:/Phoenix/intersection_%s.shp' % (i) , 'intersection_%s' % (i), "ogr")
        QgsMapLayerRegistry.instance().addMapLayers([IntersectionLayer])

I don't get any error message but each layer created contains points that are the road intersections for the layer 'road_layer_0'.

I really don't understand why it does that so I would be grateful if someone could tell me what I did wrong.

Also, whenever I simply try

import processing

self=qgis.utils
layers = self.iface.mapCanvas().layers()

for layer in layers:
  if layer.name()=='road_layer_33' :
    processing.runalg('qgis:lineintersections', layer, layer, 'id', 'id', 'C:/Phoenix/intersection_33.shp')
    IntersectionLayer=QgsVectorLayer('C:/Phoenix/intersection_33.shp' , 'intersection', "ogr")
    QgsMapLayerRegistry.instance().addMapLayers([IntersectionLayer])

I still get the intersections for the road_layer_0...

After answer by @gcarrillo it is still not working. I indeed obtain 100 layers labbeled interesection_i (i in 0-99) but they all contain only the intersections for road_layer_0.
Here is the code I used to create my road layers :

from fTools.fTools import fToolsPlugin
import shapely
from qgis.core import QgsRasterLayer, QgsCoordinateReferenceSystem, QgsGeometry

RoadLayer=QgsVectorLayer("C:/Phoenix/Project/Phoenix/Phoenix_roads.shp" ,    "Phoenix_roads", "ogr")
RoadLayer.setCrs( QgsCoordinateReferenceSystem(3857, QgsCoordinateReferenceSystem.EpsgCrsId) )
#QgsMapLayerRegistry.instance().addMapLayers([RoadLayer])
#Get the data from this layer
provider = RoadLayer.dataProvider()
fields = RoadLayer.pendingFields()
field_names = [field.name() for field in fields]

# create list of points used to create rectangles
partit=partition(10,c,d)
for i in range(0,100):
#    create new layer in which we will put the roads that are inside the rectangle
   layer =  QgsVectorLayer('LineString', 'road_layer_%s' %(i) , "memory")
   pr = layer.dataProvider()
   attrib_layer=input_layer_attrib_names = RoadLayer.dataProvider().fields()
   oldattributeList = RoadLayer.dataProvider().fields().toList()
#Generate attributes in destination layer
   newattributeList=[]
     for attrib in oldattributeList:
        if layer.fieldNameIndex(attrib.name())==-1:
           newattributeList.append(QgsField(attrib.name(),attrib.type()))
     pr.addAttributes(newattributeList)
  #    add empty attribute fields to new layer
     layer.updateFields()
  #create points left-bottom and right-top of rectangle we want to select
     x=partit[2*i]
     pointA = QgsPoint(x[0],x[1])
     y=partit[2*i+1]
     pointB = QgsPoint(y[0],y[1])
   #   create rectangle
     Rect=QgsRectangle(x[0],x[1],y[0],y[1])
   #   create request
     request = QgsFeatureRequest().setFilterRect(Rect)
     for f in RoadLayer.getFeatures(request):
    #        add geometry from road to new road
         geom=f.geometry()
         coord=geom.asPolyline()
         pt = QgsFeature()
         pt.setGeometry(QgsGeometry.fromPolyline(coord))
    #       add attributes
         atr = dict(zip(field_names, f.attributes()))
         pt.setAttributes(f.attributes())
     #        add new feature to layer
         pr.addFeatures([pt])
         layer.updateExtents()
    #        add layer to map canva
      QgsMapLayerRegistry.instance().addMapLayers([layer])

Is there anything wrong with this method that makes the use of processing.runalg('qgis:lineintersections', ... ) impossible to iterate over the layers thus created ?


texts: Map automation to iterate through attribute table using ArcPy?

texts:
I have a python code which iterates through all attribute fields and creates a map for each field automatically. However, in the code I have to specify each attribute manually which was fine so far. Now I have an attribute table with many fields and I need a way to do the same without specifying the attributes manually (as done with "Attrlist" in the code. Any suggestions?

This is the code I am using so far:

outPath = "D:\\Atlas\\test\MxdMaps\\Mxd_MapAutomation\\"
# list of attributes with the same datatype
attrlist = ["EthnicGroups_totnumberoffarmhh_txt_WNHDeAA01",
            "EthnicGroups_totnumberoffarmhh_txt_WNHDeAA02",
            "EthnicGroups_totnumberoffarmhh_txt_WNHDeAA03"]
# layer which has to change the attribute field
featurelayer = arcpy.mapping.ListLayers(mxd, "EthnicGroups_totnroffarmhh")[0]
# text element that has to change (e.g. title...)
title = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "map_title")[0]
for attr in attrlist:
        #if featurelayer.symbologyType == "GRADUATED_COLORS":
                #print attr
                # change value field
                featurelayer.symbology.valueField = attr
                #featurelayer.symbology.addAllValues()
                featurelayer.symbology.reclassify()
                #print "Break values: %s" % featurelayer.symbology.classBreakValues
                # change tilte
                title.text = attr
                # update view
                arcpy.RefreshActiveView()
                arcpy.RefreshTOC()
                # export as jpeg
                outJPEG = os.path.join(outPath, attr + ".jpg")
                arcpy.mapping.ExportToJPEG(mxd, outJPEG, "PAGE_LAYOUT", 200, 200, 300, False, "24-BIT_TRUE_COLOR", 90, False)
                # export as AdobeIllustrator .ai
                #outAI = os.path.join(outPath, attr +  ".ai")
                #arcpy.mapping.ExportToAI(mxd, outAI, "PAGE_LAYOUT", 200, 200, 300, "BEST", "RGB", "VECTORIZE_BITMAP", False)
del mxd


texts: Counting polygons intersecting with each others in the same Vector layers

texts:
I am currently working on some line string data that are very close by and have the same Start and endpoints for the most part. I want to create midpoints or average centroids that will tell me how many of those are grouped up close by.

Using mainly SQL queries i have been able to confirm that the centroids are close enough to intersect with each other using a 15 meter buffer, it is functionnal so far.

SELECT

feeder_NRO02.ca_pk_uid AS ID,
ST_centroid(feeder_NRO02.geom) AS geom,
1 AS count

FROM
feeder_NRO02

My issue is that so far, I haven't been able to count how many times those buffers intersect with each other or create an average centroid.

Any advice? I have QGIS and some basic database tools at my disposal for this.


texts: Copying code generated when you put data into GUI via QGIS Plugins?

texts:
In QGIS 2.18.6, I am using the v.distance tool based on two point layers.

Layers are:

TestPhotos
TestObjects

The inputs are:

I have a Python file called Test.py, located C:\OSGeo4W64

I want to use this file to run v.distance with these two layers as per the abovementioned input (in the image). The output is two layers (a line and point layer). The points are the photos and the lines are from photos to the object.

All I want is that data (the line and point shapefiles).

I was thinking I could just use OSGeo4W.bat or the QGIS Python console to run the file and achieve the same result as I would if I used the GUI shown above.

I figure that it would look something like this site:

v.distance [-pa] from=name [from_layer=string] [from_type=string[,string,...]] to=name [to_layer=string] [to_type=string[,string,...]] [output=name] [dmax=float] [dmin=float] upload=string[,string,...] [column=name[,name,...]] [to_column=name] [table=name] [separator=character] [--overwrite] [--help] [--verbose] [--quiet] [--ui]

So when customised for my need, it would be like this:

v.distance from="TestPhotos" from_type="point,line,area" to="TestObjects" to_type="point,line,area" dmax="-1" dmin="-1" upload="to_attr" column="OBJ_ID" to_column="OBJ_ID"

This is what I have so far in Test.py:

import sys
sys.path.append("C:/OSGeo4W64")
import qgis.core

v.distance from="TestPhotos" from_type="point,line,area" to="TestObjects" to_type="point,line,area" dmax="-1" dmin="-1" upload="to_attr" column="OBJ_ID" to_column="OBJ_ID"

How would I write the python file Test.py so I can run it with OSGeo4W.bat and get the same result as I would with using the GUI?


texts: Point shapefile shows in map but not in print composer or on export

texts: I have a point shapefile in QGIS 2.18.  I can see the points on the map but they do not appear in the print composer.  The points do not show on export to pdf.  I do have other points that display just fine and no point files have scale dependent display settings.  The points in question are symbolized with categorized SVGs.  Why will my points not show?  Note, the shapefile that is visible was created in ArcGIS and the points that are not visible were created in QGIS.


texts: How to define a point by coordinates in QGIS (Plugin)

texts:
I am developing a QGIS plugin and got currently stuck at setting the geometry of a feature based on coordinates.

I retrieve coordinates through an SQL query, create a feature and try to set the geometry.

#fetch sql query
rec = dict_cur.fetchone()

#create feature
feature = QgsFeature()
feature.setGeometry(QgsGeometry.fromPoint(QgsPoint(rec[1], rec[2])))

In the code above rec[1] and rec[2] holding the information about longitude and latitude. When requesting the type of rec[1] or rec[2] it comes back as <class 'decimal.Decimal'> and I can cast it to float. However, that doesn't change anything.

When I try to set the coordinates for QgsPoint then I am getting the following error message:

feature.setGeometry(QgsGeometry.fromPoint(QgsPoint(rec[1], rec[2])))
TypeError: 'NoneType' object has no attribute '__getitem__'

The error is directely related to QgsPoint. When I seperate QgsPoint and try to provide the coordinates I am getting the same error message.

That is for my plugin. If I use the same code in the Python console of QGIS then it works.

Python is in version 2.7.10
QGIS version is 2.8.2-Wien


texts:
I am working on an OpenLayers 2.10 map using Internet Explorer 8 (used by work), and am encountering some issues with the pan/navigation buttons disappearing after sepcific events:

If I change the layer to another background, the pan/navigation panel buttons lose their picture
If I resize the browser window, the pan/navigation panel buttons lose their picture

When one of these events occurs, the buttons show the broken image link icon.  If I refresh the page, the buttons come back.

After perusing some answers, the only similar solutions I came up with were at OpenLayers - redrawing map after container resize using the map.updatesize command.  These commands did not seem to work, though I am not entirely sure where to place them.  Also, I just started using OpenLayers a couple months ago, and a solution may exist using terminology I am not aware of.

I would like a solution that keeps the pan/navigation/zoom panel buttons constantly displayed, even if layers are changed or if the page is resized.  It is also acceptable for the buttons to refresh after an event that will cause the button to no longer load.  Thanks in advance for any guidance.


texts: Avoid OpenStreetMap loading shops, monuments, restaurants and any other POIs?

texts:
I am loading OpenStreetMap using Leaflet.

More or less I'm doing this:

var objmap = new L.Map(p_objCapaMapa, mapOptions);
var osmUrl = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png';
capaMapa= new L.TileLayer(osmUrl);
objmap.addLayer(capaMapa);

Very easy and it works. But the map is loading a lot of Points of Interest like restaurants, shops, etc A lot of little icons over the tiles of the map.

How can I avoid the map load all these points?

I need Leaflet to just load the map with its roads, squares, streets, parks, etc.

In JavaScript Google Maps API there is a way of not loading these POIs that comes with the map but I do not know a way to avoid POIs over OpenStreetMap using Leaflet.

We have searched a lot for the web but we still do not have a solution to this.


texts: GeoServer GetLegendGraphic without x-mark

texts:
I'm using the GeoServer WMS to display some raster data with a custom color ramp. When I request a legend using GetLegendGraphic, I get a legend entry with a big red X in it:

Is there any way I can get a legend graphic with this X omitted? Here's my SLD:

  <sld:RasterSymbolizer>
    <sld:ColorMap>
      <sld:ColorMapEntry color="#ffffe5" opacity="0" quantity="0"/>
      <sld:ColorMapEntry color="#ffffe5" opacity="1" quantity="0"   label="0 mm"  />
      <sld:ColorMapEntry color="#f7fcb9" opacity="1" quantity="25"  label="25 mm" />
      <sld:ColorMapEntry color="#d9f0a3" opacity="1" quantity="50"  label="50 mm" />
      <sld:ColorMapEntry color="#addd8e" opacity="1" quantity="100" label="100 mm"/>
      <sld:ColorMapEntry color="#78c679" opacity="1" quantity="150" label="150 mm"/>
      <sld:ColorMapEntry color="#41ab5d" opacity="1" quantity="200" label="200 mm"/>
      <sld:ColorMapEntry color="#238443" opacity="1" quantity="250" label="250 mm"/>
      <sld:ColorMapEntry color="#006837" opacity="1" quantity="300" label="300 mm"/>
      <sld:ColorMapEntry color="#004529" opacity="1" quantity="400" label="400 mm"/>
    </sld:ColorMap>
    <sld:ContrastEnhancement/>
  </sld:RasterSymbolizer>

(The opacity=0 ColorMapEntry is there to prevent NODATA pixels from being displayed.)


texts:
So, I'm working with a dataset of polygons generated at different levels of sensitivity to error from a citizen science project. I have been trying to visualize how my sensitivity threshold changes results, but in plotting have noticed odd behaviour in ggplot2, inserting polygons where there should be none. The comparison between plotting the native SpatialPolygonsDataSet and the tidy generated data frame is particularly telling. What's going on here, and how can I fix it?

I start with reading in the data and doing some light processing

library(tidyverse)
library(broom)
library(sp)
library(spdplyr)
library(ggplot2)
library(rgdal)
library(raster)
library(rgeos)
library(lubridate)

#read in spatialPolygonsDataFrame (shapefile)
#from https://www.nceas.ucsb.edu/~byrnes/floating_forests/ff_polys_proj.rds
classifications <- readRDS("./ff_polys_proj.rds")

classifications <- classifications %>%
  mutate(scene_timestamp = parse_date_time(scene_timestamp, orders="ymdHMS"),
         quarter = quarter(scene_timestamp, with_year=TRUE, fiscal_start=11))

(note, this all takes some time as the file is big-ish)

Here's a simple plot of the result focusing on one area filtering down to those polygons with a threshold of 6.

#plot the north coast at a threshold of 6
plot(classifications %>% filter(threshold==6),
     xlim=c(-124.001673, -122.772273),
     ylim=c(38.115699,39.176660) )

OK, now, here's a data frame instead

#convert to data frame
classifications_df <- classifications %>% tidy
classifications@data$id <- rownames(classifications@data)
classifications_df <- left_join(classifications_df, classifications@data)

which can then be ggplotted

##ggplot the north coast
ggplot() +
  geom_polygon(data=classifications_df %>% filter(threshold==6),
               mapping=aes(x=long, y=lat, group=group)) +
  theme_bw(base_size=14) +
  xlab("") + ylab("") +
  xlim(-124.001673, -122.772273) +
  ylim(38.115699,39.176660)

And now... why are there extra polygons where formerly there were none?


texts: QGIS 1.8 fails to rotate point geometries with rule-based style

texts:
I have to do some visualizations in QGIS (1.8.0) and up to now I've been very happy with the new Dialog "Layer Properties -> Style", using the rotation-field to rotate Point-Geometries visualized by a triangular Symbol of the section "Simple Markers".

For me QGIS fails to rotate Point-Geometries with rule-based style. Here what I did:

I've got a Point Layer with a Attribute table with this structure:

ID ; Place   ; orientation ; typ        ; centrality
1  ; Basel   ; NULL        ; center     ; 2
2  ; Liestal ; 120         ; depending  ; 4
3  ; Munzach ; 123         ; depending  ; 5

Under the "Categories"-Section I made two categories: 'center' and 'depending'. The Symbol of the depending Place is rotated, using the 'orientation'-column as the rotation field. Now it shows on the center place.

Now I'm trying to use the "Rotation field" in the "Rule-based"-Section which doesn't work for me. I Want to scale the Place according to their centrality.

In the 'rule-based' section I added a rule with the Filter typ='depending' and I disabled the Symbol. then I did 'refine current rule' -> 'add categories to rule'. Here I classified the Column 'centrality' and added my triangular Symbol in different sizes according to the centrality of the place. This works. Now I want those triangles to show to the center place using the -> Advanced -> rotation field -> 'orientation'.

Unfortunately QGIS doesn't apply the rotation field in Layers with 'ruled-based' Symbology.

Any ideas?


texts: How enable MrSID in Geoserver 2.9.2 [Oracle Linux Server]

texts:
On same server there are 3 Tomcats and each Tomcat has a geoserver application running on it. For some reason 2 geoservers pickup MrSID and other extensions while 1 goeserver does not. Where should I look in the configurations to enable MrSID and other extensions. All geoservers have same .jars in WEB-INF/lib.

EDIT: After geoserver restart MrSID and other coverage formats appeard in the store creation window. I have done nothing to fix it. Miracle!

Note: I have no access to the server logs so won't be asking for them anymore (already asked couple of times!!!) therefore can not give anyone more info about this issue. Also I am certain that this geoserver was restarted several times in the past, so maybe after previous restart something went wrong...


texts: ArcGIS Javascript API - Return 'normal' map extents when panning around the world

texts:
ArcGIS Javascript API maps allow the user to pan around the world, which is a great feature, however I am performing some processes that rely on map extents.  I have discovered that despite the defined projection (wkid), extents are returned in meter coordinates with an origin of lat/lon 0,0 (off the coast of Africa).  In order to get my extents in the projection I need, I use a Geometry service to project them, however when the user pans around the world, the coordinates do not reset, they just keep adding up, resulting in 'NaN' extents returned.

I saw this post, which would be one way to deal with this issue, however I would really like to keep the ability to pan around the world.  Does anyone have any ideas on how to accomplish this?


texts: Atlas - How to set in rows from table linked to polygonlayer?

texts: I have polygonlayer and a textlayer. The polygons have unique ID´numbers. The textfile also have unique ID-numbers (and therefore possible to join together, yes), but it contains several rows with the same ID-number. In Print Composers Atlas it´s no problem to make an output print for each polygon by the ID, but I want to show the rows in the table/text-document as well. I can´t figure out if I´ll have to join the text-document by the polygonlayer or try with the Expression filter?  For Q-GIS 2.0. This i s what I want as output for each output:


texts:
I am using geoserver to serve my shapefile data, and I use SLD as the style solution.

There are two fields in the data related to styling: type and status. The type field contains values like highway freeway parkway and etc. The status filed contains values like planning building completed.

Different arrangement of the type and status should result in different style, as shown in this table:

There are 9 rules for the style. Furthermore, the zoom level should be considered while styling, say I have 4 different rules related to the zoom level, there will be 9*4  = 36 rules finally.

Once there are new type values added, the number of rules will increase too.

I wonder how do you handle this kind of situation ?


texts: Using QField Relation Reference for Attribute Table?

texts:
I would like to relate an attribute class to a point class by relation reference. In QField this is supported, but I can't open the attribute class yet. By selecting the point in the graphic I get the attributes of the point class with the attribute of the relation reference. But I didn't find a way to open the attribute class via the relation reference. Only in edit mode can I add a new entry in the attribute class via the relation reference. But once the record in the attribute class is added, there doesn't seem to be any way to edit the record again. Only a new one can be generated.

Anybody have an idea how it works?


texts: Is there an open source tool for both road network analysis and geocoding?

texts: Is there any one open source software package that can perform both road network analysis and geocoding?


texts: Retrieving and comparing locations in miles

texts:
I am currently storing the locations with a SRID of 4326 (what seems to be most commonly suggested for lat, long).

I need to be able to run a query in postgres to find all locations that are within X miles.

Here's my current query:

SELECT * FROM users WHERE ST_Distance(location, ST_SetSRID(ST_Point(40, -70), 4326)) < 200;

From what I understand, this finds all points (location column in users) within 200 meters of New York (40, -70).

Would a standard conversion from meters to miles after the query or is there a better/more correct way to do this with miles?


texts: How to add points to line at every polygon/polyline intersection and have attribute table including distance of new line sements

texts:
I am working with polyline shapefile.

I need to add point data at every intersection of polyline/polygon which was no problem->simple intersection.

Problem is, points were not added to the ends (start/end) of each polyline. I am working on a basic license and am very limited to what tools are available. Ideally I would use linear referencing tools, or geometric network but neither one is available to me. I will also need the distance of each new (line segment), I do not want to break up the lines into actual line segments if possible.

Suggestions?


texts:
I have some ESRI grid data that I would like to import into Matlab.
As a typical ESRI grid data I have the two folders with the *.adf files and *.dat, .nit and.dat files.

I would like to know how can I import the information stored in those files into Matlab.

Tank you all!


texts:
I have a layer of a csv file with around 700k rows and a polygon layer, both are in wgs84.

I'm trying with the spatial query to calculate what point are in the polygons.

The problem is that the spatial query doesn't count all the points when it's finished, as in the attached picture.

It should mark all points within the polygon.


texts:
I have a spatial layer (shapefile) that has one row for each of one of the attributes in the table (representing management units). I also have a non-spatial CSV file that has wildlife inventory data for each of these management units. The problem is that there are multiple lines of inventory data (one for each year) for the management units in the CSV file, so when I do a join in QGIS only one line of inventory data is joined with the spatial layer and the others are ignored.

Is this possible in QGIS?

I think in ArcMap this is referred to as a one-to-many joins.


texts: How to create a live link from Consultation Manager to ArcGIS Enterprise?

texts:
I am working on a project that views their spatial information through ArcGIS Enterprise. They collect information in Consultation Manager and were hoping they would be able to display the Consultation Manager information on a web map in ArcGIS Enterprise.
Has anyone created a live link between Consultation Manager and ESRI products before? And if so do you have any tips for getting this to work?


texts: Are There Any OpenSource GIS Software Programs With an "Interactive" Rotate Feature Option?

texts:
![Rotate feature option][1]Sorry if the question is too simple, but are there any open source GIS software programs available with an interactive "Rotate vector feature"  option?

Eg : ArcGIS has it in Editor\Rotate.  Manifold has it in Transform toolbar (but its only angle based, not interative) .

Any help is appreciated.


texts: Generate own XYZ Tile (MBTiles) - which URL?

texts:
I want to create my own xyz tile for an own basemap on my webspace.
I want to use the own basemap for a qgis2web leaflet map.

This I have done:

processing toolbox: generate XYZ tiles (MBTiles)
I get this file: OUTPUT_FILE.mbtiles
created a new XYZ connection
the mbtiles file uploaded to
https://........../qgis2web/OUTPUT_FILE.mbtiles

Question:
Which URL for the XYZ connection is the right one?

QGIS 3.10.5 / Windows 10


texts: Select by Location and Join the Attributes of Two Line layers

texts:
I want to import the name of streets from OSM shapefile layer to another line layer based on neighborhood distance. The distances between two layers (line layers) are variable between 0 to 10 meter. The segments are not matched.

I already tried the Select by attribute tool in qgis under vector>data-management tools, but the result is not complete, and setting the parameter is tricky. I used overlap and set the distance differently but didn't t work really. for example, if I set the buffer distance higher I expect to get more matching and therefore more attributes but it is not so. Is there any solution?

Images:
1- The OSM shapefile with the name of streets
2- The Line layer which I want to import the names from OSM layer
3- The result


texts: GDAL reads NEX Climate projection datasets as netCDF on Mac, HDF5 on Ubuntu

texts:
I'm trying to work with the climate projection data that NASA released, and is hosted publicly on S3: http://aws.amazon.com/nasa/nex/. The data is available at s3://nasanex/NEX-DCP30.

The code I'm working with uses the rasterio python library to export band information.

I worked out the pipeline on my Macbook, which had the netcdf package installed by homebrew, GDAL built from source (trunk branch) and configured with python bindings, and rasterio installed with pip.

I am able to work with the netCDF files fine on my Macbook. For instance, if I pull down the file at s3://nasanex/NEX-DCP30/BCSD/historical/mon/atmos/tasmax/r1i1p1/v1.0/CONUS/tasmax_amon_BCSD_historical_r1i1p1_CONUS_inmcm4_200501-200512.nc, and run gdalinfo it, the first 3 lines are this:

Driver: netCDF/Network Common Data Format
Files: tasmax_amon_BCSD_historical_r1i1p1_CONUS_inmcm4_200501-200512.nc
Size is 7025, 3105

Great.

I can do the same setup on an Ubuntu machine, using GDAL trunk or the 1.10.1 tag, and have HDF5 and netCDF enabled (as seen by this output in GDAL's ./configure script:

HDF5 support:              yes
NetCDF support:            yes

However, when I do a gdalinfo on the same file, I get a different driver being used:

Driver: HDF5/Hierarchical Data Format Release 5
Files: tasmax_amon_BCSD_historical_r1i1p1_CONUS_inmcm4_200501-200512.nc
Size is 512, 512

I've tried this on a Ubuntu VM on my Macbook, an Ubuntu AMI on EC2, and a workstation that has Ubuntu as it's native OS. All the same.

Worse, when I try to open the file with rasterio, while on my Macbook everything runs great, on Ubuntu I get a a segmentation fault when using rasterio.open(...)

Segmentation fault (core dumped)

This only happens with the NEX dataset netCDF files. I can read another netCDF file I got from a different source just fine...gdalinfo produces the correct Driver for that netCDF file:

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: BCSD_0.5deg_tas_Amon_access1-0_rcp85_r1i1p1_200601-210012.nc
Size is 720, 278

Does anyone have any insight into why these netCDF files would be readable on my Macbook and unreadable (and considered to be corrupt HDF5 datasets) on Ubuntu?


texts: Automatically apply blending mode to raster

texts:
Currently using QGIS 2.18.24, but tested this on 2.14.18 as well.

I'm trying to load a QML file for a raster with a specific blending mode ('multiply'), however QGIS only applies it until I click 'Apply' or 'Okay', after which point it resets back to 'normal'.

I have also tested out the QgsRasterLayer().setBlendMode(6) within python, to no avail.

Is this a known bug, or am I just doing something incorrectly?

Edit: Here is some sample code:

# Apply style layer and set blending mode

layer = QgsRasterLayer(lr)
QgsMapLayerRegistry.instance().addMapLayer(layer)
layer.loadNamedStyle("Path\To\Style.qml")
layer.setBlendMode(6)
layer.triggerRepaint()


texts:
I am trying to generate an Image Collection that corresponds to a long-term 20 year daily average and apply the historical 20 year average to future daily climate projections in the collection NEX-GDDP which runs from 1950-2100.
I have figured out how to apply the 20 year average on a yearly basis, but since I am hoping to apply this average from 2018-2099 in NEX-GDDP, I can't figure out how to iterate or loop this to apply it to  each corresponding day each year.
Here is my working code:
var bandOfInterest = 'tasmin';
var NexFut = ee.ImageCollection('NASA/NEX-GDDP').map(table_bounds)
// Read in the NEX-GDDP HISTORICAL BASELINE

var nex = ee.ImageCollection('NASA/NEX-GDDP').map(table_bounds)
  .select(bandOfInterest).filterMetadata('scenario', 'equals', 'historical')
  .filter(ee.Filter.date('1980-01-01','2000-01-01'))//long term 20 year average
  //
  ;

///////// Calculate average of all model outputs for each day /////////
var days = ee.List.sequence(1,365);

// Group by day, and then reduce within groups by day();
// the result is an ImageCollection with one image for each
// day.
var byDay_20y = ee.ImageCollection.fromImages(
      days.map(function (d) {
        return nex.filter(ee.Filter.calendarRange(d, d,'day_of_year'))
                    .mean()
                    .set('day_of_year', d);
}));
//print(byDay_30y);
var day_image_20y = byDay_20y.toBands();
print(day_image_20y);

var NEX_2018 = NexFut.select('tasmin').map(table_bounds);

var startDate = ee.Date('2018-01-01');
var endDate = ee.Date('2019-01-01');

// Define a base collection.
var pcp = NEX_2018.filterMetadata('scenario', 'equals', 'rcp45').filterDate(startDate, endDate);

// Calc number of days between start and end dates.
var nDays = ee.Number(endDate.difference(startDate, 'day'));

// Enumerate a list of days.
var daysList = ee.List.sequence(0, nDays.subtract(1));

// Map over the list of days to make daily mean composites.
var NEXDayList = daysList.map(function(dayN){
  var t1 = startDate.advance(ee.Number(dayN), 'day');
  var t2 = t1.advance(1, 'day');
  var dailySum = pcp.filterDate(t1, t2).mean();
  var nBands = dailySum.bandNames().size();
  return dailySum

});

print(NEXDayList)
//var images_2019 = pcpDayList.toBands();
//print(images_2019);

// computed anomaly images in the new collection.
var anom = NEX_2018.map(function(image){
    return image.subtract(day_image_20y).set('system:time_start', image.get('system:time_start'));
});
//print(anom.size());


texts: After osm to pgsql conversion with ogr2ogr, can't add/view tables in qgis?

texts:
I'm converting an osm .pbf file to a postgres/PostGIS database with ogr2ogr. I can view the tables in psql, but when I connect to the database in qgis,  I can't see any of the tables after connecting to the database. I'm using the standard osmconf.ini. What am I missing here?

ogr2ogr command:

ogr2ogr -progress -overwrite --config OSM_CONFIG_FILE osmconf.ini -skipfailures -f "PostgreSQL" PG:"dbname=gdal_osm" osmdataset.pbf

In this image you can see that the tables 'multilinestrings', 'multipolygons', etc. are visible in psql, but not when viewed in qgis.


texts: How to remove the layer name in a legend in CARTO?

texts: Is it possible to remove the layer name within the legend (bezirke & berlin_pts). I tried all the obvious ways, but could not manage it.


texts: Batch Export feature attribute to ASCII in ArcGIS for Desktop?

texts: I have around 160 feature classes that I need to convert into ascii files. Spatial statistics tool ---> Utilities ---> Export feature attribute to ascii only allows converting one file at a time. Is there any way to export all my features to ascii in one go? Any script written?


texts: Layers not appearing when updated with UpdateLayer in arcpy

texts:
I have a script that calculates kernel density for several sets of points, adds the output to a base map, updates the added layer to use the same color scheme and resampling function of a template layer, then zooms to the extent of the new layer.

This works fine for some of the resulting maps; the density output shows on the map, the color is as it should be, and it is sampled using bilinear interpolation as I wanted. The problem is that most of the resulting maps do not show the kernel output at all; only the basemap is visible.

I have isolated the code I am using and have altered it to compute one simple test case. Here is what I have:

SubjectLyr = arcpy.mapping.Layer(r"F:\\testFiles\\output\\Population23_KDens.lyr")
mxd = arcpy.mapping.MapDocument("F:\\testFiles\\testMap.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
wPath = "F:\\testFiles\\output\\" # same as in line 1

#Add the layer to basemap
arcpy.mapping.AddLayer(df, SubjectLyr,"AUTO_ARRANGE")

#Update the layer
SubjectLyr = arcpy.mapping.ListLayers(mxd, SubjectLyr.name, df)[0] # change reference to layer in mxd
sourceLayer = arcpy.mapping.Layer(wPath+"templates\\template.lyr") # get template layer
arcpy.mapping.UpdateLayer(df, SubjectLyr, sourceLayer, True)

###Zoom to selection
df.zoomToSelectedFeatures()
lyrExtent = SubjectLyr.getSelectedExtent()
df.extent = lyrExtent
arcpy.mapping.ExportToPDF(mxd, r"F:\\testFiles\\output\\Population23_KDens_Map.pdf")
print('DONE!')

I have read that there is some existing issue with UpdateLayer in ArcMap 10.1, though I am using 10.0. As shown, my code is practically identical to the code that is a "workaround" to this problem described in the above linked forum. This code works for some, but not all of my datasets.

I have also tried using the ApplySymbologyFromLayer_management() function, which does make all the kdens output layers appera; however, this does not change the resampling function used when drawing the output, and I need to to use bilinear interpolation.

Does anyone know what could be going on? Are there any other functions I could use?


texts: List layer names from Map service using python

texts:
I am trying to get a list of layer names within our rest services for all services and folders. I have two codes that I am trying to put together and for some reason the layer code doesn't seen to work within the Map service code. Any suggestions would be great.
I am running on a 10.2 server.

import urllib2, json

server = "hcs084mcgida901"
port = "6080"
baseUrl = "http://gisdev.mcgi.state.mi.us/arcgis/rest/services".format(server, port)

def serviceLayers(URL):
    fURL = URL + "/" + "?f=json"
    openURL = urllib2.urlopen(fURL, '').read()
    outJson = json.loads(openURL)
    return outJson["services"]

if __name__ == "__main__":
    URL = "http://gisdev.mcgi.state.mi.us/arcgis/rest/services".format(folder)

def getCatalog():
    catalog = json.load(urllib2.urlopen(baseUrl + "/" + "?f=json"))
    print 'ROOT'
    if "error" in catalog: return
    services = catalog['services']
    for service in services:
    response = json.load(urllib2.urlopen(baseUrl + '/' + service['name'] + '/' + service['type'] + "?f=json"))
    print '  %s %s (%s)' % (service['name'], service['type'], 'ERROR' if "error" in response else 'SUCCESS')
    folders = catalog['folders']
    for folderName in folders:
        catalog = json.load(urllib2.urlopen(baseUrl + "/" + folderName + "?f=json"))
        print folderName
        if "error" in catalog: return
        services = catalog['services']
    for service in services:
        response = json.load(urllib2.urlopen(baseUrl + '/' + service['name'] + '/' + service['type'] + "?f=json"))
        services = serviceLayers(URL)
    for f in services:
       print '  %s %s (%s)' % (service['name'], service['type'], 'ERROR' if "error" in response else 'SUCCESS'), "Layer : {}".format(f["name"])

getCatalog()


texts:
With the help from this post Lines on reprojected sp objects with Mollweide projection, I draw a world map using Mollweide projection.

I also have some points with latitudes and longitudes, and I want to put those points on the map. What should I do to add these points?

library(sp)
library(maptools)
data(wrld_simpl)
plot(spTransform(wrld_simpl, CRS = ("+proj=moll +lon_0=0 +x_0=0 +y_0=0")))
latitudes <- c(-30, 0, 20, 50)
longitudes <- c(-20, 10, 50, 80)
# how do I plot these four points on the same map?


texts: Creating a radius of influence around points

texts:
I have a vector file of points. What I need to do is create a circle with a radius of 12km round each point.
Would it be possible to then find out what percentage of my total map is covered with circles and what area is uncovered?


texts:
I have a code for NDVI time series from an image collection. This code prints the time series in a table or a chart.
Is there a way to access the values in the table of the time series in google earth engine code? I want to calculate the slope\area between points in the time series in google earth engine enviorment.
Code snippit:
    //Export NDVI Time Series into a Scatter Chart and a CSV format
var TimeSeries = Chart.image.seriesByRegion(
S2_withNDVI, fc, ee.Reducer.mean(), 'NDVI', 3, 'system:time_start', fc_ID)
    .setChartType('ScatterChart')
    .setOptions({
      title: 'NDVI Time Series For Shift 1-2 Argentina',
      vAxis: {title: 'NDVI'},
      hAxis: {title: 'Date'},
      lineWidth: 1,
      pointSize: 4,
      series: {
}});
// Display.
print(TimeSeries);


texts: Has anyone developed a docker container for OSGeo?  Currently, the bootable image can be downloaded from here.  It is set up to run in a virtual machine such as Virtual Box, but I am looking to run it through Docker.


texts: Export Selected Features to a new Feature Class NOT VISIBLE ArcGIS Pro 2.7.2

texts:
I am exporting selected features from a feature class and creating a new feature class within the project geodatabase. The export runs and the feature class is created but the features are not visible. The attribute table is correct, it shows there are features in the file but they do not display. I have tried this with multiple feature classes and I am getting a similar result. I even tried exporting to shapefile and got the same result. I am using ArcGIS Pro 2.7.2.
Has anyone else run into this and if so is there a solution?
I've tried recalculating the spatial index, deleting and recreating spatial index, projecting the file. Nothing has worked.


texts:
This may seem like a simple question, but I cannot seem to just get the part I want from the map layer in the layout tab in ArcGIS pro.

This is a view of what it looks like,

I just want the part where the kernel density layer is plus a part of that basemap.

I try to resize it but then it does not achieve what I want.


texts: Where is the network message log tab in QGIS on Windows7?

texts:
I want to find out what web requests is QGIS exactly sending to my web server when it connects to a WMS and WFS layer.

In my Linux installation of QGIS 2.0 I can see three tabs in the Log Messages window: General, Plugins, Network.
In my Windows installation of QGIS 2.0 I can only see two tabs in the Log Messages window: General, Plugins.

How do I activate the "Network" log messages tab in QGIS 2.0 on Windows 7?


texts: Creating weighted lines on map in CartoDB

texts:
I am trying to illustrate global trade flows within CartoDb. Flows are visualized by a line connecting two lat/long points on a map. The size of the flow would be represented by the width of the line. I have a lat/long and trade magnitude file but cannot figure out, 1) how to get the points to connect in CartoDb using PostGIS functions and 2) how to adjust the line width based on the value using CartoCSS.

Pasted below is the head of the data:

id  start_lat   start_lon   end_lat     end_lon         Value
1   32.89595056 -97.0372    35.04022222 -106.6091944    100
2   41.979595   -87.9044641 30.19453278 -97.66987194    500
3   32.89595056 -97.0372    41.93887417 -72.68322833    750
4   18.43941667 -66.0018333 41.93887417 -72.68322833    1200

Here's an example of the type of visualization I am trying to make. More or less it is a spatial sankey diagram:

from 2012 GEA Global Energy Assessment


texts: Calculating all distances from single points to multiple polygons

texts:
I have two layers: a point layer defining properties (95 objects), and a polygon layer defining plots of land associated with the properties (211 objects). Each point relates to one or more polygons. What I would like to produce is a table of distances of each point to ALL of the polygons associated with that point (nearest edge or centroid of polygon -- either will do). Calculating the nearest distance is relatively easy in QGIS and ArcGIS, but those calculations omit all distant polygons, at least by the methods I have been employing. Ideally I would like an output of

Point_ID | distance to polygon 1 | distance to polygon 2 | distance to polygon 3 | etc.

Any pointers in ArcGIS 10 or QGIS 2.2+ would be most helpful.


texts: How to load GeoJSON into PostgreSQL database using Python?


texts: Web AppBuilder for ArcGIS (Developer Edition): Custom Themes Using Demo Themes

texts:
How to I navigate to the ~/client/stemapp/themes to start creating custom themes, widgets and controllers?

I have been following the protocol (numbered list below) that highlights on their site along with a youtube video, but I have yet to figure out how to enable the environment that will allow me to start editing source code and customizing as shown in the snippet below.

Go to ~/client/stemapp/themes, where the DemoTheme folder is
    located.
Make a copy of the DemoTheme folder in the same directory, and
    rename it NewTheme.
Update the following properties in the manifest.json file: Name:NewTheme

     -Note: This value must be the same as the name of the theme folder.

     -Label: New Theme

Update the localization files in the nls folder as needed. Make sure
    the _themeLabel property in the strings.js file has been updated to
    the name of the new theme.


texts: QGIS curved labels don't show up in print composer

texts:
I use QGIS 2.14.2. I set up labels names of a line vector and set it up as curved - looks like in this screenshot

But when I want to print it in print composer it doesn't show up

It only shows labels if they aren't curved.

What can I do to fix it. Curved labels are looking really nice (for compared to parallel one.


texts: Clipping a polygon based on feature from overlapping polygon QGIS

texts:
I am working in QGIS
I have two overlapping polygon layers:

A species distribution polygon named species_dist
A bathymetry polygon named depth

The bathymetry polygon attribute table has a depth field with values from 1 to 19.
I want to clip species_dist to only keep the area that overlaps with depth values 1 to 3 of depth.
What is the easiest way to do this?


texts: How to get QgsFeatureRequest in a virtual Layer?

texts:
I am wondering if/how PyQGIS can intercept the QgsVirtualLayerProvider::getFeatures() requests of a VirtualLayer, and pull out just the spatial extent from the QgsFeatureRequest to dynamically formulate a subselect for one piece of a query, to be used as an 'IN' clause of a larger, more complex query ...

I effectively want to select polygon ids from one spatial layer, and use their ids on the fly to select features from another, much larger time series table. Simply making a view joining the two by id doesn't work: for every polygon id, there may be hundreds of rows.

So I want to explicitly pick apart the QgsFeatureRequest into spatial vs attribute components, and execute well-crafted queries appropriately.

I want to be issuing queries to PostGIS that look something like:

select b.price, b.block_id, g.the_geom from hc_models.block_ts b, join tiger_latest.tabblock g on b.block_id = g.tabblock_id
where b.month = ‘2005-01-01’ and b.block_id in

(select tabblock_id from tiger_latest.tabblock WHERE "the_geom" && st_makeenvelope(-122.45845691920965237,37.75162120370001873,-122.44739741920977849,37.75793370369994761,4269)

the subselect against tiger_latest.tabblock will yield a few dozen rows, and the subsequent query against hc_models.block_ts will hit a combined month+id index.

Can this be done in a Virtual Layer? I see a tantalizing comment in the "Index Support" description for Virtual Layers:

  "For spatial indexes, a 'hidden' field named 'search_frame' is created
  for each virtual table (i.e. each referenced layer of the virtual
  layer). The bounding box of the given geometry will be used to
  restrain the query to a particular region of space. This bounding box
  is passed to the underlyng provider by a QgsFeatureRequest with a
  "filterRect".

How to use this?  Or should I just bite the bullet and start with a raw QgsPluginLayer?


texts: DGGS projections with PostGIS or pgLatLon

texts:
Important DGGS projections (like ISEA4T, ISEA4D or ISEA7H variants of ISEA) are free at https://github.com/sahrk/DGGRID , the DGGRID Software Project,
but I not see how to use it with PostGIS. Only maybe drawing and exporting each hierarchical level of the (big data) grid.
Ideal is a function pair like ST_GeoHash / ST_GeomFromGeoHash that encode/decode cell identifiers of the grid.
Some clue or workaround to "plug" DGGRID into PostgreSQL? It can be a plug  for PostGIS or pgLatLon extensions.

NOTES FOR CURIOUS READERS:

Discrete Global Grid Systems (DGGS) are the trend in GIS, and since 2017 has an standard  OGC Abstract Specification.

pgLatLon extension is a lightweight alternative for PostGIS, see  https://www.public-software-group.org/pgLatLon

DGGRID (and ISEA4T, ISEA4D, ISEA7H grids) documentation is at git/sahrk/DGGRID/dggridManualV71.pdf. About first levels of the grid, they are analog to ISEA3H, described here with general conventions for orientation here.  See also DGGRID use case examples.


texts: naturalearth sqlite db not loading into qgis 2.2

texts: I downloaded the naturalearthdatabase in sqlite-format (http://www.naturalearthdata.com/downloads/, the downloadlink goes to http://kelso.it/x/nesqlite), but qgis 2.2 (running under win 7 64bit) does not load it ("unexpected error"), but I can open the db in sqlitequi. Any ideas whats wrong?


texts: Foreign Key is not added to destination table after the creation of a relationship class

texts:
Let's say I have two tables A and B. The relationship between these two tables is as follows :

A (1,*) <---- (1,1) B

Now when I try to add a relationship class for these two tables on Arcatalog, I right-click on the .sde file in Arcatalog and then :

New ==> relationship class.

For the origin table I choose A, and for the destination table I choose B

(Is that right btw or should it be the opposite?)

For the type of the relation ship I select : Simple (Peer to Peer) relationship

Direction : Both

The cardinality of the relationship class : (1 - m) One to Many

And after that I get to select the primary key field in the origin table and the foreign key in the destination table:

As you can see for the foreign key it actually gives me the objectid of that table, now I need to understand if I have to manually add the foreign key field in that table and then select it here?

Shouldn't it be created automatically after the relationship class is created? Obviously when I click on Next and I go back to my postgresql table no foreign key is added.. and it's confusing 'cause I thought that once we define a relationship automatically the id of the origin table migrates to the destination table and becomes a foreign key!

I had to use this database in a web application using ASP.NET, with Entity Framework, and I got confused after Entity Framework has generated the model from my database (at first I added the FK manually).. So now I'm doing it all over again and I would like to do it the right way..

How can I define the foreign keys and the relationship classes from arcatalog?


texts:
My office just scanned a house using the Matterport standard camera. The results were given in OBJ format.

Does anyone know if ArcGIS Pro can read and show OBJ files in 3D?


texts: QGIS 3.14: Wrong GRASS path "boycots" the activation of OTB

texts:
I want to run OTB with QGIS on a MacOS (I work with QGIS3.14) and tried to activate OTB under Preferences > Processing > Providers. I activate OTB and set the OTB-path, but when I try to confirm I got a warning that the GRASS-Folder is not valid - see image below.
The GRASS-Folder (/Applications/QGIS 3.app/Contents/MacOS/grass7) is indeed wrong, but I didn't changed that and have no clue where GRASS is located within QGIS.
I also work a lot with GRASS-Tools and they work well, so no idea how important this path is.
The big problem is, that the OTB settings are not applied with that wrong GRASS parameter.
Any ideas how to solve that?


texts: How to make working with Postgis vector layers in QGIS faster?

texts:
I have just imported a Spatialite database into Postgresql/Postgis. The problem is: using QGIS to work with them is much slower in Postgis (spatial indexes created) than doing the same in Spatialite.

Is there any way to improve Postgis layers speed in QGIS?

EDIT:

Postgis is running on an external server while Spatialite database is on my local PC. Previously, I thought QGIS would cache all my loaded layers onto local disk, I just checked and found that I was wrong - every time I pan around in QGIS, network traffic on my Postgis database increase (right now I am the only one using this experimenting server).

Is there any way to make QGIS cache data locally to improve speed?


texts: Land cover modifications in different time periods

texts:
I have several shapefiles which contain buildings built in different periods from 1948 to 2013. My question is how could I see which buildings for example in were built 1948 but are missing in later periods? Is there a tool in QGIS, or something for automatization? I would think clicking and unclicking layers wouldn't be the best and most efficient option.

EDIT:

Data: shapefiles, aerial imagery (1948, 1971, 1981, 2009, 2013)

I use one shapefile that has a field with the dates of when the buildings were built, such as: before 1948, 1970, between 1971-1981, 1981-2009, built after 2009. They all share the same number of fields and data, the things which differs are the attributes of course. But I'm not sure about this option.

Initially, I was thinking of separating the into 5 shapefiles so as to figure out an easier way to see the differences. As an example: shp 1948 contains (buildings built until 1948), 1971 (all until 1971 including data of previous shp), so on and so forth.

The thing is I used aerial imagery from those years to see whether those buildings still exist or not, so as to update each shapefile correspondingly. Partially done while the data was under a single shapefile, as some were missing.


texts:
Really, I think this a somewhat easy question for someone with some more experience in javascript programming.

I have a Django app to map some locations using leaflet. I'm using dajax for some on mouseover event already I want to integrate some filters using the same library. Problem is, I can't manage to update the map with the filtered data. The template tag {{ geo_json|safe }}; holds the data the map is initialised with. To clarify, here is my js code from my template

var map = L.map('leaflet_map').setView({{zoom|safe}}, 13);
L.tileLayer('http://{s}.tile.cloudmade.com/439299294ed945808af82774cd1b55a0/997/256/{z}/{x}/{y}.png',
{
   attribution: 'Map data &copy; OpenStreetMap, Imagery &copy; CloudMade',
   maxZoom: 18
   }).addTo(map);

var geojsonFeature = {{ geo_json|safe }};

function onEachFeature(feature, layer) {
      layer.on('click', function (e) {

          var id = feature.properties.id
          win_str = '/vsoilmap/popup?id=' + id
          window.open(win_str,'mywindow')
      });

      layer.on('mouseover', function (e) {

          var id = feature.properties.id
          win_str = '/vsoilmap/popup?id=' + id
          Dajaxice.VsoilMap.show_vsoil(Dajax.process,{'pk':id})
          });

      layer.on('mouseout', function (e) {
          Dajaxice.VsoilMap.hide_vsoil(Dajax.process)
      });

}

var geojsonMarkerOptions = {
    radius: 8,
            fillColor: "green",
    color: "darkgreen",
    weight: 0.5,
    opacity: 0.8,
    fillOpacity: 0.8
};

var gj = L.geoJson(geojsonFeature, {
    onEachFeature: onEachFeature,
            pointToLayer: function (feature, latlng) {
               return L.circleMarker(latlng, geojsonMarkerOptions);
    }}).addTo(map);

What I'm basically trying to do is to update the variable geojsonFeature. If I use the dajax add_data method I get the error Map container is already initialized. Of course the whole point of using ajax is to NOT load the map again so my question really is: how can I change the content of a variable in js without running the whole function again?


texts: Recalculating full extent of File Geodatabase feature class?

texts:
I have a featureclass in an Esri file geodatabase. When I edit the features using the editor in ArcMap, and I delete most of my features and keeping only the ones in the middle, the zoom to layer command doesn't work as expected. Instead of zooming to the extent of the features remaining in the feature class, the extent will be the former one (containing all deleted features). When reviewing the extent values in the properties dialog for the feature class, I can clearly see the old values. So the edit session doesn't seem to alter the full extent values.

Is there a way to recalculate these values?

I am pretty sure that this problem should occur to everyone editing features in ArcMap...


texts: Do you know a fairly up-to-date ArcExplorer tutorial video for non-gis staff?

texts: Does anyone know of a fairly up-to-date ArcExplorer tutorial video to show to our non GIS staff who will be using ArcExplorer?


texts: What is missing in my spatial lag equation model?

texts:
I worked it in GeoDa but I think something is missing.

What is the right equation for a spatial lag model?

In my understanding to predict a variable Y is as follows:

Y = w_t_homi100 + CONSTANT + p_Unemploy + p_RenterOc

But, when I use it, I don´t get predicted values as in column PREDICTED. For example, my equation would be:

Y = 0.541574 -14.9127 + (783.381*p_Unemploy) + (25.614*p_RenterOc)

With this, I don't get PREDICTED values (139.59656, 67.47219, 90.87614) as in the picture.


texts: Creating a new raster based on pixel value of original raster using raster calculator but new layer is not evident

texts:
I am trying to create a new raster based on a current raster (w001001 - NVIS) which has one paletted band, with palette values ranging from 1 to 29, with each palette/value representing a vegetation type (e.g. 1 = rain forests).

I want to create a new raster (actually, preferably a vector, but I cross that bridge later) for the palette of value 3 (this one):

I have followed the answers to the questions posted in this forum, and have used the Raster Calculator to do this:

It all seems to process fine, and the new layer appears in my layer list. However, I can't see the layer at all and when I use the  nothing comes up. When I open properties for the layer, metadata tells me this:

Is there something obvious that I am missing?

This is my first post and I am a fairly new user of QGIS so apologies if anything is unclear or nonsensical.


texts: Why is a (self) join so slow with Virtual Layer

texts:
I use QGIS 3.4 on Windows 10
I have a polygon layer and a point layer. The points have a streetname attribute. I want to check if there are polygons with more then 1 streetname.
I took 2 steps, first Join attributes by location to get the streetnames as attribute in the polygons.
The second step is a self join in a virtual layer
SELECT T1.POLid, T1.Postalcode, T1.streetname, T2.streetname
FROM TestStreet AS T1, TestStreet AS T2
WHERE T1.POLid = T2.POLid AND T1.streetname <> T2.streetname;

I have to interrupt the processing after waiting very long.
My workaround is:
Save features As csv file
Import csv file in MSAccess
MSAccess gives instantaneous results with the same query.

What is it what I do wrong in QGIS, is there a better way just with QGIS?


texts:
Is there anyway to retain/easily transfer symbology from ESRI to Mapinfo? The only process I found was posted in 2010 - https://groups.google.com/group/mapinfo-l/browse_thread/thread/883223ff4ae56809 and it is rather convoluted and supposed to be unreliable.

ArcGIS 10 and Mapinfo 11 are the core programs.

This is related to Converting ArcGIS style (symbology) to font for MapInfo/OpenGeo? and Creating summary table of symbology of all datasets in multiple MXDs? as possible workarounds.


texts: Style graduated: Clases yes, Illustartion no

texts:
I joined "1980" (Fields on foto) with ¨municipal¨ (the blue one graduate) through ¨CVEGEO¨ value. Then I tried to illustrate with graduate style the  AUTO_PARTICULAR value.

Result: There is a clasification but not the illustration.

Before the aplication of graduate style, I had an another graduate style illustration of a value but of the same table of ¨municipal¨, that was good.

The problem is with the joined table. I recently imported through ¨Add Delimited Text Layer¨ with a .csvt file.


texts: Calculate Great Circle Distance between two point layers via way point

texts:
I need to create three distance matrices for around 550 points spread throughout Africa, Europe and Asia. The first matrix represents great circle distances, the second is the shortest distance along a network of trade routes, and the third is the shortest great circle distance via a way point.

I've completed the first two and was wondering if anyone knew of a method in ArcGIS 10.3 to calculate the distances between every point in  two point layers via a single point. For instance, I have a group of points in southern Africa (Layer 1) and another group of points in Kazakhstan (Layer 2). I have a waypoint in Cairo, Egypt (Layer 3). I need to calculate the great circle distance from each point in Layer 1 to the waypoint in Layer 3 and then the great circle distance to each point in Layer 2.

Right now I'm planning on using the Generate Near Table tool twice, once between Layer 1 and the waypoint in Layer 3, and then again with Layer 2 and the waypoint in Layer 3. I'll then combine the distances accordingly in excel. This will take a very long time since the distance matrix has roughly 250,000 lines, so if there is a tool that could calculate great circle distance via waypoints it would save me a ton of time.

Please let me know if you need any more information.


texts: Copying database view into an enterprise geodatabase?

texts:
Say I have DatabaseA, which isn't structured in a way that allows me to reference it as a geodatabase in ArcGIS. Fields have to be renamed, geometry added, etc. etc.

Now, it seems to me that what I should do is create a geodatabase, and then replicate the data from DatabaseA into what will be the new DatabaseB. Of course there is the Create Enterprise Geodatabase tool, but I still need to get DatabaseA information into it.

Would it make sense to create a database view from DatabaseA, have it set to automatically populate/append into DatabaseB, and then have ArcGIS connected to DatabaseB (the created geodatabase), or am I missing some of the workflow here?


texts: Geoprocessing Tools not producing an output in ArcGIS Desktop?

texts:
I've been running a few different geoprocessing tools in ArcGIS Desktop 10.5, which seem to operate as required, the various processing windows appear to give the impression that works is going on but when I get the final tick window saying the process is complete, there is no output in the specified location.

Am I missing something here?


texts: Using relationship classes with ArcGIS Online?

texts:
I am using ArcInfo 10.4 to create geodatabase with a relationship class associating a point feature to a table with the goal of using it to collect inspection data in the table of the point features through ESRI Collector.  (I am a newbie at ArcGIS but have a lot of Access database experience.)

As I have it set up right now, I can create point feature records in Collector, and can create a new inspection table record from the point feature in Collector, but I can't see why the Global ID (parent key of the GUID field type) from the site feature is not being copied into the rel_Global ID field (foreign key also GUID field type), so the inspections table has the records, but they are not associated with the point feature.


texts: Use buffer tool to draw different radii?

texts:
I am new to Arcmap-Arceditor 10 and could really use your help. I have been mapping ethanol plants (which utilize corn to produce ethanol) and corn production amounts in the counties of the USA.

My goal is to draw buffer zones (circles) around the ethanol plants (dots). However, the tricky part is that the circles must reflect what region an ethanol plant buys corn from in order to sustain its ethanol production. I know how much corn each plant needs and I know the average corn production/square mile in each county. How do I get Arcmap to calculate the radiuses of the buffer zones it needs to draw corn from? This gets especially tricky when I have an ethanol plant lying on the border of two counties and each county has a different corn production amount/square mile.

Do you know if it's possible to do this in arcmap? If so, how should I proceed? Thank you so much!

In the image attached, the dots are ethanol plants. The red regions show the county is self-sufficient in producing enough corn for its ethanol plants. The orange regions show that corn must be imported from other counties to support its ethanol plants.


texts: I am new to QGIS. I have set up several desktops with QGIS to be used by mostly GIS novices in our agency. One of the problems I have encountered is that they accidentally close off one or another of the panels (Legends TOC mostly) and then can't get it back. I would like to set up a configuration that has these panels in tabs how I want them. My question is can I lock this configuration so they cannot inadvertently close off the panels?


texts:
I have an ArcMap 10.3 MXD in which which I am trying to update the definition query of an enterprise geodatabase layer using Python:

import arcpy, os, re
from arcpy import env

inmxd = r'C:\Temp\Tenure_Shire Notification.mxd'
tenID = 'E40/01300'
env.workspace = os.path.dirname(inmxd)
mxd = arcpy.mapping.MapDocument(inmxd)
df = arcpy.mapping.ListDataFrames(mxd, "Main Map")[0]
tenLMAT = arcpy.mapping.ListLayers(mxd, "Tenement Application - Live Mining Act Tenure")[0]
oldExp = tenLMAT.definitionQuery
oldTenID = re.split(r"\'", oldExp)[1]
newExp = oldExp.replace(oldTenID,tenID)
tenLMAT.definitionQuery = newExp

My issue is that using this code, my definition query has not changed when I check the layer properties. What am I doing wrong?


texts:
I have ~1.7 GB JP2 Buckeye that I want to serve out from Geoserver; I'd like to have it tiled so that the client won't have to redraw everytime they zoom into the map. My current understanding is that the best way to serve this is by building an Image Pyramid out of this so Geoserver can dish it out--however I think the only way to build the pyramid is by using a GeoTiff.

So I have started converting that <2 GB JP2 file (using gdalwarp) and currently I'm at 30% and the file size >13 GB which is HUGE. Plus if I go ahead and build the pyramid out of that file, I feel like the total size of this one file is going to be massive.

Is there a better way that I should serve this data? Ideally I should be using around the same amount of harddrive space as the JP2 but if I'm missing something, I welcome corrections and explanations so I can understand.

UPDATE
cmd line used:

gdalwarp -t_srs "EPSG:3857" -multi source.jp2 output.tif

UPDATE per the suggestion of @iant to load up the bugger as a layer (attempted to add it as JP2K (Direct) Coverage Format Raster Data Source), i got the following error:

Could not list layers for this store, an error occurred retrieving them: Dimensions (width=111200 height=118200) are too large


texts: How to generate random points with minimum distances between points in multiple polygons

texts:
I am looking to create random points for a random stratified survey. There are three different categories being surveyed; for each, I have a single shapefile, each containing a number of polygons.

I would like to generate a number (undetermined) of random points in each shapefile, with a minimum distance of 250m between points (regardless of what category/shapefile they are in).  In other words, I would like there to be an equal number of random points in shapefile1, shapefile2, and shapefile3, with no point within 250m of another (i.e., a point in shapefile1 should not be within 250m of another point in shapefile1 or a point in shapefiles2 or 3). If possible, I would like to create the maximum number of random points possible given my minimum distance of 250m limitation.

So far, I have been able to generate random points (using ArcMap and qgis), however the minimum distance between points seems to apply only to points within each polygon, rather than all polygons in three different shapefiles.


texts: Which ESRI license is needed for >4gb rdbms for 1 user?

texts:
I want to create our >4GB geodatabases that interface with our business SQL Server databases.  Essentially create enterprise geodatabases.  Can this be done with a Standard desktop license without need for arcGIS Enterprise (Server)?

We would not need a multi-user editing environment.


texts: Obtaining the average red, green, blue values within grid cells in QGIS?

texts: Is there a way to obtain the average red, green, blue values within grid cells in QGIS?


texts: Meaning of zero "0" value in raster layer produced by Soil Texture Classification tool in QGIS?

texts: How to interpret zero "0" value in raster layer produced by Soil Texture Classification tool in QGIS? It cannot be found in classification table.


texts: Google Earth Engine: Error on point value extraction

texts:
I am trying to extract rainfall (CHIRPS) values for a set of locations but I got the following error:

  Image.reduceRegions: Unable to find a crs

This does not happens with other datasets (such as terraclimate)

The error arise when i start downloading the table from Tasks.

Here is the link to the code I have run https://code.earthengine.google.com/c27f2156e81824b0990dcfe0b0a6f455

The error should be here:

// do extraction
var ft = ee.FeatureCollection(ee.List([]));

//Function to extract values from image collection based on point file and export as a table
var fill = function(img, ini) {
  var inift = ee.FeatureCollection(ini);
  var scale = ee.Image(MM.first()).projection().nominalScale().getInfo()
  var ft2 = img.reduceRegions(pts, ee.Reducer.first(),scale);
  var date = img.date().format("YYYYMM");
  var ft3 = ft2.map(function(f){return f.set("date", date)});
return inift.merge(ft3);
};

// Iterates over the ImageCollection
var profile = ee.FeatureCollection(MM.iterate(fill, ft));


texts: Unable to connect ArcGIS Online through ArcMap/ ArcCatalog ver.10.2

texts:
I have ArcGIS Desktop 10.2.I am trying to connect to ArcGIS Online from Arc-catalog (My hosted Services) but throws me an error saying to check Network connection.

  Error : " Error Has Occurred . Check your Network connection"

But I am able to connect to Arcgis online through web browser. I have ArcGIS Global account & ArcGIS online account as well.


texts: Flattening lakes for hydro-flattened DEM using polygon vertices in QGIS

texts:
I have the vertices around lake polygons that I want to flatten for a hydro-flattened DEM. I know that ArcGIS has a trend tool that will create a sloped 'plane of best fit' for a series of points, however the best I can find in QGIS is the IDW interpolation tool, which doesn't create a completely flat surface.

Does anyone know of a QGIS tool/workflow that will do this?

Edit: Lake surfaces should be completely flat, it's large rivers that require a trend surface.


texts: Using Python ArcGIS for Developers without using Jupyter Notebooks

texts:
I want to create an ArcGIS map in Python, add layers to it based on data taken from a .csv file, and then display the map in a webpage.
Unfortunately, every Python tutorial I could find about the ArcGIS for developers needed to use something called a Jupyter Notebook. As far as I can tell, this would not work well with all the preexisting Python files I use in my project in other IDE's. I also do not like the idea of using an online IDE.
I have already successfully gotten the Python ArcGIS library in another IDE (PyCharm), and it is able to print out data and do most things normally. Except I can't find any way to actually visually display the maps or insert them into a webpage.
Am I going about this wrong, and is it impossible to avoid the Jupyter notebook? If not, how can I create an ArcGIS map and display it in a webpage in Python?
I'm very new to GIS.


texts:
I am trying to generate a normalized DSM by subtracting the original DSM (oDSM) from the DTM.

To generate the DTM I had to first use a median filter on the oDSM (remove noise) and then use a morphological filter on the cleaned oDSM. For the morphological filter I had to convert the oDSM to an 8bit from float32. Now, I need to finally subtract the oDSM from the morphologically operated DSM.

I need to know if subtracting two different data type rasters is a problem. If yes, how can I convert the 8 bit back to a GTIFF? I am working with QGIS and R.


texts:
I'm using Landsat 8 imagery in Python and want to make true color images that would look good in a publication. I have many images, so ideally this would be an automated process where each image comes out looking not too dark, not overexposed, and with good contrast.

What I've been using is a linear stretch of 1.25%, but some of the images still look washed out.

from skimage import exposure
import matplotlib.pyplot as plt
import numpy as np

r = np.random.randint(0, 255, size=10000).reshape((100, 100))
g = np.random.randint(0, 255, size=10000).reshape((100, 100))
b = np.random.randint(0, 255, size=10000).reshape((100, 100))
rgb = np.dstack((r, g, b))

def linearStretch(input, percent):
    pLow, pHigh = np.percentile(input[~np.isnan(input)], (percent, 100 - percent))
    img_rescale = exposure.rescale_intensity(input, in_range=(pLow, pHigh))
    return img_rescale

img_rescaled = linearStretch(rgb, 1.25)
plt.figure()
plt.imshow(img_rescaled)

Is there a known image enhancement technique that can be used to reliably process images automatically without much/any individual tweaking?

All of my images have had clouds removed and replaced with NaN.


texts: What is the minimum sample size when sampling spatialy from existing points

texts: I have sampled points from an existing shapefile using the Generate Near Table proximity feature to calculate the distance that I later used to select points. Depending on the distance interval I used, my sample size is 27% of the existing points (if I want to increase the sample size, I can reduce the distance interval). However, I would like to find out is there a minimum sample size for this kind of spatial sampling?


texts: ESRI Web Appbuilder developer edition vs JavaScript API

texts: What is the benefit to using one over the other besides the out of the box widgets and functionality offered by Web AppBuilder?  If you can customize the code in Web AppBuilder, what really is the difference?


texts: How to create a shapefile [polygon type] from a Geodataframe, returned from a Oracle Spatial cursor with geometry column type=cx_Oracle.LOB?

texts:
I'm trying to create a shapefile [2D polygon type] with a Geodataframe, that result from a SQL made to a Oracle Spatial database.
The SQL and the Geodataframe creation runs without errors, but when I try to use the GeoDataframe[.to_file] function, returns a few errors in the command:

geoparc.to_file('.../out/export.shp', driver='ESRI Shapefile')

Am I missing any parameters for the geometry type in .to file function?

How can I set the geometry type in the .to_file function?

geometry data type:

print(type(geoparc.geometry[0]))
 = class 'cx_Oracle.LOB'>

print(geoparc.geometry[0]) =
POLYGON ((225117.58 411950.12, 225057.26803977 411941.365038031, 225058.343749989 411940.968749974, 225057.812500075 411939.468749988))

command:

geoparc.to_file('.../out/export.shp', driver='ESRI Shapefile')

Errors:
 geoparc.to_file('.../out/export.shp', driver='ESRI Shapefile');

...Python\Python35-32\lib\site-packages\geopandas\geodataframe.py", line 413, in to_file to_file(self, filename, driver, schema, **kwargs);

...Python\Python35-32\lib\site-packages\geopandas\io\file.py", line 106, in to_file    schema = infer_schema(df);

...Python\Python35-32\lib\site-packages\geopandas\io\file.py", line 143, in infer_schema    raise ValueError("Geometry column cannot contain mutiple "
ValueError: Geometry column cannot contain mutiple geometry types when writing to file.


texts: Creating simple citygml 3D models based on 2D shapefiles (alkis) shp2gml using Python?

texts: How can I create a simple LoD1 citygml-model using Python?


texts:
I am looking for an opensource web based metadata management tool for geographic data. I have gone through a related question: Seeking Open Source ISO Metadata Editor?
First of all it was asked six years ago. And it recommends tool like Geonetwork. Geonetwork at first looked promising but has several issues while installing (tried on windows). i am unable to install it despite asking for help on several forums.

I have gone through this link listing quite a few tools : https://wiki.osgeo.org/wiki/Metadata_software

If anyone has had experience with any of it kindly share.


texts: Creating an Arcade expression within a feature service

texts:
I am trying to query a set of points in relation to a polygon feature service.
For example, I have dataset "A" within an RDMS that I want to query and get a count of features within Dataset "B" which is a polygon feature service.
I have tried a similar expression and two datasets from the same RDBMS which worked, but not with a foreign dataset.
Any suggestions?


texts: How to overlay points from MySQL served by GeoServer on an OpenLayers map?

texts:
This might be a very basic question as I am new to GeoServer.

What I am trying to do:

I have a MySQL table that stores data attached to a spatial point and have latitude and longitude columns to store the location attribute.
Use GeoServer to display these data points as layer.

What has been done:

Create a new table which is duplicate of the original table but also have POINT field which is generated using POINT(orig_lat, orig_lon).
Use the MySQL extension to create a new data store linked with this spatial_table.
Create a new layer here we need two fields called as native and declared SRS, in the demo I had seen EPSG:4326 as the declared SRS and hence I have added that SRS.

Problem:

When I view the layer preview in OpenLayers I do not see any map? Not sure if this is because of SRS or do I need to specify a map layer.
Most of the points seem like they are overlaid on top of each other and zooming and panning does not make much difference. This seems like SRS issue.
Is it possible to overlay the points on top of Google Maps using GeoServer?


texts:
Could you please help with selection problem. I have about 1000 of shapefiles. Each shapefile has points located in a rectangular net. But size of that nets (and as a result number of points) differ in every shapefile.
Only thing all shapefiles have in common are distances, net size is 10 meters.
Is there any way to select points in chessboard like pattern? Image below there is an example, I need to select all rows (not only first three as in image) in that pattern.

edit:
Second picture with labels and result of selecting odds by: MOD ("ID",  2)  = 1
 suggested in answer


texts:
I have a script that does some processing but the main part are two things

run a multiple ring buffer using arcpy.MultipleRingBuffer_analysis and four rings
do some zonal statistics using arcpy.sa.ZonalStatisticsAsTable on the four rings

I do this repetitively (through a loop) so I want to dump the results to a table called Final Results.  It makes sense to me to have the results organized in five columns (PatientID, Distance1, Distance2, Distance3, Distance4).  I making my loop off a search cursor and keeping the PatientID.

I thought to create a pivot table using arcpy.PivotTable_management where the result becomes one row and five columns (the pivot variable and four buffers) and then push that to Final Results.  However, it appears I would still need to create a search cursor for my pivot table, grab the first (and only) row, then create an insert cursor to put the data into my Final Resultstable.

Alternatively, I suppose I could just go to the table from #2 above, create a search cursor, and assign four variables (D1, D2, D3, D4) to the zonal stats results with some sort of switch-case statement, then create an insert cursor that concatenates D1, D2, D3, D4 and uses insert row.

Both seem rather inefficient.

Is there a smarter way to just push my results from zonal analysis into a table that I use?


texts: How is slope of a surface in an arbitrary direction produced?

texts:
I am currently developing a program that requires that I calculate the slope for each cell in a DEM. However when I try to understand the algorithm to be used I get stuck trying to understand exactly what is going on. The slope algorithm in question is the one detailed in the article Hillshading and the reflectance map (Horn 1981 p. 20-21).

The algorithm is also explained here

I understand how the average height factors for the x and y direction gets calculated and I also understand that the factors have to be combined if the slope of the surface in an arbitrary direction is to be established. I dread that the part which keeps troubling me is probably the easiest.

The part I do not understand is why the Pythagorean theorem is used When the two factors are combined into a single average value? Both dx and dy will be values rising towards z, so how can the Pythagorean theorem be applied to two lines along the same axis? I've tried to visualize the dx and dy values as three-dimensional vectors with dx being (1, 0, dx) and dy being (0, 1, dy). Adding these two vectors together produces the two sides in a right triangle where the hypotenuse will be equal to the answer the equation is supposed to give me. I am then left with the hypotenuse vector as the rise and that just feels wrong when I attempt to perceive it as the rise. Am I misunderstanding something here?


texts: R: Remove a Raster Attribute Table from raster for writing to Tiff

texts:
I've been given a large raster with an in-built RAT with all sorts of RGB data and classes etc etc, but it's been done a bit poorly and I would like to quickly rebuild it myself.

When i read into R and then writeRaster to GeoTiff, it's a right old mess, so i'd like to remove the RAT entirely, how is this done?

# here's a raster with RAT
require(raster)
r <- raster(ncol=10,nrow=10,xmn=0,xmx=10,ymn=0,ymx=10)
r[] <- sample(1:4,100,rep=T)

# ratify and add some attributes
r <- ratify(r)
rat <- levels(r)[[1]]
rat$letters <- LETTERS[1:4]
rat$names <- c("Tom","Tim","Jane","Jessica")
rat$vals <- c(0.265,0.456,0.037,0.756)

# set RAT
levels(r) <- rat

How do i get rid of the RAT?


texts: How to animate a sequence of maps in postscript format?

texts:
I'm working on the Atlantic ocean for my research and I need to create an animation of movement of few islands in the Atlantic from their past locations to their present locations that moved due to plate tectonics. Input will be a postscript file of that region.

Are there any tools which can do this task?


texts:
My development environment is Visual studio 2010 with ArcGIS Desktop 10.1 installed.
I need my add-in to support 10.0 AND 10.1. Right now, if I try to install it on ArcMap 10.0, it fails.

Is it possible to target an ArcMap add-in to an earlier version?


texts: Select feature by id in QGIS using Python

texts:
I have a ComboBox shows feature ids of selected features in a layer. I am trying to select one feature using ids in ComboBox currentIndex. I used setSelectedFeature(), but it didn't worked.I can't figure out how to get it.

Here is my code:

layer = iface.activeaLayer()
feats = layer.selectedFeatures()
for feat in feats:
    name = feat["id"]
    self.dlg.cb_Feature.addItem(name, feat.id())
    onefeat = self.dlg.cb_Feature.itemData(self.dlg.cb_Feature.currentIndex())
    layer.setSelectedFeatures(onefeat)


texts: QGIS Contour Extraction Malformed Path (QGIS 3.0.1 Girona)

texts:
I'm attempting to create contours from a raster DEM (provided by USGS).  I open the dialog at Raster > Extraction > Contour, choose the input level, leave the default Band number of Band 1, set an interval of 10 feet, attribute name of ELEV, offset from zero of 0.0000, and select an output directory and file using the dialog from the ellipsis next to the "Contours" field.  However, I get an error message

ERROR 1: Failed to create directory "C:/Users/Username/Documents/QGIS/Cougar/Vector/w001001.shp" for shapefile datastore.

I noted that the GDAL command is:

gdal_contour -b 1 -a ELEV -i 10.0 -f "ESRI Shapefile" "C:\\Users\\Username\\Documents\\QGIS\\Cougar\\raster\\1-3degree2014elev\\grdn45w123_13\\FEET-w001001.tif" "\"C:/Users/Username/Documents/QGIS/Cougar/Vector/w001001.shp\""
(This has been lightly edited to avoid exposing my username; otherwise, an exact cut-and-paste)

This appears to me to have a malformed path "\"C:/Users/Username/Documents/QGIS/Cougar/Vector/w001001.shp\"", and when I copy this into a command line, and edit the command to read:

gdal_contour -b 1 -a ELEV -i 10.0 -f "ESRI Shapefile" "C:\\Users\\Username\\Documents\\QGIS\\Cougar\\raster\\1-3degree2014elev\\grdn45w123_13\\FEET-w001001.tif" "C:/Users/Username/Documents/QGIS/Cougar/Vector/w001001.shp"

the command runs as expected and produces output which I can successfully add to my project in in the QGIS GUI.  What am I missing in the Contour dialog to produce proper output?  Obviously, I can work around whatever I'm doing wrong, but I'd prefer to not have to mess with a bunch of cutting and pasting and editing on the command line if I can help it.  The same thing occurs when trying to use the "Run as a Batch Process..." dialog.


texts: Industry employment and GIS: effective search methods for new grads?

texts: What are some good ways of identifying private entities that may be inclined to accept intern work or offer an entry level job?


texts: Extract attachments with increment for duplicate values

texts:
I am using Arcgis Desktop 10.4.1 and Arcgis Collector.

I have joined the globalid and rel_globalid as explained here

I almost what I want, except for the counter for the duplicated .jpg files. I want them to have a suffix like _1,_2,_3. This is my result:

How can I add a counter that does the naming in the way I want, and still extracts the attachments?


texts:
I am trying to get data from a PostGIS datatable using .NET, pass it to javascript, create a layer in OpenLayers and add the PostGIS data to the layer. I can successfully query from the PostGres/PostGIS table using the following SQL;

SELECT ""tract_name"", ST_AsGeoJSON(""geom"") FROM the_layer;

After querying the database the result ends up in the PostGres DataAdapter as a DataSet.

I was hoping to simply pass the data to javascript where I could use the following to add the data;

        var geojson_format = new OpenLayers.Format.GeoJSON();
        var vector_layer = new OpenLayers.Layer.Vector("query_POSTGIS_layer");
        map.addLayer(vector_layer);
        vector_layer.addFeatures(geojson_format.read(postGISqueryResult));

What I can't figure out is how to convert the PostGres DataSet into a GeoJSON format that can be handled by the code above.

Application stack is PostGres/PostGIS, GeoServer, and OpenLayers.