texts: Apply a dynamic style to a raster(GeoServer, OpenLayers)

texts:
I want to apply a dynamic style to a raster using either GeoServer or OpenLayers.
I have found similar examples at OpenLayers plugins but these examples work on RGB images instead of the real values of the raster.
An example of a similar tool that I want to develop:


texts:
I'm looking for a tool which compares two vector files (eg shapefiles). Eg 2 versions of a river network.
Quite a lot of the records are equal, but some are different. A possible output would be eg a shapefile with every geometry and 2 columns which refer to the id of the 2 original files, which can be missing if that geometry is not present in one of the files.

I'm thinking about a strategy to program this myself (basically creating a hashed version of the geometry of each part and compare those), but before rushing into programming (I even have a fancy name 'shapediff') I was wondering if such a tool already exists.


texts: I have an excel file with x and y. how can I use this excel file in qgis' python console as point vector layer to make polygon with all the points in the file.(it is important for me to use python)


texts: How do I order polygons so their IDs are contiguous from NW to SE?

texts: Census Tract IDs are ordered in large areas, but patchy across the state. I need to assign a new ID that flows in an order from one end to the opposite end (i.e, NW to SE) where polygons are touching. Any suggestions?


texts: Filtering a layer based on property that's an array in with Mapbox GL JS?

texts:
I've got a relational dataset broken into two layers: members and organizations. They're GeoJSON, with features that look something like this:

// Org Features
{
  "type": "Feature",
  "geometry": {
    "type": "Point",
    "coordinates": [-78.34998, 37.784]
  },
  "properties": {
    "id": 46391,
    "name": "A Big Organization",
    "lifetime_profiles": 7013,
    "present_profiles": 49
  }
}
// Member Features
{
  "type": "Feature",
  "geometry": {
    "type": "Point",
    "coordinates": [-80.99998, 47.7624]
  },
  "properties": {
    "id": 21446,
    "given_name": "John",
    "family_name": "Doe",
    "orgs": [21437, 117028]
  }
}

I want to show which members are linked to which organizations when I click on the map. I'm trying to do so using the following filter:

const handleClick = (target) => {
  const bbox = [
     [target.point.x - 5, target.point.y - 5],
     [target.point.x + 5, target.point.y + 5]
  ];
  const opt = { layers: ['orgs'] };
  const targetIDs = map.queryRenderedFeatures(bbox, opt).map((feature) => {
    return feature.properties.id;
  });
  // find clicked orgs
  map.setFilter('users', ['in', 'orgs', ...targetIDs]);
};

From my earlier research, I've gathered that the "in" filter is deprecated now. How can I accomplish this aim using current filters?


texts: Calculating area of intersection of polygon to another layer in ArcGIS for Desktop?

texts:
I am very new to GIS, especially ArcGIS.

I am currently using ArcGIS 10.2, and I have 2 layer, namely the purple feature map(this layer is an import from arcgis online, it doesn't seems to be an feature map) and voronoi diagram layer.

I am currently doing a project which requires me to find the intersection area of the purple map and each slices of voronoi polygon (for example, the highlighted polygon in the sample image) respectively.

I will expect answer like(if it is possible):
Polygon 0: 120km^2
Polygon 1: 80km^2
Polygon 2: 20km^2
.....

Since there is only 9 polygon slice, I don't expect it to require python scripting. I hope someone can point me to the correct toolbox/feature to achieve the objective required.


texts:
I'm trying to analyze set of samples into vector grids in order to generate average data in corresponding grid.I'm working with 50m separated grids but the size of grid may change over time.That's why we need to generate grids on the fly.There is a snapshot of data distribution below.

There are two important things for us,one of which creating grid,as you may know we can create grid from minimum bounding rectangle (MBR) of our data.The trick is let say we have 1000 rows of distributed samples (points) for a large area,then the processing of this MBR plus dividing into grids will be too long.Since,what we want is intersecting grids alone,we must generate them from samples which means we must count as the length of our sample data.That's what I'm asking,is there an easy way to do it without processing all grids?
The other is bound to first where our samples contain various data that we should aggregate into the single grid.How can we do it with or without processing grid at the same time ?
NOTE: This question may require programming or querying which we are wiling to migrate to solution in any ways.Our samples are in text format that we can transform them to any other solutions you can suggest.


texts:
I'm trying to add a WMS layer in EPSG:2180.

I get the definition for EPSG:2180 from here and I am using this guide on GitHub to
make a new CRS definition, but I think it is corrupted. It is not working properly.

I made an example at jsfiddle, so if someone can have a look at it I will be grateful.

I want to know how to display wms layers with leaflet in EPSG:2180

Updated example with only one map with corruped CRS definition.

http://jsfiddle.net/udsE2/4/

Old example with 3 maps, on second map is same layer like in third but in EPSG:4326

http://jsfiddle.net/udsE2/3/


texts:
Is there a way to query basemap to extract all coastal coordinates? Say user provides lat/lng and the function returns true/false if the coordinates are within 1km from the coast?

Links are more than welcome.


texts: How to list all of feature class in SDE?

texts:
I am getting a problem: my SDE contains some feature classes, but after using the statement:

var namesEnum = workspace.DatasetNames[esriDatasetType.esriDTFeatureClass];
var tempName = namesEnum.Next(); // get value null

The value of tempName is Null, but it is ok if I use esriDatasetType.esriDTFeatureDataset.


texts: ListBookmarks does not output any features

texts:
I'm trying to use ListBookmarks in Python to output my bookmarks to features. I directly copied the code in the sample provided. It "succeeds" and creates the feature class I wanted, but it is empty with no features. Did anybody have the same problem and solve it? Using ArcMap 10.4, Python 2.7.10, and Windows 8.1.

Here's the code I used, but it's essentially the same as the sample:

import arcpy, os

# The map with the bookmarks
mxd = arcpy.mapping.MapDocument(r"D:\GIS\Projects\ACE\Easement Maps\UtilityOutgrants2.mxd")

# The output feature class to be created -
# This feature class will store the bookmarks as features
outFC = r'D:\Default.gdb\Bookmarks'

# A template feature class that contains the attribute schema
# Including a "Name" field to store the bookmark name
template = r'D:\Default.gdb\Template'

if arcpy.Exists(outFC):
    arcpy.Delete_management(outFC)
arcpy.CreateFeatureclass_management(os.path.dirname(outFC),
                                    os.path.basename(outFC),
                                    "POLYGON", template,
                                    spatial_reference=template)

cur = arcpy.da.InsertCursor(outFC, ["SHAPE@", "Name"])
array = arcpy.Array()
for bkmk in arcpy.mapping.ListBookmarks(mxd, "Layers"):
    array.add(arcpy.Point(bkmk.extent.XMin, bkmk.extent.YMin))
    array.add(arcpy.Point(bkmk.extent.XMin, bkmk.extent.YMax))
    array.add(arcpy.Point(bkmk.extent.XMax, bkmk.extent.YMax))
    array.add(arcpy.Point(bkmk.extent.XMax, bkmk.extent.YMin))
    # To close the polygon, add the first point again
    array.add(arcpy.Point(bkmk.extent.XMin, bkmk.extent.YMin))
    cur.insertRow([arcpy.Polygon(array), bkmk.name])
    array.removeAll()

print "done"


texts: Assigning fishnet attributes by row instead of by column with ArcPy?

texts:
I created a fishnet grid and assigned numeric attributes 1-16 to a new field using the arcpy.da.UpdateCursor function. It however assigned the attributes north to south (see image).

How can I change it so that the attributes are assigned from top-left to right i.e row by row as opposed to column by column? I know this could be achieved by the grid index tool but I used the fishnet tool because it allows me to draw the grid starting at the northwest corner drawing to the east which the grid index tool cannot achieve.

with arcpy.da.UpdateCursor("myfish1",("OID","switch")) as cursor:
for row in cursor:
    row[1]=row[0]
    cursor.updateRow(row)


texts: RuntimeError: Object: Error in executing tool when running Union

texts:
I have a list of a 2 layers (nfhl_fld_haz_features, prelim_fld_haz_features) and a feature class (master_merge) that I want to perform a union on:

 wkflw2_union = arcpy.Union_analysis([master_merge, nfhl_fld_haz_features, prelim_fld_haz_features], wkflw2_union_output)

I keep getting the following error:

  Traceback (most recent call last):

  File "", line 1, in
      wkflw2_union = arcpy.Union_analysis([master_merge, nfhl_fld_haz_features, prelim_fld_haz_features], wkflw2_union_output)

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\analysis.py", line 734, in
  Union
      raise e

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\analysis.py", line 731, in
  Union
      retval = convertArcObjectToPythonObject(gp.Union_analysis(*gp_fixargs((in_features,
  out_feature_class, join_attributes, cluster_tolerance, gaps), True)))

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing_base.py", line
  496, in
      return lambda *args: val(*gp_fixargs(args, True))

  RuntimeError: Object:  Error in executing tool

I don't see any issue with using a list, but could this error be due to using a feature class mixed with 2 layers? Or could this be something else (the main files come from a network drive)


texts: Can .lyr files be used without shapefiles?

texts:
I've been sent some .lyr files, with no associated shapefiles.
I'm using QGIS, so it seems I need to try some conversions to get these to work anyway - but I just wanted to check before I embark on that journey, if .lyr files are supposed to be used independently of shapefiles?
Or if it would be a normal request to ask for the associated .shps!


texts: How to open most recent directory using PyQGIS?

texts: For example, when I add a vector layer, next time that I want to add another layer it opens up where I chose my last layer from. How would one implement this in Python and/or PyQgis? I'm writing a plugin, and I want user to select layers in the same way as above.


texts:
I am trying to insert a value for different group of percentile eg (0- 10 percentitle) value = 1.

I saw an example over here and went to try it out How do I Calculate Grouped Percentiles or Batch percentiles in ArcMap?

However, i was thrown this error. I am using ArcGis 10.3.1 for desktop. and the file is Shapefile Feature Class.

Runtime error
Traceback (most recent call last):
File "<string>", line 30, in <module>
File "<string>", line 10, in CalcPercentile
File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\lib\function_base.py", line 3096, in percentile
return _compute_qth_percentile(sorted, q, axis, out)
File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\lib\function_base.py", line 3132, in _compute_qth_percentile
return add.reduce(sorted[indexer]*weights, axis=axis, out=out)/sumval
TypeError: unsupported operand type(s) for *: 'numpy.ndarray' and 'numpy.ndarray'

The source code is over here. (using python)

 import arcpy
 import numpy as np
 import os

  #loop through all Shapefile in a folder and call the CalcPercentile method
 def CalcPercentile(inputFeatureClass):
  #to create 3 rank for example
    print inputFeatureClass;
    arr = arcpy.da.FeatureClassToNumPyArray(inputFeatureClass, 'NEAR_DIST_')
    p1 = np.percentile(arr, 33)  # rank = 0
    p2 = np.percentile(arr, 67)  # rank = 1
    p3 = np.percentile(arr, 100)  # rank = 2
    #use cursor to update the new rank field
    with arcpy.da.UpdateCursor(inputFeatureClass , ['NEAR_DIST_','PerRank1']) as cursor:
    for row in cursor:
         if row[0] < p1:
             row[1] = 0  #rank 0
         elif p1 <= row[0] and row[0] < p2:
              row[1] = 1
         else:
              row[1] = 2
         cursor.updateRow(row)

workspace = "D:\HDB_accessibility\MergeAllDistance"
walk = arcpy.da.Walk(workspace, datatype="FeatureClass")
for dirpath, dirnames, filenames in walk:
 for filename in filenames:
     featureClass = os.path.join(dirpath, filename)
     #First add the percentile Rank Field
     #arcpy.AddField_management(featureClass, "PerRank", "DOUBLE")
     CalcPercentile(featureClass)


texts: Attribute table trailing zeroes missing in Print Composer

texts:
I have QGIS 2.18.3 installed. I'm trying to create a printable .pdf map with QGIS Print Composer. It needs to have X & Y coordinates listed in an attribute table, with 3 decimals accuracy.

Print Composer has a separate button for "Add attribute table". I added the table using this feature, but found out, that it looses the trailing zeroes after my coordinate decimals.

For example, the original data has X-coordinate listed as "6818424.020", but Print Composers "Add attribute table" feature lists that same coordinate like this "6818424.02".

It is required for my work to display coordinates with 3-decimals accuracy. I could not find any option to add trailing zeroes at any precision, not to mention with 3 decimals.

Is there any way to change how (with what accuracy) trailing zeroes are displayed in Print Composers attribute table? I know this is possible in Arcmap (which I have no access to at the moment), but cant find this in QGIS.


texts: Filter ImageCollection to Images with non-masked coverage within AOI in Earth Engine

texts:
My goal is to Filter an ImageCollection to contain only Images with uniform coverage within a given AOI (Images that have to 'holes' / masked pixels within an AOI). I'm trying to do it by mapping the collection over ee.Reducer.allNonZero() setting isNotEmpty flag and then filtering images by corresponding flags.

This is a script I created (EE code version):

var NO2_coll = ee.ImageCollection("COPERNICUS/S5P/OFFL/L3_NO2")
                .select('tropospheric_NO2_column_number_density');

var aoi = ee.Geometry.Polygon(ee.List([
  [9.18804289906566, 45.468566496899996], [10.203754688320892, 45.468566496899996],
  [10.203754688320892,45.695567257444594], [9.18804289906566, 45.695567257444594],
  [9.18804289906566, 45.468566496899996]]));

var date_from = '2020-01-04';
var date_to = '2020-01-05';

function setEmptyFlag (image) {
  var isNotEmpty = image.reduceRegion({
    reducer: ee.Reducer.allNonZero(),
    geometry: aoi,
    scale: 7000,
  }).values().get(0);

  return image.set('isNotEmpty', isNotEmpty);
}

var coll_filt = NO2_coll.filterBounds(aoi).filter(ee.Filter.date(date_from, date_to));

var coll_with_zero_flag = coll_filt.map(setEmptyFlag);

var coll_filt_clean = coll_with_zero_flag.filterMetadata('isNotEmpty', 'equals', 1);

var single_scene = coll_filt_clean.first();

Map.addLayer(aoi);
Map.addLayer(single_scene);
Map.centerObject(aoi);

But ee.Reducer.allNonZero() always returnes '1' no matter the presence of masked pixels within AOI.

How do I fix it?


texts: Wrong number of parameters in PyQGIS tools

texts:
I want to use the "clip raster by extent" tool of PyQGIS.

The doc Clip qgis doc says that there are five parameters :

processing.runalg('gdalogr:cliprasterbyextent', input, no_data, projwin, extra, output)

But when I run it, it answers :
ALGORITHM: Clip raster by extent

INPUT <ParameterRaster>
NO_DATA <ParameterString>
PROJWIN <ParameterExtent>
RTYPE <ParameterSelection>
COMPRESS <ParameterSelection>
JPEGCOMPRESSION <ParameterNumber>
ZLEVEL <ParameterNumber>
PREDICTOR <ParameterNumber>
TILED <ParameterBoolean>
BIGTIFF <ParameterSelection>
TFW <ParameterBoolean>
EXTRA <ParameterString>
OUTPUT <OutputRaster>

RTYPE(Output raster type)

0 - Byte
1 - Int16
2 - UInt16
3 - UInt32
4 - Int32
5 - Float32
6 - Float64

COMPRESS(GeoTIFF options. Compression type:)

0 - NONE
1 - JPEG
2 - LZW
3 - PACKBITS
4 - DEFLATE

BIGTIFF(Control whether the created file is a BigTIFF or a classic TIFF)

0 -
1 - YES
2 - NO
3 - IF_NEEDED
4 - IF_SAFER

So 14 parameters, even so the doc is for 2.18 version and I have 2.18.9 version.

My questions is : How can I know how to fill this field that are not in the doc ? There are some help for rtype, compress and bigtif but not for the others.

I didn't name my topic with Clip tool name as that is not the only tool that make that. Last time,for create grids, I finally followed the advises of python console, but this time I would like to understand more generally how to solve the problem.


texts: Rounding raster pixel values to one decimal place in QGIS

texts:
I have single band raster whose values somehow were set to something like 1.399999999 instead of 1.4, 1.29999999999 instead of 1.3 (I'm not sure about the number of 9's there.), 1.2000047683 instead of 1.2. While some pixels do have a value 0.5.
I actually need all values to one decimal place like0.5, 1.2, 1.3, 1.4. For this I searched around and bumped into one answer in stackexchange. I did not find gdal.calc as suggested but found r.mapcalculator (dialog box looked similar). I then used the function round(A,1). But it did not work.

How can I round the pixel values to one decimal place ?

When I used raster calculator raster@1 * 100 the new raster formed has values like 50, 120, 130, 140 and then again used newraster@1 / 100 to obtain 0.5, 1.20000047683, 1.299999999, 1.399999999.


texts:
When I create a web app using the template it points to cos-t420:8399/arcgisoutput but I need it to be cos-t420.web.com.sg:8399/arcgisoutput
How can it be done? am using the java version.


texts: Expression for random numbers without duplicates

texts:
Using only the QGIS field calculator, how to create an expression that generates random numbers without duplicates?

rand (1,10) and randf (1,10) generate random numbers but could be repeated in a set of 10 records.

this expression generates 10 random numbers without duplicates

array_distinct (array_agg (expression: = rand (1,10))

but it's an array I need a single number for each record.


texts: Creating point from lat / lon values stored in two fields?

texts:
I'm trying to create a point from latitude & longitude values stored in two fields (lat, lon). I need to perform a transform to get the point from its source coordinate system into the coordinate system of the shape field. When I hardcode the lat/lon values the SQL statement works. When I try to concatenate the field names the SQL throws an error.

-- hardcoded values works:

INSERT INTO table_pt (objectid, shape) VALUES (1, sde.st_transform(sde.st_geometry('point (-90.012345 45.012345)', 1), 2, 3));

-- reading values from fields does not work:

INSERT INTO table_pt (objectid, shape) VALUES (1, sde.st_transform(sde.st_geometry('point ( ' || lon || ' ' || lat || ' )', 1), 2, 3));

and results in this error:

SQL Error: ORA-00984: column not allowed here

How do I create a point from lat / lon values stored in two fields?


texts: If I have an arbitrary polygon, and wish to hatch it with lines at a given angle and spacing, is there a simple low-level algorithm which will return the set of hatch lines inside the polygon?


texts: Set band value as property in an image collection Google Earth Engine

texts:
I have function which calculates monthly means from an image collection

var UK_SST = ee.ImageCollection('NOAA/CDR/SST_WHOI/V2').select('sea_surface_temperature')
    .filterDate(ee.Date('1990-01-01'), ee.Date('2019-01-01'))
    .filterBounds(ROI)

// month data for dataset

var startyear = 1990;
var endyear = 2019;

var startmonth = 1;
var endmonth = 12;

var startdate = ee.Date.fromYMD(startyear, startmonth, 1);
var enddate = ee.Date.fromYMD(endyear , endmonth, 1);

var years = ee.List.sequence(startyear, endyear);

var months = ee.List.sequence(1, 12);
var monthlySST =  ee.ImageCollection.fromImages(
  years.map(function (y) {
    return months.map(function(m) {
      var w = UK_SST.filter(ee.Filter.calendarRange(y, y, 'year'))
                .filter(ee.Filter.calendarRange(m, m, 'month'))
                .mean();
      return w.set('year', y)
          .set('month', m)
          .set('system:time_start', ee.Date.fromYMD(y, m, 1));

    });
  }).flatten()
);

I'm looking to set the band value as a property so I can export the image collection as a csv with the band values contained. However, when I try to set this

  .set("sst_mean", w.get("sea_surface_temperature"))

In this part of the function.

  return w.set('year', y)
          .set('month', m)
          .set("sst_mean", res.get("sea_surface_temperature"))
          .set('system:time_start', ee.Date.fromYMD(y, m, 1));

Nothing gets set. Any ideas?


texts: OpenLayers Style Text incorrectly aligned and moving on map resize

texts:
I have an OpenLayers map which I want to display SVG markers on with some text. I have got it to display the marker and it looks fine but when the map is resized the alignment of the text moves by a few pixels. See images below.
This is the marker style I have used:
    var markerStyle = new Style({
      image: new Icon(/** @type {olx.style.IconOptions} */({
        anchor: [0.5, 0.5],
        anchorXUnits: 'fraction',
        anchorYUnits: 'fraction',
        opacity: 0.9,
        scale: 1,
        color: '#253679',
        src: 'assets/images/map-marker.svg'
      })),
      text: new Text({
        offsetY: -5,
        text: "AA",
        scale: 1.2,
        fill: new Fill({
          color: "#fff"
        })
      })
    });

When the marker is first placed it looks like this:

If I resize the window as you can see the text visibly moves as I resize and is incorrectly aligned:

It sounds like a small problem but it is making my mapping application look unprofessional.


texts: Using a formula to set the offset of a line pattern on QGIS 2.18

texts:
I’ve been trying to fill city boundaries polygons with a 3 colors stripes pattern (red, orange, green), whose respective thickness is proportional to 3 variables, different for each city, but whose sum is always 104 (see figure for example).

I set a line pattern with 3 lines, red, orange and green. Fixing the distance to 4mm and using a field calculator expression, I succeeded in making the width of each color vary:

Variable_1/104*4    #for red
Variable_2/104*4    #for orange
Variable_3/104*4    #for green

But the expressions to set the offset does not produce any offset change:

Variable_1/104*2                                            #for red
Variable_1/104*4 + Variable_2/104*2                         #for orange
Variable_1/104*4 + Variable_2/104*4 + Variable_3/104*2      #for green

Is there something I did not set correctly?


texts: Using Viewshed Analysis in ArcGIS Desktop to Find Points Where Entire Area is Covered?

texts:
I have a DEM of an entire state and I would like to do a viewshed analysis. But instead of me using existing points to determine the viewshed, I would like to find all points where the entire state would be covered by a viewshed.

So for example, I want to find suitable locations for antenna placements. Instead of me just guessing trying to figure out which location would give me the best results, I would like the GIS to optimally find points where the entire state would be covered with antenna coverage.

I have access to all tools/extensions in ArcMap/ArcGIS Pro.

Is doing something like this possible?


texts: Using python to get extent and split coordinates into list

texts:
I have the following code:

for fc in fcs:
    desc = arcpy.Describe(fc)
    ext = desc.extent
    ext.split(" ")
    print ext

I'm getting error: AttributeError: Point: Get attribute split does not exist

I know my fcs list is okay and the extent property is working correctly. I want each fc in fcs to output a list with the coordinates like this:

[-179.168566355952, 51.2134728616834, 179.748708978836, 71.3896098800921, NaN, NaN, NaN, NaN]


texts: R scripts in QGIS : Get input vector layer path

texts:
When running R scripts in QGIS, upon loading my vector layer with ##Layer=vector it  automatically runs Layer <- st_read("path/to/shape.shp", quiet = True, stringsAsFactors = FALSE) on my input layer
This prevents me from being able to use the raw filepath later on in my scripts. Is it still possible to get the raw filepath for an input layer from an R script in QGIS?
The output I want when I type Layer is:
/path/to/myFile.shp

The output I get is:
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 581560.5 ymin: 5649045 xmax: 583123.1 ymax: 5649831
Projected CRS: NAD83(NSRS2007) / UTM zone 10N
id geometry
1 1 POLYGON ((583085.2 5649232,...

I'm not very familiar with st_read and I'm not sure if it stores the filepath or if I can retrieve it after st_read has been run


texts: Use variable in Focal Statistics Neighborhood Settings using ModelBuilder

texts:
Is there any way to use a ModelBuilder variable in the neighborhood settings area for Focal Statistics?

I would like to have this value be a parameter to the model.

I cannot type anything but integers in the boxes. For the other input boxes you can use %variable_name% as a placeholder but not for the neighborhood settings.


texts: TopologyException issue when doing intersection in Geopandas?

texts:
When I do overlay processes like intersection, I mostly get:

  :TopologyException: found non-noded intersection between LINESTRING
  (229971 4.39944e+06, 229971 4.39944e+06) and LINESTRING (229971
  4.39944e+06, 229971 4.39944e+06) at 229971.08857010119 4399438.0928708706

while during the same process in Arcgis it gets done. After trying many workarounds like setting buffer to no result most of the time(I don't even remember if it ever worked) I think that it has to do with high sensitivity that Geopandas may have to such mistakes in the topology whilst Arc just let slip.

In one intersection for example although it produced such an exception, it eventually produced a result that after checking it seemed fine.Other times it gives the exception and no result. Plus, intersection takes a long time in Geopandas. What should I do with the exception and possibly time?

Also sometimes in intersection it says:TopologicalError: This operation could not be performed. Reason: unknown

I hope you can help with this because it is very confusing.Thanks


texts: I want to use the QGIS Map Composer with Python to automatically generate maps. As a base map of these maps I want to use tiles generated with Mapnik. Currently I am serving the tiles with mod_tile just like the regular osm tiles. I know I can simply use my url of my tiles in QGIS, but the quality of the tiles is made for web maps, and I need higher resolution. So my question is, if it is possible to somehow export a georeferenced image with Mapnik to use it in the QGIS Map Composer. Or if there is an other way to create some sort of high resolution tile with Mapnik and use it in QGIS?


texts: QGIS Export Webmap with QGIS2Web only using Python Console?

texts:
I want to automatically export my QGIS map to Webmap using only pyQGIS commands. By now I just know how to run the plugin with commands on python console:

import qgis.utils
qgis2web = qgis.utils.plugins["qgis2web"]
qgis2web.run()

but I dont know how to export the Webmap.


texts: Constrain Moveable DND Dojo Popup Window to Map Element

texts:
I have a ArcGIS Enterprise (Portal) Web AppBuilder template that I've customized to allow for my pop-up windows to be dragged around the screen, largely using the example located here.

The issue I'm having, that I've been unable to resolve, is that the popup is not constrained to the map window.  I've tried, but not succeeded, in using Dojo's parentConstrainedMoveable.  I've found an example here that keeps a floating pane (with id = "searchbox") within the confines of a parent pane (with class = "container").  Essentially, I want the "map" element of my Web AppBuilder template to be the parent that tells my identify pop-up window that it can't be moved outside the confines of the map element.

Is there an easy way to modify the JavaScript in link 1 to incorporate the concepts of link 2?


texts:
Can QGIS do automated spatial ETL?  ... especially the automated point

I'm aware of tools like FME, Talend, etc. but price is a concern, so I'm looking for FOSS alternatives that provide spatial ETL that can be automated/part of a workflow

Also aware of ogr2ogr, but whilst powerful, it's seems to be just a file format converter, rather than a true ETL

An examples of the sort of ETL's we're looking at are ...

connect to back office application using its SOAP web service and extract list of lampposts, including USRN (Unique Street Reference Number)
query against LLPG database to get corresponding OS NGR for the USRN
transform various lookup code field into end-user English
remove sensitive data fields
insert into SQL Server spatial table
alert to configured end points on any issues

And these ETL's should be run on

scheduled days, times, etc.
response to external trigger


texts: I am in need of around 10000 and odd Indian latitude and longitude values for my openlayers project. I have to load them in a postgis db. so a sql file or atleast text file would be better than a link.


texts: QGIS KML layer won't show up in Google Maps API

texts:
I've trimmed US Census shapefile data to show census blocks in just three states, by editing the attribute table.

I then color coded the states based on values in the table by changing the style in qgis.

Now I want to export the new map as a KML or GeoJson layer. However each time I've tried, the layer will not show up on the Google Maps API basemap.

I believe the API code is correct, because when I use a test kml file from online, it shows.

I've tried exporting using various attributes (no symbology, feature symbology, symbol layer symbology, wsg 84,etc). I've tried hosting the file on dropbox, and my server for my website. I've tried geojson and kml and the layer simply won't show.

The file is 300+ MB. And Inspect Element shows no console errors

When I open the file in a text editor, things seem fine. The tags look closed.

I'm not sure what exactly is causing this problem. Is there an alternative to using QGIS? How can I get the Census data as a layer on my Google API project?


texts: How can I label the pie charts in QGIS so the labels are written under the diagrams? My pie-charts are representing localities: such as villages, towns, cities. So, I want that their names appear under the diagrams. I want them to look like this image


texts: Field Filter and Value Pick List for a Python Toolbox ValueTable Parameter

texts:
In an ArcGIS Python Toolbox (pyt file), how do I specify a filter for a Value Table parameter? In my value table I have two columns.  The first is a a field name derived from a selected table view and the second is a argument that should pass through with the field name.  I would like to filter the list of available fields by data type, and then define a pick list of values for the second field.

I've searched high and low, but setting a filter for a value table is not obvious, or it may not be possible.  When I attempt to specify a filter list for the field datatype in  the usual way the toolbox fails to load with an error message of "line 118, in getParameterInfoAttributeError: 'NoneType' object has no attribute 'list'".  Where line 118 is where I try to specify the filter as in the code below.

I can specify list of values to populate the table, but I really need this to be user selected.

field_map = arcpy.Parameter(
        displayName='Output Field Map',
        name='out_field_map',
        datatype='GPValueTable',
        parameterType='Required',
        direction='Input')
field_map.columns = [['Field', 'Field Name'], ['GPString', 'NVEL Attribute']]
field_map.parameterDependencies = [trees_table.name]
#specifying a filter list fails
#field_map.filter.list = ['Short', 'Long', 'Float', 'Double']


texts: Show popup on the map for the selected item in the side list using the mapbox gl js

texts:
I followed this tutorial "Filter features within map view" and managed to adapt it to my data, which are polygons instead of points.
However, the popups on the side menu are not displayed, when I place the mouse over one of the locations, the popup does not appear on the map indicating the location as it should.
The problem I believe is that because it is a polygon and not a point (a pair of coordinates) the code cannot display the location, the error in the console is:
"Uncaught Error:` LngLatLike` argument must be specified as a LngLat instance, an object {lng: <lng>, lat: <lat>}, an object {lon: <lng>, lat: <lat>}, or an array of [<lng>, <lat>] "

In the map popup I changed it from ".setLngLat (feature.geometry.coordinates)"
to ".setLngLat (e.lngLat)" and it worked, but not in the popup on the side menu.
I'm still starting to learn javascript and I couldn't find a solution, I thought about extracting the centroid from the polygon, but I don't know how to do it.
the code that generates the popup:
var prop = feature.properties;
var item = document.createElement('a');
item.href = prop.wikipedia;
item.target = '_blank';
item.textContent = prop.NM_UF + ' (' + prop.SIGLA_UF + ')';
item.addEventListener('mouseover', function () {
    // Highlight corresponding feature on the map
    popup
        .setLngLat(feature.geometry.coordinates)
        .setText(
            feature.properties.NM_UF +
            ' (' +
            feature.properties.SIGLA_UF +
            ')'
        )
        .addTo(map);
});
listingEl.appendChild(item);
    });

Here is a codepen like the code:
https://codepen.io/hugonbgg/pen/PoNwWVj and the debug: https://codepen.io/hugonbgg/debug/PoNwWVj/


texts: Export QGIS Print Composer PDF with searchable text labels

texts:
I have point features with text labels on my map. When I export as PDF, I can search elements that were constructed in the print composer, like the legend, scale bar, etc. However, labels of features on the map are not exporting as searchable text.

I've tried exporting in SVG and PDF. I'm not exporting as raster. Just to be sure I've isolated the problem I've also tried to export a map with just my text labels and get the same issue (text labels not searchable).


texts: Get closest road-segment to segmentized LineString points

texts:
So I have a table that stores a bunch of traffic counts for various road segments across the country. I'm regularly getting real-time location data from vehicles that are moving around a given region and I'm trying to match the segments that correlate to the path it traveled. In order to fill in some gaps in where the vehicle has traveled, I'm using Mapbox's Matching API to snap to roads. This is good for getting a general route that the vehicle may have traveled, and is returned to us in an encoded polyline. In PostGIS I'm able to decode this polyline, use ST_Segmentize to fill in any gaps between 'snapped-to-roads' points (to avoid missing any small segments) and loop over these points to find the closest segment geometry they match to. This is fine, and it technically works, but it's really slow. I'm hoping to find a faster/better way of doing this.

Edit: In case anyone was wondering, the reason we break the segments into distances of no-greater-than .000032 is to make sure we capture even the smallest road-segments when breaking up the encoded polyline.

This is what I have so far (apologies for bad formatting):

    CREATE OR REPLACE FUNCTION findASegment(polyline text)
      returns table(t_index float4, t_xdsegid TEXT, t_distance DOUBLE PRECISION) as $$
    DECLARE
      node     text;
      iterator float4 := 1;
    BEGIN

      CREATE TEMPORARY TABLE temp_table (
        t_index    float4,
        t_xdsegid  TEXT,
        t_distance DOUBLE PRECISION
      ) ON COMMIT DROP;

      FOR node IN (SELECT ST_AsText((dp).geom) as node
                   FROM (SELECT (st_dumppoints(
                         ST_Segmentize(st_linefromencodedpolyline(polyline),
     .000032))) AS dp) AS foo)

      LOOP

        INSERT INTO temp_table (t_index, t_xdsegid, t_distance)
    select iterator as t_iterator, xdsegid, distance
        FROM (SELECT bar.xdsegid, st_distance(geom, st_geomfromtext(node)) as distance from alabama as bar) as baz
        ORDER BY distance ASC
        LIMIT 1;
        iterator := iterator + 1;
      END LOOP;

      RETURN query SELECT t.t_index, t.t_xdsegid, t.t_distance FROM temp_table t;

    END; $$

    language plpgsql;

    SELECT *
    FROM (SELECT distinct on (t_xdsegid) t_xdsegid, t_index
          FROM findASegment('{hhkEzz|pOGQq@cC')
          order by t_xdsegid, t_index) as p
    ORDER BY p.t_index;

Here are some pictures as per-request:

So we initially are given two sets of lat/long coordinates. We use MapboxMatching to fill in points on a road that a vehicle could have traveled to get from point A to point B. Mapbox returns us an encoded polyline that looks like this on a map:

But, the issue with this, is that sometimes those points in the encoded polyline are too far apart. They may miss some of the smaller road-segments along that path. So we segmentize it to interpolate additional points into it.

Now it looks more like this:

Then I want to get all the possible segments it passed through from the geometry table. Each row in the table looks like this:

The reason we use distance to get the nearest segment from the table is because the snapped-to-roads geometries don't necessarily match up to the stored geometries. The geometries in the database may veer off the road slightly, making it to where contains isn't viable (to my knowledge). Hope this helps


texts: How to get nearby tiles from a location?

texts:
I want to get the tile in which a particular location lie and also its nearby tiles.
Is there any google map api which could return this information.


texts: 3-band RGB raster with 0 value to null -ArcGIS Pro

texts:
I've got a 3-band RGB raster layer where I'd like to convert the value of a pixel to null if b1, b2 and b3 are 0.
Any idea how this can be done in ArcGIS Pro?


texts: Mapfish 3 print service, GeoServer 2.15, and Jasper Report set up

texts:
I've been going through this article, this Mapfish documentation and this Geoserver doc (apparently out of date) in attempt to get the Mapfish3 print service going in order to eventually use in my web app (using GeoExt3 and Openlayers 4, but this probably is extraneous in terms of this post).

I was hopeful when I saw this question, but it is for Mapfish2, which is quite different from Mapfish3. Also note that I am not trying to use the Geoserver print plugin mentioned in the Geoserver doc, but I referred to it as it is how I deduced to copy the jrxml and yaml files (potentially incorrect).

I have done the following:

Deployed the test print module to my tomcat web server
Created a simple jrxml file and a simple config.yaml file (copied into the GeoServer2.15.1/data_dir/printing folder (contents reflect the mentioned article's models exactly except they point to a wms on my hosted geoserver),
Written up specs to be pasted into the test print input window

Upon pressing Create and Get Print a pdf is generated. However, it is a blank white page with the title "Map" on it, which does not reflect any of the three files I've established (e.g. default title should be "Pays" and map is not displayed). I assume I have not put these files in the correct location, or am missing something else.

According to this page of the Mapfish documentation under reportTemplate:

  The path to the report template used to render each row of the data.
  This is only required if a subreport needs to be compiled and is
  referenced in the containing report's detail section. The path should
  be relative to the configuration directory.

Is the configuration directory the geoserver/data_dir as I've assumed?

My config.yaml:

throwErrorOnExtraParameters: true
transparentTileErrorColor: "rgba(78, 78, 255, 125)"
opaqueTileErrorColor: "rgba(255, 155, 155, 0)"

templates:
  A4 landscape: !templates
    reportTemplate: Map_A4_Landscape.jrxml
    attributes:
      title: !string
        default: "Pays"
      map: !Map
        maxDpi: 400
        width: 780
        height: 330
      processors:
      - !reportBuilder
        directory: '.'
      - !createMap

The specs copied into the test app:

{
  "layout": "A4 landscape",
  "outputFormat": "pdf",
  "attributes": {
    "map": {
      "projection": "EPSG: 3297",
      "dpi": 72,
      "rotation":0,
      "center": [-149.3, -17.7],
      "scale": 130000000,
      "layers": [
        {
          "baseURL": "http://xxx.xx.xxx:8080/geoserver/test/wms",
          "opacity": 1,
          "type": "wms",
          "layers": ["workspace:layername"],
          "imageFormat": "image/png",
          "styles": ["stylename"],
          "customParams": {
            "TRANSPARENT": "true"
          }
        }
      ]
    }
  }
}


texts:
I have two DSM rasters, one is 10km radius photogrammetric data, the other is 2km radius LiDAR data, both centred around the same point.

I need to combine these together, essentially inserting the 2km LiDAR raster into the other and removing that portion of the photogrammetric data, the idea being that the LiDAR is much more detailed in the area closest to the centre and the site that I'm working on.

I've resampled the photogrammetric raster to match the LiDAR so the cell sizes are the same.

What is the best way to go about doing this? I think it's using 'Mosaic to New Raster' but can't get my head around whether to use FIRST or LAST


texts: Dynamic Display doesn't respect my annotations when my layers are in a BasemapLayer

texts:
I have 2 annotated layers which I load in a BaseMapLayerClass.  When the map is static the annotations are visible.  When I enable DynamicMap I lose all annotations.  I would like a way to annotate layers in a BaseMapLayerClass while dynamic is enabled.

An interesting note about the problem.  If I don't have my layers in a BaseMapLayerClass, dynamic display does not affect my annotations.

However I would like them to be in a BaseMapLayerClass because would make it much easier for me to save layers and restore them across sessions.

I'm using ArcEngine C# VS2010

I think the problem is somehow related to extents.

This post: http://briannoyle.wordpress.com/2008/03/18/esri-dynamic-displaysome-limitations/

Led me to the limitations page
http://edndoc.esri.com/arcobjects/9.2/NET/0b062497-a3fc-4c56-a688-0fb9b2fcb4ea.htm

which states :

Labels will be cut off on custom layers that don't return the correct extent

And actually it's weird.  Sometimes the annotations will draw and then I'll move the map up a few pixels with the hand tool, and then they will no longer draw.  Then move it back a few pixels and it will draw again.

It's helpful there is a consistency there, but I haven't a clue how to implement a workaround.


texts: Adding bike infrastructure data to OpenStreetMap?

texts:
Is there a best practices guide for adding bike infrastructure data to OSM? Not simply how to add a bike path, but a list of best practices for actually fleshing out sometimes disjointed bike infrastructure information for a city.

Some big issues:
How to split road segments to mark on-road cycling infrastructure and when it's appropriate and when it's not. How to deal with streets with a separated cyclepaths that becomes a bike lane and then sharrows.
How to mark separated Sidepaths, shared use bike / pedestrian paths, paths through parks, different surface types (such as unpaved rail trails)
How to denote official US federal / state / county / local bike routes.
When is it appropriate to add mountain bike trails and how to mark them and differentiate them from commuter trails.

I'm looking for any official or unofficial guides to point people to that will help map editors fill in this information.  If this is more suited to multiple questions, I'm more than willing to have this closed and reopen separate questions.


texts: Using linear referencing for real time traffic in ArcGIS Network Analyst?

texts:
I'm using real time customized traffic for Network Analysis in ArcGIS 10.2.2.

The traffic display as segment based. In other words, when i want to show free flow for a road i have to assign one or more segment(s) as free flow.  But i want to assign a part of a segment as free flow not entire segment.

Is it possible to assign traffic as part of a segment?


texts: Joined field not recognized as number in QGIS 2.4

texts:
I have a weird problem. I've loaded in QGIS a point shapefile and a CSV table. Then I performed a join between the 2 layers and all works fine.

But then the joined columns (made of integers or real numbers, it doesn't matter) are not recognized as numbers but as strings.

Even if I save the shapefile (original + join fields) and I reload it in QGIS the columns are always seen as strings.

BTW: I saved the CSV as a DBF table but same behavior.

Am I missing something?


texts: Error exporting shapefile from QGIS into PostGIS?

texts:
I know this error below has been posted a couple times on this forum, but none of the solutions seem to work for me.

I am trying to export a shapefile from QGIS to a PostGIS server using DB Manager.

If I check "create single-part geometries instead of multi-part" it works but obviously not as a multi-part geometry type.  Everything uploads fine in PG Admin III, but I want to be able to do everything from QGIS.

If I try to edit a layer already on the PostGIS server I get the same error (second error message below).  What's the source of this problem?

If it's the shapefile why does it work in PG Admin III?

Also if I sometimes take the same file as a shapefile, I can upload it through DB Manager and work with it fine.  No problems.  Makes all sorts of multipolygons.

But if I try the same thing with a file uploaded through PG Admin III an they work on th layer through DB Manager it doesn't work.

Note: some layers nothing works.

Error 7
Feature write errors:
Creation error for features from #0 to #0. Provider errors was:
PostGIS error while adding features: ERROR:  Geometry type (LineString)
does not match column type (MultiLineString)`

Errors: ERROR: 1 feature(s) not added

Provider errors:
PostGIS error while adding features: ERROR:  Geometry type (LineString)
does not match column type (MultiLineString)

Key Question: Why does the exportation of a shapefile from QGIS to PostGIS using DB Manager create errors with multi-part polygons and polylines?


texts: QGIS2web page does not display curved labels

texts:
I have designed a non-real map, using QGIS. Until now, my map was very basic (a few layers with some polygons, dots, and lines). This map is a "decal" from a raster image.

I recently add some custom labels, which are displayed into a certain zoom range, and also some labels (using invisible dots or lines, if I need curved labels).

Here is how a part of my map looks in QGIS :

Here is what the qgis2web export give me for the same area :

How can I export the map to the web in the most accurate way possible?


texts: 90 angle between coordinate breaks when applied to world

texts:
I have a point A and direction DIR and a distant DIST

Let's use

A: `longitude: 139.33435143, latitude: 35.5385006, altitude: 0`
DIR : `longitude: 139.33547146, latitude: 35.54034106, altitude: 0`
DIST : `500`

I would like using those value, to, create two point B/C so that BAC is a 90degree angle. and that AB is giving the starting vector on which apply the rotation.

Basically like this :

Sandcastle

The problem I have is the following, event tho, my long/lat that I calculate are forming a 90 angle :

origin = A.longitude, A.latitude
xAxis = A.longitude - DIR.longitude, A.latitude - DIR.Latitude
yAxis = cos(Math.PI / 2) * xAxis.longitude - sin(Math.PI /2) * xAxis.latitude,
         sin(Math.PI /2) * xAxis.longitude  - cos(Math.PI /2) * xAxis.latitude

So this would give me the following which have a 90degee angle when drawn.

origin = 139.33435143,35.5385006
xAxis = 139.33547146,35.54034106
yAxis = 139.33251097,35.53962063000001

When I apply those coordinate and convert them to Cartesian3 in cesium, the resulting angle is not 90degree anymore but something smaller (as you can see in the link)

Since I believe converting to Cartesian3 apply the circular shape of the globe to the coordinate is there some additional conversion I am missing ?


texts: how do I map from a mapbox map or tile to the lat lon extents of the map?

texts:
A web app I use displays map tiles from mapbox and then uses leaflet to overlay that map with an svg path representing a route.

The problem is that the area given to the map is very tiny, it's short and wide.
I'd like to create an extension that takes the map, examines the route, and then creates a kml file representing the route.

How do I map from mapbox tiles to their real world lat and lon?


texts: Applying SLD to WMS layer with OpenLayers

texts:
I am attempting to adjust the symbology of a loaded WMS raster layer (img file through MapServer).  It is the NLCD dataset.  The WMS layer loads find with

nlcd = new OpenLayers.Layer.WMS(
    'nlcd',
    'http://server /fswms/html/rlayers',
    {
        layers: 'nlcd-2006', transparent: true,
        projection: new OpenLayers.Projection("EPSG:4326")
    },
    {
        isBaseLayer: false
    }
);

And I have a SLD file that I have created in which I am attempting to make all of the pixels with the value of 11 opacity of 0.

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
    xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
    xmlns="http://www.opengis.net/sld"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:xlink="http://www.w3.org/1999/xlink"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<NamedLayer>
    <Name>nlcd</Name>
    <UserStyle>
      <Title>ncld test</Title>
      <FeatureTypeStyle>
        <Rule>
          <RasterSymbolizer>
            <ColorMap>
              <ColorMapEntry opacity="0" color="#000000" quantity="11" />
            </ColorMap>
          </RasterSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

The WMS layer ends up being map.layers[3] in my test application so I attempting to apply the SLD by the following code but it is not working.

//this is in the init:
OpenLayers.Request.GET({
    url: "sldcode.sld",
    success: complete
});

//function that attempts to apply SLD
function complete(req) {
    mySLD = format.read(req.responseXML || req.responseText);
    map.layers[3].mergeNewParams({sld_body:mySLD});
}

Here a link to the dev site: http://rain.nemac.org/~derek/ol/sldTest.html
Any help would be greatly appreciated.


texts: Is it possible to layer stack all 13 Sentinel 2 bands or can I only layer stack those that are the same resolution (so for example the blue, green red and NIR bands). I have a very large area of interest so will have many sentinel images. In order to exploit more band combinations, it would be great to layer stack all bands for each image and then create a mosaic, but I am not sure if this is possible? Do I need to rescale the coarser bands prior to layer stacking?


texts: What does it mean when in certain publications they write population density (1000 people/km2)? Does it mean that they calculate the population density by dividing the number of people by the area of land in square meter and then they divide the result by 1000?


texts:
I have a Postgres/PostGIS DB with two geometries in a single table. I want to export them as a GeoJSON GeometryCollection (see here) but I can't seem to get the GEOMETRY_AS_COLLECTION layer create option to work properly in ogr2ogr.

ogr2ogr -skipfailures -progress -s_srs EPSG:4326 -t_srs EPSG:4326  -f "GeoJSON" .\myFile.json PG:"host=server user=postgres dbname=name password=password" -sql "select * from table" -nln test -nlt GeometryCollection -lco GEOMETRY_AS_COLLECTION=YES

How do I get this working?


texts:
I am having issues reading spatial tables in QGIS with a select-only "reader" user.  I can create the reader user to access regular postgres tables like this:

GRANT USAGE ON SCHEMA test TO reader;

ALTER DEFAULT PRIVILEGES IN SCHEMA test
    GRANT SELECT ON TABLES
    TO reader;

What do I need to add to allow the user to view tables in QGIS?


texts: Is there a website that will let me see add markers and bounding boxes on top of a street map?

texts: Is there a website out there that will let me markers (lat/lon) and bounding boxes (minimum / maximum lat/lon) and let me see it overlayed on top of Google Maps (or OSM, or Bing, etc)? I'm not really looking to create my own webpage for this... I'm hoping someone's implemented something like this for a general use case or a demo and it's online somewhere.


texts: Checking which application Python toolbox was launched from?

texts:
Esri's documentation on Python toolboxes states, 'When the script is run as a script tool, arcpy is fully aware of the application (such as ArcMap) it was called from.'

How can I check which application launched the Python toolbox tool?

I have a parameter allowing the user to use the display extent as the area of interest, but it only works in ArcMap. I would like to disable the Boolean parameter if launched from ArcCatalog.


texts: Calculating the greatest length of a polygon using QGIS

texts:
I'm new to QGIS and I need to know the easiest way to calculate the greatest length of a polygon. I already know the "rule" tool but I want to do it automatically because I have a huge amount of polygons to process. I only know the basics for the moment.
My polygons are watersheds and they look like this:

or even that:


texts: Odd performance of PostGIS Select Query with ST_AsText?

texts:
I have this query:

    select ST_AsText(geo) as wkt, city, street_number, street_name, street_type,
           street_direction, street_unit
           from parcels_union
           order by id
           LIMIT 1000 OFFSET 18000

I takes 5500ms to run.

This query without the ST_Astext runs in 22ms. So I would expect it is taking some time to AsText my geometries.

Now here is the weird part. This query also runs in 5500ms even though it is returning 5 times as many rows:

    select ST_AsText(geo) as wkt, city, street_number, street_name, street_type,
           street_direction, street_unit
           from parcels_union
           order by id
           LIMIT 5000 OFFSET 18000

This tells me there is a fixed overhead to using AsText which I do not find particularly odd in itself, but 5.5 seconds?

Any ideas on how to bypass this overhead or otherwise optimize this?


texts: Consume WFS Stored Query with Parameters using ArcMap

texts:
I am working on providing data from our application in a WFS friendly format. I have set up a stored query with parameters that is to be used to return a feature collection. My stored query is within our own web application, not through ArcGIS. A simple call to my stored query with just one parameter looks like:

http://localhost/AWQMS/api/wfs?service=wfs&request=GetFeature&storedQuery_id=GetFeaturesByParameters&OrganizationIdentifiersCsv=BISHOP_PAIUTE_TRIBE_WQX

I want to be able to consume this using ArcMap. So far in ArcMap I can set up an Interoperability connection with my url http://localhost/awqms/api/wfs and my capabilities seems to work alright...at least my Feature Types list populates correctly in the Constraints section of the connection's parameters window.

I only have a single FeatureType and was planning to rely upon the stored query's parameters to narrow down which features are actually returned.

Is this a reasonable approach or should I be thinking differently about this?
How can I consume my stored query with parameters using ArcMap?


texts: How to change PostgreSQL layer encoding for QGIS Server

texts:
I'm actually stuck trying to update encoding of my PostgreSQL layers to be able to set my project as QGIS server project.

In the layers properties, I can't change the layer's encoding.

Any hints on how to proceed?


texts: Is there a way to filter results in the attribute table of a point shapefile by the extent of a second polygon shapefile in QGIS?

texts:
I am currently trying to temporarly filter data with in the attribute table to a confined area defined by a second polygon shapefile. The data set to be filtered consists of 100s of points that extent outside the areas I am intrested in. I have multiple seperate area that I am examining so would prefer not to clip data multiple times for each individual area. Is there a way to simply filter by an expression in the field calculator using the second polygon as the area extent?
I am using QGIS v.3.14.0


texts: QGIS's "Points in Polygon" counts zero in every census tract

texts:
I am mapping the location of incidents of asthma using QGIS. I would like to count how many points are in census tracts to better understand the areas where asthma occurs most frequently. My map canvas presently displays the points of asthma from a CSV file, as well as census tract polygons from U.S. Census Bureau shapefiles.

I should be able to run the Points in Polygon function to help me. The function seems to work as I enter my CSV as the points, and the census tracts as the polygon. The function then returns a shapefile with an added column displaying "PNTCNT." Great so far. But, the attributes table indicates all zeros under the "PNTCNT" column. (This is obviously incorrect as I can plainly see that there are a many points in the census tract polygons.)

I tried saving the CSV file as a shapefile, and running the function again, but the result is the same. I'm at a loss.

Why is Points in Polygon giving me all zeros?

Here's the process I've used. Maybe someone can tell me where this is going wrong.

First, I add census tracts as shapefiles to my map. Second, I add a CSV that includes latitude and longitude for about 400 points. The CRS is set to WGS84 (EPSG:4326) for the points, and WGS84 / Pseudo Mercator (EPSG:3857) for the shapefile. Everything looks great on the map canvas.

Next, I select the points in polygon under the vector analysis tools. I select the shapefile under "input polygon vector layer." I select the CSV file under "input point vector layer." I input a file name and then hit okay. QGIS adds the new layer to my map canvas. However, when I open the attributes table the PNTCNT column reads all zeros. I have also tried the same process described above with the CRS set for both the shapefile and the CSV at WGS84 (EPSG:4326) (instead of using Pseudo Mercator for the shapefile).


texts:
I'm relatively new to service area analysis. I use ArcMap 10.2

I made the following travel time contour maps by performing service area analysis.I made the first one first. It looks good. Then I add more roadways to make the second map. As you see, for unknown reason, the lower half of the contour shrinks (picture 3). I can't figure out why.

The dataset I used for the second map is built based on the dataset used to make the first map. I just added a little more local roads. So both maps basically use the same data. They calculate the travel time towards the "star".

I rebuilt the network when adding more roadways. Actually, the network I built to do service area analysis is rough. There is no turning, direction (so no loop back problem?), elevation restrictions. The only attribute I used from the dataset is travel time calculated based on links.I used network analyst to identify disconnected links and corrected them. And I tried topology check, it seems there is no self-overlap lines. Not sure what else I can check with, there are only a few options for topology check of lines.

I'm seeking for suggestions for potential reasons caused the failure of the second map.
Thanks!


texts:
I am trying to quantify the number of months required to reach a specific value at the grid level. For this, the monthly raster values were then ranked, and the time length was defined as the number of months required to reach the total value.

I have built a function using raster::overlay to solve this, but I don't know how to include the order
of the pixels in function.

Any help here much appreciated.

Here is my code.

library(raster)

# Rasters at monthly basis (represent 3 months)
r <- raster(nr=100, nc=100)
r[] = sample(ncell(r))

r2 <- r
r2[] = sample(ncell(r))

r3 <- r
r3[] = sample(ncell(r))

# Total value raster
r4 <- mean(r, r2, r3)*0.8

# Rank pixels by month
rm <- stack(r,r2,r3)
rm_ord <- calc(r_m,function(x, ...) order(x,decreasing=T) )

# Function for count n months to reach the total value
count_m <- function(x1,x2,x3,x4) {
  ifelse(x1 >= x4, 1,
   ifelse(x1 + x2 >= x4, 2,
      ifelse(x1 + x2 + x3 >=x4, 3, 0)))
}

# Compute time length in n months
s <- stack(rm, r4)
test <- overlay(s, fun=count_m)

plot(test)


texts: pgsql2osm? How to generate osm files from osm2pgsql database?

texts:
I  want a method to export osm files from a Postgres database. I already have an osm2pgsql database, an apidb database, and a nominatim database in which I have imported osm files. Now I need to export osm files from a  database with an osm2pgsql schema and use it with the nominatim and api-db(openstreetmap-website) database.
I know that the osm2pgsql is a tool for populating a PostGIS database that allows us to run the Mapnik rendering software, Nominatim search software, or our geospatial queries and have a different database table structure to that required by the openstreetmap-website rails app.
But, I need to create for use in future imports using osmosis and osm2pgsql tools and for backup purposes.
Can anyone help me with how to do this?


texts:
I am trying to combine two different raster files using a loop, with the following script.

import  arcpy
import os

thectr = 0
ws = arcpy.env.workspace = r'D:\BRB Snow Cover\Sample\TerAqua'
Terra_Rasters = arcpy.ListRasters("MOD*")
Aqua_Rasters = arcpy.ListRasters("MYD*")
output = "D:\\BRB Snow Cover\\Sample\\output""\\"+ str(thectr)+'.tif'
matches = 0
for item1 in Terra_Rasters:
    for item2 in Aqua_Rasters:
        print item1
        print item2
        if item1.endswith(item2[32:51]):
            print "Match"
            arcpy.gp.Con_sa(item1,item2,output,item1,"\"Value\" =8")
        else:
            print "No"
            matches = matches+1

I need to use a loop for the output files name, I faced the following error.

  Runtime error  Traceback (most recent call last):   File "",
  line 7, in    File "c:\program files
  (x86)\arcgis\desktop10.5\arcpy\arcpy\geoprocessing_base.py", line
  510, in      return lambda *args: val(*gp_fixargs(args, True))
  ExecuteError: ERROR 000872: Output raster: Dataset D:\BRB Snow
  Cover\Sample\output0.tif already exists and cannot be overwritten
  since the Overwrite existing datasets option is disabled.


texts: Not able to parse Latitude & longitude to ogr_geometry

texts:
I have a table (Table1) of locations with Latitude and Longitude values.
Loc_Nm || geocode_latitude || geocode_longitude
   A   ||    43.43044      ||    -80.0910793
   W   ||    43.9114523    ||    -80.0911329
   X   ||    43.9209011    ||    -80.091323
   Y   ||    43.9083203    ||    -80.0913039
   Z   ||    43.914577     ||    -80.091198

I downloaded shapefiles for "Federal Electoral Districts - Canada 2015 (English)" from the link given below.
https://open.canada.ca/data/en/dataset/737be5ea-27cf-48a3-91d6-e835f11834b0
The, Shaepfile was in "3347 - NAD83_Statistics_Canada_Lambert" projection and as I dealing with Logitudes and Latitudes, I converted it to "4326 - GCS_WGS_1984" projection. I achieved this through ogr2ogr command:
ogr2ogr -f "ESRI Shapefile" -t_srs EPSG:4326 -s_srs EPSG:3347 FED_CA_2_2_ENG_Modified.shp FED_CA_2_2_ENG.shp

Then, I loaded the resulting Shapefile to Local SQL Server DB using the ogrogr command:
ogr2ogr -f MSSQLSpatial "MSSQL:server=RBANSAL\SQLEXPRESS;database=MySpatial;trusted_connection=yes" C:\Users\rbansal\Desktop\ShapeFile\FED_CA_2_2_ENG_Modified.shp -nln "Spatial_Data" -progress

Now, I need to group the Table1 locations as per the Federal Ridings loaded into the "Spatial_Data" Table above. I am using below given query for that purpose. It runs for quite a while but, eventually completes with 0 records.
select * from MySpatial.dbo.Spatial_Data s, MySpatial.dbo.geocodes_raw g
where ogr_geometry.STContains(geometry::Parse('POINT(' + CAST(g.geocode_longitude AS VARCHAR(20)) + ' ' + CAST(g.geocode_latitude AS VARCHAR(20)) + ')'))=1;

It runs for quite a while but, eventually completes with 0 records.
What is possibly wrong here?


texts: How to split and orthomosaic into a grid / index polygon

texts:
I would like to know is there a way to split a large Orthomosaic of an area into equal grids of 10 meters each?

Sofar i have tried the " split raster" tool under " raster processing "[, but this does not work, it indicates that the output was successful, but the folder is empty.

I am using ArcGIS desktop 10.5.


texts: Extract polygon data for each neighbourhood from shapefile

texts:
I have a dataset that has four types of files with the following file extensions - .DBF, .PRJ, .SHP, .SHX

My goal is to get all the polygon coordinates for each neighbourhood and draw the boundary outline on the Google map when a user is on a specific neighbourhood page. Currently the data source for the site is a MySQL database.

I have tried converting the shapefile to GeoJSON but then parsing that file is time taking.

I have also tried using tools such as QGIS and shp2sql etc.

I would like to know what is the best way to get the main attributes from the dataset along with each neighbourhood's polygon coordinates?


texts:
I would like to calculate the total number of events within each feature of a polygon using a dbf file and a polygon shapefile.

Is there an easy way to automate this so I can then produce charts per feature area?

Thanks,

Haskett


texts: How to find out which polygon contains the most points using Python OGR?

texts:
I have two shapefiles: point (railway stations) and polygon (regions). My task is to write a script using OGR module in Python that can determine in which region the biggest amount of railway stations is.

I don't have any problem loading shp and reading its content, it's just the mechanics of solving the task that I can't seem to figure out.

I thought if I iterate over the regions (kind of abusing SetAttributeFilter), get their geometry, then in another loop get geometry of the points and finally do intersect (or maybe within, I am not sure, seemed to un/work the same), I get the result but I must be doing something wrong. Funny thing is, the first result I get is correct but following on it doesn't work as expected.

You can check the crucial part of my code (until here it works fine), I am a beginner so it may be ugly.

filters = []

for code in codes:
    filter1 = str("KOD_KRAJ='" + code + "'")
    filter2 = str('"'+filter1+'"')
    filters.append(filter2)

stations = 0

for filter in filters:
    regions.SetAttributeFilter(filter)

    for i in range(0,numberRe):
        regionFeat = regions.GetNextFeature()
        regionGeometry = regionFeat.GetGeometryRef()

        for r in range(0,numberRa):
            railwayFeature = railway.GetFeature(r)
            railwayGeometry = railwayFeature.GetGeometryRef()

            if railwayGeometry.Within(regionGeometry):
                stations = stations + 1
        print filter, str(stations)

Note: KOD_KRAJ is the header in shapefile's attribute table and variable "codes" is a list containing its instances (administrative code of each region). SetAttributeFilter should be used like this ("NAME='something'") which is why the weird first loop.

I just want to get it working to the stage when I get the numbers for each region, after that it will not be hard for me to do the rest.


texts:
I need to keep a record of the OpenStreetMap objects that I have edited. But I can't see any IDs attached to road segments & the URL doesn't change when I select them.

How can I uniquely identify road segments?


texts:
I'm using QGIS raster calculator to convert a panchromatic image of digital numbers (DNs) to radiance. The equation I need to use goes:
0.923*"15JAN12-P2AS"*(7.973540e-02/2.896000e-01)+(-1.700)

The raster DN values range from 0 to 1143, but when I apply the raster calculation above the result gives an image of NaNs and statistics as follows:

STATISTICS_MAXIMUM=-1.7000000476837
STATISTICS_MEAN=-1.7000000476837
STATISTICS_MINIMUM=-1.7000000476837
STATISTICS_STDDEV=1.2508972339943e-05

Any ideas as to why this calculation fails?


texts: How to create points grid on a 3d surface with geometry polygonz in postgis?

texts: I am using 3dCityDatabase that imports CityGML buildings to postgresql as polygonz surfaces. In order to carry out an analysis I need to create points on these building surfaces. How can I do that? Which postgis functions or sql statements should I use for this purpose?


texts: cross platform Mobile tracking app to report to custom server?

texts:
Does any one know of IOS, Andriod apps for asset tracking that I can point to my server?

I can then write some code to intercept the messages from the mobile app.

ideally a solution which works on multiple platforms

I have a tracking web app which current works with dedicated GPS hardware, I now want to add mobile devices to the web app.
instead of developing my own tracking software on every mobile platform, i'm wondering if there are any apps already in existence
I need the mobile app to to send gps tracking positions to a nominated IP address, it can be very simple, and all it will do is send a message to the IP with the Lat and LON every x secs

Is such a thing already out there?


texts:
-A layer that contains multiple features
- I want to draw new features with new OpenLayers.Control.DrawFeature and to add them to my layer.

   var control = new OpenLayers.Control.DrawFeature(layer,OpenLayers.Handler.Polygon);

control.events.register('featureadded', control, function(f) {
            f.feature.state = OpenLayers.State.INSERT;
            f.feature.attributes.Geo_ID = 794;
            saveStrategy.save();

       });

When the feature is added, a xml file is sent to the geoserver with the following content:

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<wfs:Insert>
<feature:Graeber xmlns:feature="myserver.de/Ahaus_neu">
<feature:the_geom>
<gml:Polygon xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:31466">
<gml:exterior>
<gml:LinearRing>
<gml:posList>363183.0625 5770740.875 363192.0625 5770755.375 363177.5625    5770764.375 363168.5625 5770749.875 363183.0625             5770740.875</gml:posList>
</gml:LinearRing>
</gml:exterior>
</gml:Polygon>
</feature:the_geom>
<feature:Geo_ID>794</feature:Geo_ID>
</feature:Graeber>
</wfs:Insert>
</wfs:Transaction>

The coordinates are not saved in the geoserver, only the feature:Geo_ID Element.

If I choose an existing Element to insert in the geoserver Database instead:

  control.events.register('featureadded', control, function(f) {

            layer.features[0].state = OpenLayers.State.INSERT;

            f.feature.state = OpenLayers.State.INSERT;
            f.feature.attributes.Geo_ID = 795;

            saveStrategy.save();

       });

I get the following XML in the Request:

<wfs:Insert>
<feature:Graeber xmlns:feature="myserver.de/Ahaus_neu" fid="Graeber.1">
<feature:the_geom>
<gml:MultiSurface xmlns:gml="http://www.opengis.net/gml"     srsName="EPSG:31466">
<gml:surfaceMember>
<gml:Polygon>
 <gml:exterior>
<gml:LinearRing>
<gml:posList>363273.29099999997 5770863.855 363273.177000001     5770864.300999998 363274.1561700017 5770864.5512810005 363274.2701699999 5770864.105281 363273.29099999997 5770863.855</gml:posList>
</gml:LinearRing>
</gml:exterior>
</gml:Polygon>
</gml:surfaceMember>
</gml:MultiSurface>
</feature:the_geom>
<feature:Geo_ID>795</feature:Geo_ID>
</feature:Graeber>
</wfs:Insert>

And the coordinates DO get inserted in the geoserver DB.

I see one major difference between these two XMLs:

the insertion  XML of an existing Element of the layer contains two more tags:

<gml:MultiSurface xmlns:gml="http://www.opengis.net/gml"     srsName="EPSG:31466"> and
<gml:surfaceMember>

How do I make the application know that the features added by drawing a polygon are to be treated the same way and put in the same format as the features that already exist in the Database?


texts:
I made a plugin and when I launch it it gives me this problem.

How can I fix this?


texts: Search Cursor Return Field Name If All Rows Contain a Value

texts:
I'm trying to print the field name to a list if all rows in that field are not NULL. Here is the logic I have so far.

fiberCable = r'Orlando\Orlando_FIM_prep\FIBERCABLE'
list = []
with arcpy.da.SearchCursor(fiberCable, "inventory_status_code") as cursor:
     for row in cursor:
         if row[0] is not None:
             list.append(str(row[0]))

Using the logic on the table above it prints out this list:

['Preliminary Designed', 'DesignComplete', 'DesignComplete', 'DesignComplete', 'DesignComplete']

How do I get the search cursor to look at the row of values in the field and if there is a value present in all rows, print the field name, which in this case would be  'inventory_status_code'.

The outcome I would be looking for here is:
['inventory_status_code']


texts:
Does anybody know how to create a new folder with arcobjects?
This is what I want to do: insert a text in a TextBox and, when I run the aplication, it will create a folder with the name wroten in the TextBox.

Thanks,
Alex


texts:
I would like to use r.grow.distance on a raster map. I rasterised a shapefile (about 50 cities in a region) and we can see them on the raster map.

Yet, when I launch r.grow.distance, I have a distance raster map from the entire original raster, and not from the features in my raster map (the cities). I have a big dimension raster map with completely useless distances.

I tried many solutions (resolution, rasterization,...) and nothing worked...

If anyone has a suggestion, it would be great.


texts: Ignore Fields in Target Table When using FieldMapping in Append Tool

texts:
I am transforming one feature class into another (currently empty) feature class.

The input feature class and output feature class share the same fields, but the output feature class has some additional fields. I'm aware I can do this using a Data Acess module InsertCursor, but I'm trying to learn how to use the FieldMappings module and got stuck here.

 staging = 'C:/Path/To/Input.gdb'
 gdb = 'C:/Path/To/Output.gdb'
 input_fc = os.path.join(staging, 'Hazard')
 output_fc = os.path.join(gdb, 'Hazard')

# Get a list of common fields
common_fields = [f for f in input_fields if f.name in [g.name for g in output_fields] and not f.required]
# Get a list of fields found in the output, but not in the input
uncommon_fields = [f for f in output_fields if not f.name in [g.name for g in input_fields] and not f.required]
# Get a list of fields found in the input, but not in the output
unmapped_fields = [f for f in input_fields if f.name not in [g.name for g in output_fields] and not f.required]
print [f.name for f in common_fields]
print '*' * 20
print [f.name for f in uncommon_fields]
print '*' * 20
print [f.name for f in unmapped_fields]

This returns

[u'hazard_type', u'inspection_staff', u'comments', u'inspection_date', u'pp_survey_cycle', u'pp_timestamp', u'record_uid', u'record_timestamp', u'flagged', u'reason_flagged', u'marked_for_deletion', u'data_quality_assured', u'data_qa_staff', u'data_qa_comments', u'data_qa_timestamp', u'gdb_timestamp']
********************
[u'pp_azu', u'pp_city', u'pp_state', u'pp_zip', u'pp_survey_group']
********************
[]

What I want to happen is for fields in common_fields to get mapped over to output_fc and fields in uncommon_fields to remain Null.

Here's the code I'm using, which runs without error but doesn't append anything (grr ESRI...). Do I have to specifically remove fields from the FieldMappings object? How do I do that?

# Specify the mappings using a dictionary. We're just mapping fields that have the same name to one another here.
custom_mappings = {}

for f in input_fields:
    if not f.required:
        custom_mappings[f] = next((g for g in output_fields if g.name == f.name), None)

# Create field mapping object
fms = arcpy.FieldMappings()

# Specify the field map for each field and add it to the FieldMappings object
for field in input_fields:
    if not field.required:
        fm = arcpy.FieldMap()
        fm.addInputField(input_fc, field.name)
        fm.outputField = custom_mappings[field]
        fms.addFieldMap(fm)

# Append to the output_fc using the field mappings
arcpy.Append_management(input_fc, output_fc, 'NO_TEST', fms)


texts: How to map places with zip codes that have leading zeros using CartoDB?

texts: I have a large set of data that has several places on the East coast that have leading Zeros in the Zip.  How do I map these when leading zeros are automatically removed?


texts: Cannot Edit Python Script tool using its context menu from ArcMap/ArcCatalog?

texts:
I am unable to edit scripts in ArcMap or ArcCatalog.

I navigate to the script, right-click on it and select "Edit", then nothing happens. I am using notepad++ as my editor (set in geoprocessing options). I have looked into this issue (http://support.esri.com/fr/knowledgebase/techarticles/detail/38581) and have attempted the steps therein. However I am still getting the same result, nothing happens.

I am using ArcGIS for Desktop Basic Version 10.3.1 and Python 2.6

What steps should I attempt to solve this issue?


texts: Is it possible to create animated markers in ArcMap?

texts:
I remember reading somewhere last year that arcgis 10 would allow the creation of animated markers (alike animated .gifs) that could be either exported to pdfs or seen inside arcmap.

Does such feature actually exist?


texts: How to upload data from QGIS to a SQL Server?

texts: I'm new at QGIS. I'm wondering if there's a native way or plugin to allow me to upload a shapefile into my SQL Server directly from QGIS. So far the only solution that I could find is the Shape2SQL tool.


texts: Merging TIFFs while retaining details in GDAL

texts:
I am currently trying to merge overlapping raster TIFFs into one single, continuous image. Just using Merge results in a raster with borders, and it is unsatisfactory :

However, using a pixel function to choose the maximum value present on overlapping parts results in a smooth image, but much of the detail is lost.
How do I merge them such that they look more or less smooth but maintain the details? I can use any Python or GDAL solution


texts: Area calculations don't agree after raster to polygon conversion

texts:
I have a landsat 8 image [UTM 19s] (30m resolution, so 900m sq per pixel).
I need to calculate the total area (hectares) from a land cover type.
I take the count of pixels from the attribute table, calculate field using the ([count]900).0001 to get total hectares.

As a check, I convert the raster to a polygon and do not simplify so that it returns polygons in a stepped fashion and shaped just like the pixels.This allows me to use the calculate geometry tool to populate an area field. This tool returns different values than the straight calculation from the raster itself. It is off by a few hundred hectares which is far too much for my reporting purposes.

How can I "convince" the two calculations to agree? is one method more precise than the other? I would be inclined to go with the straight conversion from the raster itself rather than relying on the arc tool to produce the areas...Thoughts or suggestions?


texts:
Is there an introduction or guidelines how to read the QGIS API documentation?

For example, I want to write a vector layer. A quick search brings me to the API documentation page of QgsVectorFileWriter.

This is a very long page with an inheritence diagram, different classes, public classes, public membership functions,... but I don't fully understand what this means.

I can figure out that I need the QgsVectorFileWriter::QgsVectorFileWriter function, but also in this part I have some questions:

here are three identical(?) functions listed?
which variables are mandatory and which are optional?

I've found some useful resources on pyqgis at:

Learning resources for PyQGIS?
Understanding the QGIS API documentation to write python code

but afaik these resources don't give a gentle introduction on how to read the API docs.

I think understanding how to read the API documentation would really help me to write my own code rather than searching for examples.


texts:
I have just one table of street segments of Berlin. I have a query that outputs all segments of category > 1, that are within a distance of 2000m of those segments which are of category 1.

This works fine:

SELECT DISTINCT s.id, s.geom
    FROM table s
    INNER JOIN table m
    ON ST_DWithin(s.geom, m.geom, 2000)
    WHERE m.category = 1
    AND s.category > 1;

But now i want to exclude from the result set those segments of s that intersect with m. So i just appended the following expression:

...
AND NOT ST_Intersects(s.geom, m.geom);

With this additional line the result set is the same as without, though i know there are segments that meet this condition. If i omit the 'NOT', asking for segments that DO intersect, i get about 400 segments.
I also tried to use Disjoint, with the same problem.

Any idea?


texts: Google Maps directions link doesn't work properly in the Explorer for ArcGIS mobile app?

texts:
I have an AGOL web map for finding beach accesses in Florida. Within my pop-up I have a Google Maps link to get directions to the current access. Everything works fine on desktop, but the directions link doesn't work properly in the Explorer app. When I click the directions link it opens Google Maps with the destination as the Gulf of Mexico and the Current Location is way off (for example I am in Tallahassee,FL and it said I was at a Redbox in Palm Harbor,FL about 330 miles away). When looking at the same accesses on my desktop the directions work perfectly.

Directions Link in popup:
https://www.google.com/maps/dir/Current+Location/{Y_LATITUDE},{X_LONGITUD}

I have also tried:
http://maps.google.com/maps?saddr=Current+Location&daddr={Y_LATITUDE},{X_LONGITUD}

Some context:
I am using Explorer app on IOS 8.3 and the web map is out of the box no custom code.The Lat/Long is populated from fields in the attributes.My location settings allow for Safari, Google Chrome and the explorer app to have access to my location and all three accurately show my location when I use their location features outside of this instance. This error occurs whether I am on WiFi or 4G.

When in the Explorer app if I copy the link instead of clicking, this is what is copied to my clipboard(looks ok except the lat/long values are truncated):
http://maps.google.com/maps?saddr=Current+Location&daddr=29.95,-85.42

This is the URL that is shown when the link is clicked and Google Maps opens in Safari:
https://www.google.com/maps/dir/Redbox,+207495-3414+Tampa+Road,+Palm+Harbor,+FL+34684/29.95,-85.42/@29.3194888,-83.8341354,9z/data=!4m6!4m5!1m2!1m1!1s0x88c2ed70bc996c9f:0xbf08cf8cc97df53f!1m0!3e0

When compared to the same access on the desktop:

Copied to clipboard:
https://www.google.com/maps/dir/Current+Location/29.9460010011,-85.4170700002

When clicked and GM is opened in the browser:
https://www.google.com/maps/dir/30.4305062,-84.2542195/29.946001,-85.41707/@30.2169314,-85.3848336,9z/data=!3m1!4b1!4m4!4m3!1m1!4e1!1m0

I apologize for the long explanation, if there is feedback in comments I will edit as necessary to clarify. I was unsure what info is useful so I added what I thought may be. I can also add screen shots if needed. I will also try adding the URL with a blank start location as this may be better from a UX standpoint.


texts: How to emulate a specific vector grid in QGIS?

texts:
I would like to analyze data included in a pdf file. The file contains information on characteristics of 16km grid cells for a specific study area. The file also includes a map displaying the grid. I know the coordinate system and datum that has been used to create this map.

I would like to create the exact same grid in QGIS as to merge the data from the pdf file with further geographical information on the study area. I am a total QGIS noob and I don't know if what I want to do is possible at all.

What I have done so far is: I copied, digitized and and georeferenced the map contained in the pdf. I then set the CRS according to the information I got from the pdf. I then tried to "reproduce" the grid with 16km cells using the QGIS vector grid tool.

Now, here's my problem: the side length I enter into the vector grid tool is in degrees not km. I have measured the length of the grids in degrees (using the QGis distance measurement tool). If I try to create the vector grid using the degree measure for the side length I get cells that are approx. the size of the original map at some parts of the map but not at others. Finally, even if I got a grid of the same size as the original, I am not quite sure how to adjust the location of the grid I created to the exact same location as the original one.


texts: esri rest service ?f=sitemap returns no services

texts:
I have some code that queries the esri webserver for a list of services.

I start off with:

https://MyEsriUrl/arcgis/rest/services?f=sitemap

That code worked for years.

We recently attached a portal to our esri webserver.  We're now running v10.5 of both.

Now, the esri admin can issue this command in the browser and it still works, but if I issue it I get an empty list back.

I thought it might be a permissions issue except for three facts:

(1) If I go to https://MyEsriUrl/arcgis/rest/services I can see everything I ought to see.

(2) I've found out that I can now use https://MyEsriUrl/arcgis/rest/services?f=json and that also works correctly.

(3) Esri admin added me to the esri admin role and I still can't see anything with ?f=sitemap.

Obviously, I can rewrite the code to work, but that app is scheduled to be replaced and I would rather spend time working on the replacement app...

Any ideas on how to get ?f=sitemap to work properly?


texts: Removing cloud cover using Sentinel Toolbox?

texts:
According to the ESA website it would appear possible to remove cloud cover form Sentinel data using the Seninel-2-Toolbox.

Is this really possible or are they refering to the method of mosaicking images from different dates?


texts: QGIS - line between ordered layer of points - take attributes of both points

texts:
I have a series of locations in time in a point layer. I want to make a line layer of lines between consecutive points that takes the attributes (in my case, a timestamp) of the 1 points it connects.
Is there an existing plugin or script that does this? The plugins I found makes it a long line.

A solution I thought of was to use PointtoPaths plugin, then explode the resulting line, and use a spatial join to somehow get the attributes of the points.


texts: Converting census data to block groups with regard to their area

texts: I have population data for census blocks and want to convert them to block groups. I know that block groups has an extra digits, but how it will be possible to break down census data between blocks weighted by their area!I have little experience with qgis and have to use it for this purpose.


texts: Tearing in ggmap (R) after using spTransform and random point generation using splancs

texts:
I am using R and the package ggmap to try and generate some random points within the city limits of Chapel Hill, NC. I have downloaded the shapefile for the town of Chapel Hill from the Orange County website: http://web.co.orange.nc.us/gisdownloads. The city limits for Chapel Hill and several other nearby cities are in the file labeled city.zip. I saved this file to my hard drive and unzipped it. My code is as follows:

library(ggmap)
library(rgdal)

map1 <- get_map(location='Chapel Hill, NC', zoom=12, maptype='roadmap')

setwd('<Path to unzipped city.zip files>')
city <- readOGR('city.shp','city')
sub.city <- city[city$CITYCODE == "CH",] # Take only the part of the shapefile we need (the area around Chapel Hill)
city1 <- spTransform(sub.city, CRS('+proj=longlat +ellps=WGS84'))

If I simply plot city1 the polygon looks as it should:

plot(city1)

but if I try to plot it over my street map of Chapel Hill, I get some weird tearing in the northwest part of the city:

ggmap(map1)+
  geom_polygon(data=city1, aes(x=long, y=lat), color='blue', fill='blue', alpha=.3)

Why is this and how could I fix it?

The second thing I want to do is generate some random points within this polygon. The package splancs claims to do this with its csr() function, but so far I cannot get this to work.

library(splancs)
poly1pl <- coordinates(city1) # Extracts long and lat into a matrix
set.seed(2913)
randpoints <- data.frame(csr(poly1pl, 20))

The last command of this gives me an error: Error in runif(ngen) : invalid arguments
I get the same error if I try using csr() on sub.city as well, so the error is not a result of using spTransform. What does this error mean and how can I get this to work? Or is there an easier way to generate these points?


texts: Does anyone have a categorized Geographic coordinate system List?

texts:
We are making a custom desktop application which will project a selected GDB or shapefile to a selected coordinate system. This will be a custom app since most of the users will not have access to GIS or are not even GIS people. The projection will run on a server, which will have the tools.
Problem is, we need the user to specify the projection. Although the list of projections are available (https://developers.arcgis.com/rest/services-reference/projected-coordinate-systems.htm), it is not categorized... Imagine having a dropdown with 5000 items... So we are thinking 3 levels, like in ArcGIS... But I can't find those categories anywhere. I was hoping an excel file with 4 columns - The WKID, WK Text, Parent, System. This could work well I believe. Unless there is a library somewhere available?


texts: Reduce space in between subplots containing choropleth maps

texts:
(I posted this on Stackoverflow, but thought this may be a better venue)

I am trying to generate a figure containing six choropleth maps in R. Currently my code looks something like this:

library(rgdal)
library(GISTools)

# code that loads data needed to generate choropleth maps not shown

par(mfrow = c(3,2), mar=c(0,0,0,0), oma=c(0,6,6,0))

choropleth(mn_county, mn_county$emppct.tert.2006, shades)
mtext( 'Secondary', side=3, line=2)
mtext( '1998', side=2, line=2)

choropleth(mn_county, mn_county$emppct.tert.2006, shades)
# choropleth maps are repeated in this post for simplicity
mtext( 'Tertiary', side=3, line=2)

choropleth(mn_county, mn_county$emppct.tert.2006, shades)
mtext( '2006', side=2, line=2)

choropleth(mn_county, mn_county$emppct.tert.2006, shades)
choropleth(mn_county, mn_county$emppct.tert.2006, shades)
mtext( '2014', side=2, line=2)

choropleth(mn_county, mn_county$emppct.tert.2006, shades)

The problem is, even though I set margins mar to be zero throughout, there are still big gaps between the columns (see below). Not just between the choropleth maps, but also between the row headings (1998, 2006, 2014) and the maps.

How can I reduce this wide space?

P.S. I've been trying to find a GIS shapefile which I can directly pull off from the internet (I don't think I can attach shapefiles in posts), but can't find one. Any idea how I can circumvent this problem?


texts:
I have a CSV text file which i got from GBIF containing a species distribution.
How would I arrange it so that I could display it in DIVA GIS?


texts: Spatialite 'SELECT UpdateLayerStatistics()' returns "0"

texts:
I have got several sqlite-databases (some geneated from within qgis 2.14.0/2.16.2). When running SELECT UpdateLayerStatistics() to the databases some results are "0" and some are "1".
A PRAGMA integrity_check; to the databases with "0"-result haves the result "OK".
Running the SELECT UpdateLayerStatistics('myTableNames') against all tables within the "0"-result-db have the result "1".

So why is there a "0"-result at several databases when running SELECT UpdateLayerStatistics() against the whole db - althrough running the command against every single table returns "1". What can I do to fix the problems and finally get a "1"-result ?


texts:
I am attempting to teach myself how to make a basic GetFeatureInfo request on a WMS (specifically from Geomet, GetCapabilities here). The following returns 400 Bad Request, and I'm not sure if it is a simple syntax error, or something more systematic.

Request:

http://geo.weather.gc.ca/geomet/?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&LAYERS=GDPS.ETA_GZ&STYLES=&SRS=EPSG%3A3857&BBOX=-20037508.342789244%2C5009377.085697313%2C-17532819.799940586%2C7514065.628545966&WIDTH=256&HEIGHT=256&QUERY_LAYERS=GDPS.ETA_GZ&X=100&Y=100

Parameters:

SERVICE:WMS
VERSION:1.1.1
REQUEST:GetFeatureInfo
LAYERS:GDPS.ETA_GZ
STYLES:
SRS:EPSG:3857
BBOX:-20037508.342789244,5009377.085697313,-17532819.799940586,7514065.628545966
WIDTH:256
HEIGHT:256
QUERY_LAYERS:GDPS.ETA_GZ
X:100
Y:100

On the other hand, the following similar GetMap request works:

Request:

http://geo.weather.gc.ca/geomet/?SERVICE=WMS&REQUEST=GetMap&VERSION=1.1.1&LAYERS=GDPS.ETA_GZ&STYLES=&FORMAT=image%2Fpng&TRANSPARENT=true&HEIGHT=256&WIDTH=256&SRS=EPSG%3A3857&BBOX=-17532819.799940586,5009377.085697313,-15028131.257091932,7514065.628545966

Parameters:

SERVICE:WMS
REQUEST:GetMap
VERSION:1.1.1
LAYERS:GDPS.ETA_GZ
STYLES:
FORMAT:image/png
TRANSPARENT:true
HEIGHT:256
WIDTH:256
SRS:EPSG:3857
BBOX:-17532819.799940586,5009377.085697313,-15028131.257091932,7514065.628545966

The GetCapabilities document lists the GDPS.ETA_GZ layer's SRS as 4326, but I've attempted using lat/long coordinates, as well, to no avail.


texts: ECW files no longer display or open in ArcCatalog or ArcMap?

texts:
When I closed and reopened ArcCatalog yesterday ECW files will no longer display in the preview window and will not open in ArcMap.  I am using ArcGIS Desktop 10.4.

I have checked supported raster formats and they are still checked.  Have done a repair installation to no effect.

Has anyone else struck this issue?


texts:
In GeoServer, I published a polygon shapefile with default polygon SLD, while in preview when I click on feature it shows the error rendering process failed null.

Please explain what is this error, what I have yo do to overcome this error?


texts: "Points To Raster Tool"- Output pixel values differ from original point values

texts:
I am trying to convert an evenly spaced set of points to a raster using the "Points To Raster" tool in 10.1.
Problem: Generally, the raster makes sense, but in some areas the cell values differ quite significantly from the original point values (e.g. point= 0, raster cell= 87).
I have already tried different cell sizes and assignment types (although I would like the MEAN).
Would anyone know why that is and what I need to do in order to obtain "meaningful" values (cell values that represent the original point values).


texts:
When I used QGistoThreejs plugin in QGIS and run after adding a tif file... the output generated with the files created in the specified location....

Let's say if I name it as sample.html when I run the plugin... the following files will be generated with the specified location

dat-gui(folder)
threejs
dat-gui_panel.js
Qgis2threejs (css file)
Qgis2threejs.js
sample.html
sample.js

But if I want to change the DEM dynamically to display in HTML page that was created, I observed that we need to change sample.js file... (i.e, javascript file with specified name) should be overridden for the DEM selected.

So, is there any way to override that javascript file using web?


texts: Adjusting condition in field calculator in QGIS

texts:
I need to calculate a score based on an area field using these ranges:

if >= 30000 (square meters) then 4 (score)
if < 30000 AND >= 10000 then 3
if < 10000 AND >= 5000 then 1
if < 5000 then 0

I've tried with the CASE function but it only works for lines with a single operator (i.e. >=30000 then...) Lines with two operators (where I wrote AND) are seemingly not supported. Here is one of the syntaxes I've tried.

CASE
WHEN Surface >= 30000 THEN 4
WHEN Surface < 30000 >= 10000 THEN 3
WHEN Surface < 10000 >= 5000 THEN 1
ELSE 0
END

The result of this exact command in the score column are only 4 or 0. Which makes me understand that the first line works but not the others.

I know SQL has its BETWEEN operator, which would probably be a good solution to my syntax issue, but I don't see it in the field calculator.

So I'm wondering if another function could help me here or if I just need the correct syntax.


texts: How to expand a polygon to a minimum allowable distance?

texts:
I am looking for a method to expand certain areas of a single polygon to meet a minimum distance requirement. For example, the screenshot shows a polygon in green where I want to expand areas to a minimum distance of 2.5 km. The solid black polygon is a rough example of what I hope to accomplish. Buffering everything by X distance will not work because I am trying to preserve as much of the original shape and area as possible. How can I expand certain areas to meet the minimum 2.5 km requirement and minimize any increase in total area?

I am interested in any vector or raster based approaches and I am using ArcGIS for Desktop with a Basic license and the ArcGIS distribution of Python.


texts:
I've prepared my Shapefile in ArcMap 10.7.1 by Raster to Polygon. When I am opening this shapefile in ArcMap its showing a rectangular polygon around the shapefile like this:

When I select this rectangular shapefile from Open attribute table, its select all other part of shapefile. I just want to remove the boundary level shapefile, so that, I able to calculate area of shapefile.
Output looks like this