texts: PostGIS 2.3 Splitting multiline by points

texts:
I am struggling to make work the ST_Split function, since it seems to be the way to go. I have this situation, and I want to end with the line split in three on according to the points:

And here is my code, the cable2's geometry is MultiLineString and the boite :

WITH brs AS (SELECT row_number() over(), boite.geom FROM boite, cable2
    WHERE st_intersects(boite.geom, cable2.geom)
    ORDER BY st_line_locate_point(st_linemerge(cable2.geom),boite.geom)),
    cables AS (SELECT (st_dump(st_split(cable2.geom, brs.geom))).geom FROM brs, cable2)
SELECT st_astext(geom) FROM cables WHERE st_endpoint(st_linemerge(geom)) = (SELECT geom FROM brs WHERE brs.row_number = 4)

Here I am just trying to get the last geometry, but I get this as a result, where the first and last line are the same, and anyway the endpoint is always the same :

"LINESTRING(949758.692587021 6567419.6350433,949790.867103036 6567425.81704638,949825.149120101 6567402.63453484,949835.265125137 6567380.857024)"
"LINESTRING(949825.149120101 6567402.63453484,949835.265125137 6567380.857024)"
"LINESTRING(949790.867103036 6567425.81704638,949825.149120101 6567402.63453484,949835.265125137 6567380.857024)"
"LINESTRING(949758.692587021 6567419.6350433,949790.867103036 6567425.81704638,949825.149120101 6567402.63453484,949835.265125137 6567380.857024)"

Would you happen to know what I am doing wrong? I tried to look in other similar requests here but unsuccessfully.

**EDIT ** Using the solution proposed by Redoute, and since I am trying to implement that in a trigger, I fail to do so, and here is my code :

CREATE OR REPLACE FUNCTION public.split_cable()
  RETURNS trigger AS
$BODY$
DECLARE compte integer;
DECLARE i integer := 2;
BEGIN
compte = (SELECT count(*) FROM boite WHERE st_intersects(boite.geom, new.geom));

WHILE i < compte LOOP
    WITH brs AS (SELECT row_number() over(), boite.geom FROM boite, cable2
        WHERE st_intersects(boite.geom, new.geom)
-- here the ORDER BY serve to get the "boite" objects in a specific order
        ORDER BY st_linelocatepoint(st_linemerge(new.geom),boite.geom)),
        brs2 AS (SELECT st_union(geom) AS geom FROM brs),
        cables AS (SELECT (st_dump(st_split(new.geom, brs2.geom))).geom FROM brs2)
    INSERT INTO cable2 (geom) VALUES (
    SELECT st_multi(cables.geom) FROM cables WHERE st_startpoint(geom) = (SELECT geom FROM brs WHERE brs.row_number = i));
    i = i + 1;
END LOOP;

new.geom = (WITH brs AS (SELECT row_number() over(), boite.geom FROM boite, cable2
        WHERE st_intersects(boite.geom, new.geom)
        ORDER BY st_linelocatepoint(st_linemerge(new.geom),boite.geom)),
        brs2 AS (SELECT st_union(geom) as geom from brs),
        cables AS (SELECT (st_dump(st_split(new.geom, brs2.geom))).geom FROM brs2)
        SELECT st_multi(cables.geom) FROM cables WHERE st_startpoint(geom) = (SELECT geom FROM brs WHERE brs.row_number = 1));
RETURN new;
END
$BODY$
  LANGUAGE plpgsql;

EDIT 2 As suggested by Redoute, I asked the question here : https://stackoverflow.com/questions/46444147/postgis-2-3-split-a-line-by-points


texts: I'm running into an issue where my MBtiles files are exceeding the max file size for the file system when created with gdal. Is there a known method of creating a set of mbtiles files that don't exceed a set file size? or will I need to create something?


texts: How to add Georeferenced PDF as layer to QGIS 2.0?

texts: I have some georeferenced pdf's that I would like to open in QGIS. I have read that GDAL version 1.8 and higher will support pdf, and I have GDAL 1.10 installed with QGIS 2.0 on Windows 7, but when I add raster layer then choose my pdf file, I get an error that my file is not a supported source.


texts: ArcGIS/ ArcPy - XY coordinates of feature class does not match calculated geometry XY

texts:
Tested on: different workstations, different Windows versions; ArcGIS Desktop 10.3, different license levels

Background:

I use the Google Maps DistanceMatrix to calculate drive times between two point feature classes (FC). Therefore I use the XY coordinates with simplejson (instead of street and city names). I get this XY coordinates from the attribute table of a point FC. This part of the workflow is working well.

But before I can do that, I have to geocode adresses from an Excel file to get the needed FC. I do this with the Geopy module. In a part of my script I use this also to calculate the X and Y coordinates of the feature and put it in the attribute table of the FC, because it´s needed for the calculation of drive times. It´s also working great. My results are looking like 47.5925516, 7.6168828, my FC is in WGS84 (4326).

Problem:

Most of the features are geocoded correctly, but some points - because they are not so exactly how I need them - need to be moved by hand.
My idea for that was easy: If I have to move a feature, my XY coordinates also need to be updated. Simple click on calculate geometry (attribute table), set it to decimal degrees and finished. But if I do that without moving a feature, the calculation result - not for all features but a lagre part of them - differs from my original XY coordinates, which I get from Geopy/ GoogleMaps in WGS84. For example the original XY coordinate 47.5925516, 7.6168828 were calculated to 47.5923835, 7.61109560000006. Only the first few decimal places are correct.

Because I need the XY coordinates for my drive time matrix and I am working with very short distances, the result could be a significant deviation. If I use a aerial image or a basemap to check the points, the position of the original XY coordinates/ the features are exact but the calculated position could be more then 500 m away.

What I have tried:

What I am doing wrong, this could not be an typical projection inaccuracy,or?! I tried different fieldtypes (double, text) and precisions. I changed my layer coordinate system to WGS84, WGS84 auxiliary sphere and DHDN_3_Degree_Gauss_Zone_3. Nothing worked for me.

*******UPDATE********

I geocoded the adresses with the "original" Google Maps API (same FC, in WGS 84, too) and it works. If I put both FC (Geopy and original Google Maps API) in the same Layer and calculating the XY coordinates, the old ones are still calculated wrong, the Google API ones are calculated right.

They are both in the same coordinate system (WGS84)
visually the point FCs have the same location on the map


texts:
I have a spatial postgis database named "testdata" having two schema named "afgdata" and "pakdata" with common table names viz. building_a, places, railways, roads, water_a, waterways of different area.

I am not familiar with both Postgis and FME.

I want to merge both schema in to new schema  named "mergedata" within same database. Can anybody help me out to do it with FME or direct in PostgreSQL?

Is it possible to merge both schema and output in new database in public schema?


texts: Rendering ol.style.Text as HTML in OpenLayers 3?

texts:
I am trying to display some text in my feature using the vector style settings.
I am doing as follow:

new ol.style.Style({
    text: new ol.style.Text({
      offsetY: -10,
      font: '18px Calibri,sans-serif',
      text: feature.get('label'),
      fill: new ol.style.Fill({
        color: '#ec971f'
      }),
      stroke: new ol.style.Stroke({
        color: '#fff',
        width: 3
      })
    })
});

in order to display a label on my map. The problem is that the html is rendered as text.
for example: 10km2 is rendered as 10km2

Any idea how to do this?


texts:
I am trying to mosaic together two smaller mosaics with two Landsat scenes each in ArcMap 10.2 but have come across an issue. As background, I am mosaicking the scenes this way since each of the smaller mosaics' bands were stretched initially and mosaicking all four scenes before stretching produces a poorly displayed raster. Stretching the bands on that initial mosaic did not fix the issue. When attempting to mosaic these smaller ones together, the following larger mosaic is produced:

Overlaps do exist between all four scenes but when combining all of them initially into a single mosaic, the collar was not in between the scenes like it is now. Removing the collar from the larger mosaic does not remove the collar seen in between the two smaller mosaics:

How can I go about making a cohesive mosaic of these two smaller mosaics? Is there a setting or tool within Raster toolbox that may be of use?


texts: PyQGIS - QgsSpatialIndex.intersects returns unexpected results

texts:
I have a polygon and a line. The line starts very close to the polygon but the end vertex is further away. I create a rectangle around start and end vertex. The spatial index of the polygon is used for intersecting both rectangle. I expect only one rectangle to intersect but I get intersection for both. Why?
The screenshot shows the extent of the rectangles - only one intersects with the polygon.

QGIS 2.18, Win7.

Data source: SHP files.

# -*- coding: utf-8 -*-

# is start/endpoint of section close to a buildung within tolerance?
tolarance_in_meters = 0.5

# Get the active layer
layer1 = QgsMapLayerRegistry.instance().mapLayersByName("line")[0]
layer1.setProviderEncoding(u'UTF-8')
layer2 = QgsMapLayerRegistry.instance().mapLayersByName("polygon")[0]
layer2.setProviderEncoding(u'UTF-8')

# create indices
index1 = QgsSpatialIndex(layer1.getFeatures())
index2 = QgsSpatialIndex(layer2.getFeatures())

# create a memory layer with two points
polyTemp =  QgsVectorLayer('Polygon', 'temp' , "memory")
polyTemp.setProviderEncoding(u'UTF-8')
pr = polyTemp.dataProvider()

def getRectangle(point, size):
    x = point.x()
    y = point.y()
    ll_x = x-size
    ll_y = y-size
    ur_x = x+size
    ur_y = y+size
    return QgsRectangle(ll_x,ll_y,ur_x,ur_y)

iter = layer1.getFeatures()
for feature in iter:
    # feature is a line (not a multi - line)
    geom = feature.geometry()
    # we need first and last vertex
    vertices = geom.asPolyline()
    qgis_point_vFirst = geom.vertexAt(0)
    qgis_point_vLast = geom.vertexAt(num_vertices-1)
    rect1 = getRectangle(qgis_point_vFirst, tolarance_in_meters)
    intersection_vFirst_ids =  index2.intersects(rect1)
    rect2 = getRectangle(qgis_point_vLast, tolarance_in_meters)
    intersection_vLast_ids =  index2.intersects(rect2)

    # add the rectangles
    pt = QgsFeature()
    r =  QgsGeometry().fromRect(rect1)
    pt.setGeometry(r)
    pr.addFeatures([pt])
    pt = QgsFeature()
    r =  QgsGeometry().fromRect(rect2)
    pt.setGeometry(r)
    pr.addFeatures([pt])

   numInter1 = len(intersection_vFirst_ids)
   numInter2 = len(intersection_vLast_ids)

   print("num intersections: {0}, {1}".format(len(intersection_vFirst_ids), len(intersection_vLast_ids)))

polyTemp.updateExtents()
polyTemp.commitChanges()
QgsMapLayerRegistry.instance().addMapLayers([polyTemp])


texts: Reclassification and Con(0 not working arcpy

texts:
I have a raster that has three pixel values=0,1,2. I want all pixels of value 2 to be 1. The raster is in tif format, 32 bit signed type. Below is my arcpy script-

inRaster1="path address/input.tif"

arcpy.CheckOutExtension("Spatial")

#Reclassifying
outReclass=Reclassify(inRaster1, "Value", RemapValue([[0,0],[1,1],[2,1]]))
print "Reclassification is complete"

outReclass.save("address path/output.tif")
print "Saved. Check"

arcpy.CheckInExtension("Spatial")

However, I get the 9999 Error and I do not know how to fix it! So instead, I used the Con tool as below and saved the output.

outCon = Con(inRaster1==2, 1, inRaster1)
print "Reclassification is complete"

BUt the pixel values do no change i.e. the saved raster still has pixel values=0,1,2 instead of 0 & 1. And this was completed with no errors.
I then modified the code line to

outCon=Con(inRaster1>1,1,inRaster1)

And this time all values were reclassified to 1 i.e. even the pixels with value 0.

What is going wrong with the Reclassify() and why is Con() not reclassifying 2.How do i reclassify pixel value=2 to 1


texts: Changing pixel's dimension (raster) on QGIS

texts: I would like to change the size of the pixels from 2.5 to 10. How can I achieve this?


texts: Calculating block statistics from land cover raster using R?

texts:
I have a raster (5 x 5 m resolution) representing land cover. Each cell contains a integer value (from 1 to 10) representing a total of ten land cover classes (e.g. 1-forest; 2-agricultural areas; 3-urban areas;(...)10-scrublands).

From this land cover raster, I would like to generate 10 rasters with a resolution of 10 x 10 m representing the proportion of each land cover class per cell. In another words, I would like to produce a raster representing the proportion of forest per cell, another raster representing the proportion of agricultural areas per cell, etc... (please see the figure below that shows what I want to accomplish).

I know this task can be accomplished in ArcGIS Spatial Analyst by generating 10 binary rasters (1- presence of a specific land use class; 0-absence of a specific land use class) and then using the "Block statistics" option. However, because I do not hold a license for ArcGIS, this is not a valid option.

Instead, I would like to know how to do this in R (version 3.40). I have checked previous posts and I managed to find one where an user intended to do something similar (Calculating Proportion of Land Cover Classes with moving window around point in R?). However, I do not want to calculate proportion of land cover with a moving window.

Below find the code to generate a dummy RasterLayer object with five land cover classes coded as integers.

    r<-raster(ncols=100, nrows=100,xmn=100,xmx=1100,ymn=0,ymx=1000, crs="+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

    x<-0:5

    r[]<-sample(x,10000,replace=T)


texts: How to combine (join?) two overlapping polygons to determine combined area they cover?

texts:
How do I combine (join?) two overlapping polygons to determine the combined area they cover in ArcGIS 10?

I am sure the Geoprocessing tools can do it, but just how eludes me.


texts:
I have a MapServer map file that gets data from PostGIS tables. One of the tables holds several linestrings that I set as points in the map file. They show up based on the zoom-level in my map.

Is there any way that I can tell the map file to display all points from the table at the same time regardless of zoom-level?


texts: QGIS cannot save to SpatiaLite after vector difference?

texts:
In QGIS I have a spatialite DB full of roads. I'd like to do a Difference operation (Vector > Geometry Tools > Difference) with another layer in order to keep roads only available in a given zone. The difference works fine but once I try to save that new layer to a Spatialite DB I get OGR error saying that it can't insert GEOMETRYCOLLECTION where a MULTILINESTRING is expected.

I don't know how to fix that and what to do.

Following the @she_weeds' answer it was suggested to try saving that layer to another format and then to Spatialite. Saving it to GeoJSON worked, I could convert it to Spatialite after that and import it to my QGIS project (even if the auto-import didn't worked at first).

Even if everything works now using this workaround, I will update QGIS to latest version too.


texts: Why are there null geometry lines after using the contour tool in Spatial Analyst?

texts:
A few months ago, I ran the Contour tool on two separate DEMs (one was smoothed and the other wasn't).  The contour interval was set to 1m.

Today, I noticed that there were several features created with NULL geometry.

I created a field called partCount and populated it with the following python expression:

!SHAPE!.partCount

and then I applied the following query to my dataset:

partCount =0

You can see, there are several features with no parts and zero length.

However, they did get an elevation value.

My question is: why would the Contour tool create features with null geometry?

Note: the DEMs did contain No Data areas, but I'm not yet convinced this is the cause of the null geometry.

I created the contours in ArcGIS 10.3.1.


texts: How to export data from saga or tas into excel?

texts:
I want to derive terrain attributes from DEM with saga and whiteBox GAT (SAT). I have problem  in exporting data into excel.
would you help me?


texts: viewshed analysis for observers wearing visors?

texts:
I'm starting with a DEM, and looking to identify every point that can be directly "seen" from any point in a particular polygon on an overlaid layer defining an area on that DEM - like any other viewshed analysis, but with the additional stipulation that the observer points can only "see" at certain downward angles.

To help motivate what I'm looking for, my real goal is to identify run-out zones for avalanches.  It's easy to identify "slideable slopes": just everything that's, say, 30-45 degrees.  But if I want to identify the areas which snow would slide from the slideable areas into, the above quasi-viewshed starts to look like what I'm after.

My exploring the surface and hydrology tools (and help files for the tools) in Arc hasn't turned up anything that I think will do the trick.  Any ideas - either in Arc or otherwise?

Thanks in advance for any insight.


texts:
I am relatively new to python and I am trying to build a tool within python toolbox with GPSQLExpression as parameter.

I am having trouble with populating its field and value list. Below is part of my code, I put it in updateParameters.

Feature classes in geodatabase were made from OSM data.

def getParameterInfo(self):
    """Define parameter definitions"""

    p1 = arcpy.Parameter(
        displayName = "In gdb",
        name = "in_gdb",
        datatype = "DEWorkspace",
        parameterType = "Required",
        direction = "Input")

    p2 = arcpy.Parameter(
        displayName = "Geometry type",
        name = "userinput",
        datatype = "GPString",
        parameterType = "Required",
        direction = "Input",
        multiValue = False)

    p2.filter.type = "ValueList"
    p2.filter.list = ['pt','ln','ply']

    p3 = arcpy.Parameter(
        displayName ="Selected fc",
        name = "fcs_of_gdb",
        datatype = "GPFeatureLayer",
        parameterType = "Required",
        direction = "Input",
        multiValue = True)

    p4 = arcpy.Parameter(
        displayName = "Pick counties",
        name = "colli",
        datatype = "GPString",
        parameterType = "Required",
        direction = "Input",
        multiValue = True)

    p4.filter.type = "ValueList"

    p5 = arcpy.Parameter(
        displayName = "Select Feature Layer",
        name = "fl_gdb",
        datatype = "GPFeatureLayer",
        parameterType = "Required",
        direction = "Input")

    p6 = arcpy.Parameter(
        displayName = "SQL Expresssion",
        name = "expres",
        datatype = "GPSQLExpression",
        parameterType = "Required",
        direction = "Input",
        multiValue = False)

    p6.parameterDependencies = [p5.name]

    p = [
        p1, p2, p3,
        p4, p5, p6
    ]
    return p

Below I attached 2 screenshots from ArcGIS Pro.
When I type for ex. highway or any other field name from OSM data scheme -  it is fine, no errors, but when I try to write whole SQL expression, error shown below appears and also I cannot compose SQL expression using drop down list, because it is empty.

Any ideas how to write whole SQL expression?

Below error message:

  ERROR 000800  The value is not a member of OBJECTID | Shape | highway | building | natural | waterway | amenity | landuse | place | railway | boundary | power | leisure | man_made | shop | tourism | route | historic | aeroway | aerialway | barrier | military | geological | OSMID | osmuser | osmuid | osmvisible | osmversion | osmchangeset | osmtimestamp | osmMemberOf | osmTags | osmSupportingElement | osmMembers | Shape_Length.

EDIT:

Removed meaningless code, fixed bugs.
Now code is working properly.


texts:
I'm using the spatial extension for Talend (Ex-TOS SDI) and i'm experiencing a blocking issue. I'm trying to make what ESRI call a Spatial Join and I can't find a proper solution to do that ( transferring attributes from one shape to another based on the location of the features).

Does anyone tried (and succeed!) to make a spatial join on Talend?


texts:
I want to test if two raster images for the same scene taken in different moments are aligned. I would like to do this automatically as I will be working with a lot of scene pairs. The process would not need to be perfect I would like to spot obvious miss-alignments (i.e. more than 2 pixels) and flag them for manual inspection.

Is there a standard way to do so?

I am using python + GDAL


texts:
I want to show the attributes of the column "PHK_text" as a row on top of the table. "PHK_text" has ten different terms and all of them should be shown on top of the table as a new column. Now there are 339 objects in the table and in the end there should be 51 because the column "ortsteile_b" has 51 different terms. The new table should show all ten terms of "PHK_text" for all 51 terms of "ortsteil_b".
Which tools can I use to get a result?


texts: Saving vector memory layer to file using PyQGIS gives empty file?

texts:
I am learning to write standalone scripts with Pyqgis.

I am trying to read a file with a big shapefile of polylines (64MB), create a memory layer, copy the features to, and save into another file all the features of the memory layer.

The code looks ok, no errors, but I get an empty (no features) output file from memory layer.

Why ?  I can list the features into the memory layer, but they are not write to the output file.

This is the code:

#!/usr/bin/env python

from qgis.core import *

QgsApplication.setPrefixPath("/usr", True)
qgs = QgsApplication([], False)

qgs.initQgis()

# Origin vector layer into file
vlayer = QgsVectorLayer("BigShapefile64MB.shp", "bigShapeMB", "ogr")

feats = [feat for feat in vlayer.getFeatures()]

mem_layer = QgsVectorLayer("Lines?crs=epsg:4326&index=yes", "duplicated_layer", "memory")

# Add field and populate to memory
attr = vlayer.dataProvider().fields().toList()
mem_layer_data = mem_layer.dataProvider()
mem_layer_data.addAttributes(attr)
mem_layer.updateFields()
mem_layer_data.addFeatures(feats)

mem_layer.updateExtents()

# Check mem layer
print " Mem Layer features:", mem_layer_data.featureCount()
print "fields:", len(mem_layer_data.fields())
e = mem_layer_data.extent()
print "extent:", e.xMinimum(), e.yMinimum(), e.xMaximum(), e.yMaximum()

mem_layer.commitChanges()

crs=QgsCoordinateReferenceSystem("epsg:4326")

# Save memory layer to file
error = QgsVectorFileWriter.writeAsVectorFormat(mem_layer, "outShapefile.shp", "UTF-8", crs , "ESRI Shapefile")

if error == QgsVectorFileWriter.NoError:
    print "success! writing new memory layer"

qgs.exitQgis()

All ok, but outShapefile.shp has 0 features, is created, but empty, why ?

Note: (last line: qgs.exitQgis() give me a segmentation fault, but this not the problem, I think. )


texts:
First time, map loads and the locations are plotted properly.

Second time, when the new set of locations are pulled from the app, map reloads and gets the new set of data. But it is not refreshed. The old plotted points are not getting removed.


texts:
I have a line layer with segment lines, some of them mean a mainstream system that I'd like to keep. Others are isolated from the mainstream which I'd like to delete.
The mainstream is made up by splines, each line segment is connected to another line segment. The issue here is that I don't have any attribute/field that can be used to divide them.
I tried "select by location" to select the ones have common ending points, but that only delete part of the isolated segments. Because some of the isolated segments are connected by several lines.
Then I tried making a buffer of the lines, but "dissolve" never worked for me.
Is there a way to unite all connected lines to one?
Then I can select my mainstream by a single click.

'Polyline to raster' to create a raster.
'Region group', set neighbor number as 8, pixel size as 100 (very small)
Use a color ramp to check on the new region group raster.

But now I am lost.
How can I extract the mainstream part from this raster?
The attribute table looks like this:
Sort by Count

Sort by Value


texts:
Is it possible to import a .gfs file into PostGIS?

I can open the .gfs file in QGis, so I assume it is a geographic file of some kind: but I'm not sure how to get it into PostGIS.

I have Googled but can't find any examples of doing this with ogr2ogr. Do I need to convert it into some intermediate format first?


texts:
How can I estimate kernel density, taking barriers (i.e. a priori zero-probability regions) into account?

At the moment, regions where the population cannot occur are clipped out for display. However, I imagine this is biasing the results, since observations near the zero-probability regions don't have as much impact as they might.

I have several thousand observations, so a tool or programmatic solution is a necessity. Bonus points for an implementation in ArcGIS or R.


texts: `NameError: name `qgis` is not defined` in QGIS Plugins

texts:
In the QGIS (2.99) Python Console, the name/variable qgis seems to be automatically defined.

However, in writing a plugin, there is an error whenever I refer to qgis. For example, in I have qgis.iface.mapCanvas() in my Plugin's `initGui(self)' function, I get an error

NameError: name `qgis` is not defined

How can fix this error?

By the way, I am trying to put the following code from this GIS.SE answer(https://gis.stackexchange.com/a/45105) into my plugin

tool = PointTool(qgis.iface.mapCanvas())
qgis.iface.mapCanvas().setMapTool(tool)

The code works fine inside QGIS main window (in Python Console), but generates the above error when moving it into the Plugin. I can't figure out how to let the PointTool respond inside a plugin. The code for PointTool class from the linked answer is:

class PointTool(QgsMapTool):
    def __init__(self, canvas):
        QgsMapTool.__init__(self, canvas)
        self.canvas = canvas

    def canvasPressEvent(self, event):
        pass

    ...


texts: How to make QGIS Categorical Symbology Classify only the values that are present (i.e. not provide a category for "All other values")?

texts:
I am making a choropleth map to show the value, for each polygon, of an ordered categorical variable (a decile, in fact). So, each polygon already has a specific singular value of this variable, between 1 and 10. There are no missing or NA values. I am using an ordered palette (viridis).
When I use "Classify" within the Symbology tab, I get 11 categories, because QGIS appends an "All other values" (AOV) category. This category is allocated the far end of the viridis palette (yellow), and all other palette colours within the ramp are slightly shifted accordingly (as the palette is split into 11 rather than just 10). For consistency with other maps, I would like my 10 categories to be allocated evenly along the viridis palette. I do not want to have an "All other values" category - it has no use for me as I don't have any other values not already included in the range 1-10, or any missing values.
In order to get the proper 10-step colour ramp (for consistency with other maps), I have been deleting the AOV category and then manually amending all the other colour values to what they would have been if the palette had been divided up into 10 originally. This is time-consuming and frustrating.

Is there a way to prevent QGIS from appending the AOV category?
Or, is there a way to get QGIS to allocate a 10-split viridis colour ramp
to my values 1-10 and to allocate an entirely different shade (say,
grey) to the AOV category, so I could simply delete it without affecting the rest of the ramp?

If I make my map with the tmap package in R, I do not get this issue.
For those thinking I should use the original data and use the Graduated tab and create my own deciles from that, imagine that all I have is the already-calculated decile number for each area. While the source data would have been a continuous variable, once it is converted into quantiles, it becomes ordered categorical (or factor) data - unless I am mistaken? (A polygon is either in decile 5 or 6, say, it's not 5.78)


texts: How can I Transfer mysql wkt directly into a postgis table?

texts:
I have wkt geometry data in mysql stored in blob format.
I need to get this data into a postgis db and am looking for the quickest-easiest way to do it. I have managed to select the data "ASTEXT" but can then only copy it into QGIS -> export to shape -> import to postgis....but this is probably not the easiest way.

Is there a way to directly import a table from mysql into postgis?

thanks for any help,


texts:
I am working with GRASS and I am trying to do a supervised classification. The issue here is that I get an error for NaN values.

  ERROR: Not enough data points

Now I know that is an error that you get, most of the times, when the projection is different and because I know that I re-projected all of my data in GGRS87 (I have both vector and raster data) but I still get the same error. What should I do? Should I try another projection or is there another way to do supervised classification.

P.S I tried with this tutorial https://grasswiki.osgeo.org/wiki/Image_classification#Supervised_classification

and I tried to create a raster map from CORINE2000

g.region vect=corine align=raster
v.to.rast in=corine out=corine_training use=attr col=CODE labelcol=CODE2 --o

And now I can't generate statistics using

i.gensigset group=OB subgroup=raster_maps sig=corine_smap training=corine_training


texts: Create Mercator map with arbitrary center/orientation?

texts:
I want to see a Mercator map where I choose what latitude/longitude is
in the center and what direction is "north". In some sense, a Mercator
map of a rigid rotation of the globe.

Does such a program/site exist? I'm familiar with Mercator Rotator,
but this appears to be vaporware. I'm working on a project like this
myself, and want to see if I'm duplicating anything.

Same general question for other projections.


texts: How can I transfer data from one PostGIS database to another?

texts:
I have quite a bit of GIS data in a old PostGIS Database (version "1.1.6") and I need to transfer it to another PostGIS Database ("POSTGIS="1.5.3").

I tried creating a backup and restoring that on the new Database, but it throws several hundred errors.

What is the correct way of transferring the data?


texts: How can I find out if there's a traffic light in the route given an array of lat, long (GPS positions)?

texts:
I've the minutely (60 seconds) GPS trail of a vehicle travelling in USA. I want to find out if in the particular 60 second GPS trail there's a traffic LIGHTS/SIGNS or not.

Can I get this information from OSM database?

How do I go about doing this?


texts: Comments not getting saved in layer attribute

texts:
I have to send the shp files for the filling of the attributes to the field agents. To clarify the type of attribute to be collected I want to provide details in comment column of the Fields in Layer Properties.
The comment details are not being saved in project or qml.
I want to save comment in file other than project (like qml) so that I can send it along with the shp data and it guides the person filling its attribute. We are using QGIS for our project.


texts: Imposm does not create planet_osm_* tables

texts:
I'm trying to follow the Tilemill OSM-Bright tutorial on a fresh install of OSGeoLive 8.0.

I followed the instructions in the tutorial, and when I run tilemill, it can't access planet_osm_polygons.

Postgis Plugin: ERROR:  relation "planet_osm_polygon" does not exist
LINE 2:   FROM planet_osm_polygon                 ^

...

Looking over the console from imposm, no planet_osm* tables are mentioned. Was imposm supposed to create this table?

~/Downloads$ imposm -U user -d osm -m ../osm_bright/mapbox-osm-bright-5f237ac/imposm-mapping.py --read --write --optimize --deploy-production-tables new-york_new-york.osm.pbf
loading ../osm_bright/mapbox-osm-bright-5f237ac/imposm-mapping.py as mapping
password for user at localhost:
[14:23:43] ## reading new-york_new-york.osm.pbf
[14:23:43] coords: 5144k nodes: 102k ways: 734k relations: 5k (estimated)
[14:27:26] coords: 8216k nodes: 18k ways: 1355k relations: 3k
[14:27:29] reading took 3 m 45s
[...
[14:42:33] ## optimizing tables
Clustering table osm_new_mainroads
Clustering table osm_new_buildings

...


texts: Exporting shapefile from ArcGIS Desktop into Autodesk Ecotect?

texts:
I am trying to export a shapefile from ArcGIS 10 into Autodesk Ecotect Analysis 2011 for solar analysis. I have a shapefile with polygons, and the heights of buildings is stored in the attribute table. I tried saving it as a KML file but Ecotect does not import KML files. Instead, Ecotect would accept any DXF or 3D files (such as 3D Metafile and 3D Studio).

Is there a way of saving 3D files from GIS other than KML?

I tried importing KML into Sketchup and Autocad without success.

Any other idea for exporting shapefiles into a software that can extrude them and save as 3D would help too!


texts: OpenLayers Map Full Screen Option Not working

texts:
I want to add a full screen option in my map implemented using OpenLayers. I tried adding it with:

var map1 = new OpenLayers.Map({
        div: 'map',
        projection: new OpenLayers.Projection("EPSG:3857") ,
        controls: OpenLayers.control.defaults().extend([
        new OpenLayers.control.FullScreen()
        ]),
    });

and also by adding:

map1.addControl(new OpenLayers.control.FullScreen());

But it didn't work.

here is my fiddle:

http://jsfiddle.net/expedio/gmy0rpam/


texts:
I have videos with gps track (normally gpx file) and i need to put them together in the same final product, there in an add on "VIDEO UAV" in QGIS but it seems that he is not functioning well on the latest versions.
Can you help me please?


texts: Exporting data to a text file using a template

texts:
I've build a toolbox in my project environment to plan "Aids to Navigation".
The vessel (ship) that does the actual work in the field has a navigation computer that can only import waypoint.
These waypoints are stored in a *.wpt file, that is just a text file with a few lines per waypoint.

So to make this file: I've prepared a layer that grabs the geometry from my planning tool and displays the data as a HTML maptip in the correct manner. I show all the waypoints I need (1 by 1), select the text, copy the text and paste it into a text-editor. After pasting all the waypoints in one file I save it as a .wpt file.
This isn't a lot of work when designing small projects. But with 30+ mutations it isn't very efficient.

I am trying to come up with a manner to export all the data in one go to a text file in the *.wpt format. Can't find a plugin in QGIS or come up with a way to get this done


texts: Issue dsplaying GeoJSON file in OpenLayers 3.9.0

texts:
I'm trying to get vector feature from geojson file displayed in openlayers 3.9.0.

{"type": "FeatureCollection","crs": {"type": "name","properties": {"name": "EPSG:3857"}},
"features": [{"type": "Feature","geometry": {"type": "Point","coordinates": [0, 0]}}]}

Does anyone see why this code does not do it?

<!doctype html>

    <title>2</title>
    <link rel="stylesheet" href="http://openlayers.org/en/v3.9.0/css/ol.css" type="text/css">

Note that reading the geojson from an internally defined object works fine:

  var geojsonObject = {
    "type": "FeatureCollection",
    "crs": {"type": "name","properties": {"name": "EPSG:3857"}},
    "features": [
      {"type": "Feature","geometry": {"type": "Point","coordinates": [0, 0]}}
    ]
  }
  ...
  new ol.layer.Vector({
    source: new ol.source.Vector({
      features: (new ol.format.GeoJSON()).readFeatures(geojsonObject)
    })
  })


texts: Can an invidual sell ArcGIS Pro add-ins at ArcGIS Marketplace?

texts:
I have an add-in that I would like to sell to ArcGIS pro users. I checked the ESRI website and they have only mentioned ESRI Network partners as eligible providers. Can an indvidual not list products at ArcGIS Marketplace?

Who are ArcGIS Marketplace program providers?
Providers are Esri Partner Network partners, including Startup and NGO partners and Distributors, in good standing. Contact ArcGIS Marketplace Resource Center for information and to apply to become a provider. https://doc.arcgis.com/en/marketplace/provider/faq.htm#anchor2


texts: Joining 2 layers by location doesn't join desired fields from attributes tables

texts:
I have 2 layers that I'm trying to join by location: 1 is a point layer (CRS 4326) and the other is a shapefile layer (CRS 4326). The shapefile layer has a field in the attribute table whose values I want assigned to the points that fall within/along the boundaries of their respective shapes. (i.e. each point is a lat/lon coordinate, and the shapefile layer is a bunch of property parcels--some of which contain the points).

But when I join them by location the output layer doesn't include any of the fields from the shapefile layer. I just get new fields with sum, mean, max, min etc.

Is there a way to make this happen? I've tried this with QGIS 2.18 and QGIS3 without success.


texts: Unable to get Latitude in the Y field in ArcMap

texts:
I have an issue, and I don't know to solve it. I've downloaded a CSV file which contains among other metadata latitude and longitude. I want to display them on the map. But I want to display that by right-clicking and selecting display XY; I get only longitude in all X and Y fields. I've tried to rearrange the columns manually by dragging or by choosing the fields options in the table properties.

The reference system is a GCS WGS 84.

Is there any other way to fix it?


texts:
I'm using QGIS and MMQGIS and I get this error:

Traceback (most recent call last):
              File "C:/Users/Ana Eljall/.qgis2/python/plugins\mmqgis\mmqgis_dialogs.py", line 956, in run
                message = mmqgis_geometry_convert(self.iface, layername, newtype, splitnodes, savename, 1)
              File "C:/Users/Ana Eljall/.qgis2/python/plugins\mmqgis\mmqgis_library.py", line 2398, in mmqgis_geometry_convert
                if (feature.geometry().wkbType() == QGis.WKBPoint) or \
            AttributeError: 'NoneType' object has no attribute 'wkbType'


texts: Im working on project for calculating pedestrian and bike accessibility to public services.For pedestrian I already could produce Isochrone maps ( by applying the street slops) now I want to calculate and produce the service area (Isochrone) for bikes by applying both street slope and also street intersections.I mean applying street intersection as travel cost. I don't know how to do it.


texts: Requested JSON parse failed while getGetFeatureInfoUrl

texts:
I am trying to get the feature attributes from WMS service published in GeoServer and showing the attributes as a popup in the .NET application.

When I make an Ajax call it gives "requested JSON parse failed" exception, However, the URL works and returns JSON data in a browser.

PFB code block:

 var url = wmsSource.getGetFeatureInfoUrl(
        evt.coordinate, viewResolution, 'EPSG:3857',
        { 'INFO_FORMAT': 'jsonp' }); //text/html

      if (url) {
          var parser = new ol.format.GeoJSON();
          var lookup = {};
          $.ajax({
              url: url,
              dataType: 'jsonp',
              method: 'GET',
              contentType: 'application/json; charset=utf-8',
              success: function (response) {
                              alert("hI");
                          },
              error:  function(jqXHR, exception) {
                          if (jqXHR.status === 0) {
                              alert('Not connect.\n Verify Network.');
                          } else if (jqXHR.status == 404) {
                              alert('Requested page not found. [404]');
                          } else if (jqXHR.status == 500) {
                              alert('Internal Server Error [500].');
                          } else if (exception === 'parsererror') {
                              alert('Requested JSON parse failed.');
                          } else if (exception === 'timeout') {
                              alert('Time out error.');
                          } else if (exception === 'abort') {
                              alert('Ajax request aborted.');
                          } else {
                              alert('Uncaught Error.\n' + jqXHR.responseText);
                          }
                      }
          });

      }

The JSON response on browser from the URL is:

{"type":"FeatureCollection","features":[{"type":"Feature","id":"village_maps.115","geometry":{"type":"MultiPolygon","coordinates":[[[[79.13895583,25.3233091],[79.13482866,25.32499247],[79.13317419,25.32360682],[79.12904055,25.32355175],[79.11112382,25.32624913],[79.11054308,25.32213731],[79.10696928,25.32399947],[79.1034709,25.32324635],[79.09974985,25.31942675],[79.09800786,25.32022119],[79.0979736,25.31867895],[79.09166887,25.31635774],[79.08927284,25.31821429],[79.08661964,25.31754376],[79.07432635,25.32064035],[79.07417549,25.32250699],[79.07728205,25.32362325],[79.08025674,25.32743095],[79.08831236,25.32775063],[79.09032527,25.33283685],[79.08757208,25.33551219],[79.07363746,25.34013379],[79.07406066,25.34194949],[79.07257536,25.34266333],[79.07427301,25.34407714],[79.07571542,25.35252829],[79.08648018,25.35161346],[79.09057063,25.35709824],[79.09046858,25.35497181],[79.08648332,25.34991794],[79.09656147,25.34616603],[79.09976983,25.35526709],[79.10152538,25.3574017],[79.10472238,25.35755957],[79.10617664,25.36085029],[79.11197018,25.35874088],[79.11217002,25.35463277],[79.11563481,25.35019114],[79.11465262,25.348675],[79.11888327,25.34534672],[79.12336584,25.34908317],[79.12612299,25.34570295],[79.13045243,25.34409387],[79.13712699,25.34887448],[79.13983827,25.34651654],[79.13916554,25.34453145],[79.14239623,25.34280575],[79.14273195,25.33875767],[79.1448726,25.33784395],[79.14887941,25.33979993],[79.1535527,25.33947772],[79.15312606,25.33435852],[79.15454125,25.33390489],[79.15492136,25.330539],[79.15069001,25.3216973],[79.14512387,25.32357109],[79.13895583,25.3233091]]]]},"geometry_name":"geom","properties":{"objectid_1":115,"objectid_2":115,"objectid":115,"shape_leng":0.313447905145,"pipeline_l":2.72376,"shape_le_1":0.313447907911,"shape_area":0.00194710986264}}],"totalFeatures":"unknown","numberReturned":1,"timeStamp":"2019-08-21T05:33:08.387Z","crs":{"type":"name","properties":{"name":"urn:ogc:def:crs:EPSG::3857"}}}


texts: Same behaviour for minzoom for google maps api and leaflet

texts:
I need to implement the same behavior for leaflet and google maps for the case of minzoom. Google maps duplicates all the objects from the main map to the world`s copies, however leaflet doesn't since leaflet range is not limited between -180 180 . So I need either prevent coping on google maps or clone every objects on leaflet.

There is an option optimized for markers in the google api, but other objects types (polygons, polylines) do not support this option.
And for the leaflet there is a way to set multiple coordinates manually (or may be some plugin exists?).

What is a suitable solution for this issue?


texts: openlayers 4: force refresh tiles in ol.source.VectorTile on interval

texts:
I'm using openlayers 4, and have a ol.source.VectorTile which uses MVT tiles. I want to refresh the tiles every 5 seconds without refreshing the browser.
I've used source.changed(), source.refresh(), layer.changed() and map.renderSync()

But none of them not working for me, I mean tiles not refreshed after calling any of them.

Can anyone help me please.


texts: How to create sequential field based on sorted data? (or extract pop-up data?)

texts:
I downloaded a kmz field and converted it to a GDB fc on arcmap.

However, the field I need to label the layer by is stored in the pop-up html field...

It looks like this:

<META http-equiv="Content-Type" content="text/html">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">

Benchmark: 100

Point ID
100

Type
Vertical

Northing
2119938

Easting
533366

Elevation (NAVD88)
6.1 ft

Elevation (NGVD29)
7.18 ft

All I need is the benchmark number (Benchmark: 100) from each field

Is there a way to either extract fields from the popup field OR create a new field based on the sorted popup field?

Tried a bit of the code here: How to populate Unique ID field after sorting?

but did not get anywhere.


texts: Re Exporting complete data from Google Earth to QGIS

texts: I am trying to export Google Earth searches as a KML file to QGIS. How come all existing locations of a search (e.g. Walmart) in the United States only shows select locations. I'm doing a project on food deserts and want to include all major chains and their locations in the U.S.


texts: how to do horton analyis on single catchment in grass gis

texts:
I have generated streams using r.stream.extract and used r.stream.order
then I have tun r.stream.basin using inputs from r.stream.order. I am getting multiple basins. How to select a single basin and calculate the horton statistics on single catchment area.


texts: Conditional Statement (CASE-WHEN) in QGIS for slope aspect

texts:
After processing a DEM with a mask to get the aspect values of some polygons I would like to create a new field classifying those values as N, NE, E, SE, S, SW, W, NW.

I've written this but I can't get what I'm looking for:

CASE
WHEN  "Aspect"  < 22.5  THEN  'N'
WHEN  (22.5 >= "Aspect" AND "Aspect" < 67.5) THEN  'NE'
WHEN  (67.5 >= "Aspect" AND "Aspect" < 112.5) THEN  'E'
WHEN  (112.5 >= "Aspect" AND "Aspect" < 157.5) THEN  'SE'
WHEN  (157.5 >= "Aspect" AND "Aspect" < 202.5) THEN  'S'
WHEN  (202.5 >= "Aspect" AND "Aspect" < 247.5) THEN  'SW'
WHEN  (247.5 >= "Aspect" AND "Aspect" < 292.5) THEN  'W'
WHEN  (292.5 >= "Aspect" AND "Aspect" < 337.5) THEN  'NW'
WHEN  "Aspect"  >= 337.5  THEN  'N'
END

Is it possible to do this using CASE-WHEN, or is it better to use another kind of conditional?

I don't know how to write ranges using conditional sentences without getting an error except this way.


texts: Returning a new vector layer with additional setup from a processing algorithm in QGIS 3

texts:
I am working on a processing script in QGIS 3, which does the following:

Take a few inputs
Run a native algorithm on these inputs, saving the resulting new layer in memory
Set up the symbology for the layer (specifically, setting the renderer via layer.setRenderer(renderer_instance))
Return the layer, retaining its symbology, as the primary algorithm output

The key here is that I want to return a layer after setting up its symbology.
This is problematic because I'd like to use the common input parameter QgsProcessingParameterFeatureSink, since this lets the user choose how to output the results (memory layer, temporary file, saved in a format of their choice, etc).

However, as far as I can tell, the ...ParameterFeatureSink parameter is primarily used with QgsProcessingAlgorithm's parameterAsFeatureSink method, which returns an interface to a new layer, which does not allow setting up symbology. Also, it appears that this layer is what the algorithm expects to be returned at the end, and that returning another layer (like the one  described in the list above) in its place causes an error.

My question is this: how do I utilize the strengths of QgsProcessingParameterFeatureSink, while still being able to return a layer with its symbology set up as described?

The ideal solution here would be to use another of the parameterAs methods to get a place or name for the layer I'll eventually return. I could feed that name to the internal algorithm as its OUTPUT parameter, making it return a layer that my algorithm could subsequently return without generating an error.


texts:
I downloaded the 3.10 version of QGIS onto my Mac (version 10.15.6). There are no Plugins available in the programme. Do I have to download plugins from somewhere else?
Also, in the Raster Menu, the only 2 options I have are Raster Calculator & Align Rasters, no other options. I'd like to build a virtual raster.


texts: How to choose pixel values according to each land cover class?

texts:
I would like to do following steps using two types of images (precpitation and Land cover image):

select pixel values of precipitation image within each land cover classes
calculate total amount of precipitation of each land cover classes

Could you please tell me the tools I use and the steps I need to follow in ArcGIS, ENVI, or ERDAS.


texts: Understanding difference between GeoServer and OpenLayers?

texts:
I recently installed GeoServer and downloaded OpenLayers also. I am a bit lost about which one to use when it comes to add functionality in my web map.

For example, I want to add layer switcher, shall I try to do it from GeoServer or through OpenLayers?

Another example is opacity, shall I try to do it from GeoServer or through OpenLayers?

Does GeoServer stop when layers are imported and displayed as WMS and then OpenLayers take full responsibility for Controls, opacity, thematic maps, base layer etc?


texts: Why is the heatmap color not showing in TileMill?

texts:
I created a heatmap on QGIS (2.8) and styled it using singleband pseudocolor.  I want to turn it into map tiles, so I added the TIFF as a layer in TileMill. However, when it's in TileMill, it loses the style settings from QGIS and becomes black and white. I'm not entirely sure if this is an issue with something I'm doing on QGIS or on TileMill.

I just started using QGIS a few days ago.


texts: Importing raster files warning and extracting covariates error with 'raster' and 'amt' packages in R

texts:
I am using 'amt' package in RStudio, trying to assess habitat selection of a species in an area with telemetry data. My two preliminary problems are:

When I import the land data rasters (elevation and land cover) I get this warning message:

Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition

Is this something to worry about? Since all the individuals' shapefiles and the land data files are in epsg:2100 (Greek Grid) already before importing them in RStudio? I also tried importing them in epsg:4326 (WGS84) but I get the same warning, and also transforming them with transform_coords() function but again get the same warning. Moreover, I tried using the proj4 for the projection that I use, but I still get the same warning.
landcover = projectRaster(landcover, crs="+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs")

#or

landcover = crs("+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs")

When I try to extract covariates from the raster files, I get this error, which I have not figured out how to fix.

Error in UseMethod("extract_covariates", x) :
  no applicable method for 'extract_covariates' applied to an object of class
 "c('RasterStack', 'Raster', 'RasterStackBrick', 'BasicRaster')"

I have already checked the .tif files in ArcMap and they have identical dimensions, and the pixels have values (elevation and land cover of each pixel).
Any suggestions?
Code:
library(amt)
library(raster)
library(tidyverse)
library(lubridate)
library(sf)

#Import the land_cover and elevation rasters
landcover = raster("landcover.tif")

#Warning message:
#In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
#  Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition

elevation = raster("elevation.tif")

#Warning message:
#In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
#  Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition

map <- stack(landcover, elevation)
plot(map) #I can plot and see the maps
extract_covariates(map, where = "both") #Also tried to extract covariates without using stack()

#Error in UseMethod("extract_covariates", x) :
#  no applicable method for 'extract_covariates' applied to an object of class
# "c('RasterStack', 'Raster', 'RasterStackBrick', 'BasicRaster')"


texts: Using MapProxy to convert WMTS to WMS or TMS

texts:
I am trying to configure MapProxy to enable me to access a WMTS service via WMS or TMS.  Here is the source WMTS service:

http://tiles.arcgis.com/tiles/ovln19YRWV44nBqV/arcgis/rest/services/ScottIA_2014Aerial/MapServer

Here is my MapProxy instance:

http://104.197.62.82:8080/demo/

Here is my MapProxy configuration file:

services:
   demo:
   tms:
   wms:
layers:
 - name: scott_layer
   title: Scott
   sources: [scott_cache]
caches:
   scott_cache:
      grids: [scott_grid]
      sources: [scott_source]
sources:
   scott_source:
      type: tile
      grid: scott_grid
      url: http://tiles.arcgis.com/tiles/ovln19YRWV44nBqV/arcgis/rest/services/ScottIA_2014Aerial/MapServer/tile/%(z)s/%(y)s/%(x)s
grids:
   scott_grid:
      origin: 'nw'
      bbox: [2341500.5208333135, 533999.4791666567, 2517000.5208333135, 662999.4791666567]
      srs: 'EPSG:3426'
      res: [
         #  res            level     scale @90.7 DPI
          868.0555554189, #  0     3100198.41221019
          434.0277777094, #  1     1550099.20610510
          217.0138888547, #  2      775049.60305255
          108.5069444274, #  3      387524.80152627
           55.5555555468, #  4      198412.69838145
           27.7777777734, #  5       99206.34919073
           13.8888888867, #  6       49603.17459536
            6.9444444434, #  7       24801.58729768
            4.1666666660, #  8       14880.95237861
            2.0833333330, #  9        7440.47618930
            1.0416666665, # 10        3720.23809465
            0.5208333333, # 11        1860.11904733
      ]

If you use the MapProxy demo link above, you will see that nothing is showing up on the map.

Any idea what I am doing wrong?

Update:

I ultimately want to display the source map in Google Maps.  I can do that if the source map is a tile service that uses the standard Google Maps tiling scheme.  I can also do that if the source map is WMS.  The challenge here is that this source map uses a custom tiling scheme (see http://tiles.arcgis.com/tiles/ovln19YRWV44nBqV/arcgis/rest/services/ScottIA_2014Aerial/MapServer/WMTS/1.0.0/WMTSCapabilities.xml).  Therefore, I cannot find a way to display the tiles in Google Maps.  I was hoping to use MapProxy to convert this WMTS source map with a custom tiling scheme into WMS or TMS (TMS uses a tiling scheme compatible with Google Maps).


texts: Displaying single label for multiple features with same value in ArcMap?

texts:
I'm working on ArcMap 10. I have a layer of points which I categorized based on a field that we can call "Name". Thus, as you can see in the picture below, the layer is now composed of 24 colors, each of which for points with the same value of "Name".
I want to show the labels for this field, but instead of having a label for each point (which as you can assume wouldn't be appropriate given the amount of features) I would have a label for each group of points with the same color (same "Name" value), maybe in the barycenter of the points or linked to the barycenter by an arrow.
Is there a way to accomplish this?


texts: How to set a "no data" value to use in Zonal Statistics?

texts:
I have a raster with the thickness of an aquifer confining unit and I'm using zonal statistics to calculate the average thickness under a buffer zone around well points. However, in areas where the confining unit isn't present, there is just no data rather than a 0 value.

How can I change this raster so that zonal statistics will include the area with no thickness data as a "0" datapoint?


texts: How do I Spatially Partition a whole lot of points?

texts:
I have location data of customers, with over 130 million records, spread all across the country. This is a Spatial Table in PostGIS

I want to now partition the entire country into 'service areas', such that each service area has:

Atleast 30 Customers
Is at the most 1000 sq km.

I have researched ArcGIS's Aggregate tool, but it does not meet my criteria, since:

It does not take minimum count into consideration
I only have Qgis & Grass at my disposal.

What tool/process can I use to Achieve this?


texts:
I have a checkbox button, and I am trying to change the variable getStyle dynamically, when jQuery change event is triggered. Though style property doesn't get updated when the event is triggered.

var style1= new ol.style.Style({
      fill: new ol.style.Fill({
        color: 'khaki'
      }),
       stroke: new ol.style.Stroke({
          color: '#880000',
          width: 2,
          opacity:1
        })
    });

    var style2= new ol.style.Style({
        fill: new ol.style.Fill({
        color: '#5755d900'
      }),
       stroke: new ol.style.Stroke({
          color: 'yellow',
          width: 2,
          opacity:1
        })
    });

     var getStyle = style1; //default style

    $("#layerbg").change(function () {
      if ($(this).is(':checked')){
        getStyle = style2; //change style when checkbox is ticked
      }
      else{
        getStyle = style1;
      }
    });

    var vector = new ol.layer.Vector({
      style:getStyle, // this variable (getStyle) is not updating
      source: new ol.source.Vector({
      url: 'JS/getjson.json',
      format: new ol.format.GeoJSON()
       }),opacity: 1
    });


texts: Defining polygons based on thickness maps using ArcGIS for Desktop?

texts:
I want to find and create polygons that demarcates three different inflection points in a thickness raster: where the thickness start to increase, where the thickness is at its max and where it starts to flatten again (i.e. increase from the opposite direction). Please see the different interpolation lines in the attached figure. Warm colors define thin areas.
I want to do this quantitatively, instead of qualitatively drawing polygons by hand as in the example.

I'm using ArcGIS 10.2


texts: Getting corners of polygon in latitude and longitude using PostGIS?

texts: How can I get the corners of a polygon -a square/rectangle, to be precise- in latitude-longitude format by using its geom value -or some other way, using PostGIS?


texts: End execution of script and send email if schema has changed

texts:
I would like to append the data from one feature class to another after confirming that the schemas match in both feature classes.

I can append the data from one FC to another by using this simple code:

    arcpy.DeleteFeatures_management(outFeatureClass)
    arcpy.Append_management("FC", outFeatureClass, "TEST", "")

...but what I would like to do is stop execution of this code completely if the schemas have changed...in my example above, the data has already been deleted by the time that the schema change parameter ("TEST) fails.  I am using 'Append...' b/c these feature classes are connected to map services.

Ideally, I would like to first compare the schemas, and then continue to delete/append if they match.  If they do not, I would like to send an email, and end the script.

I have tried using the FeatureCompare tool, which will print a message if the schemas do not match, but I haven't figured out if or how I can use the data in this message as a variable in my code:

 compare_result = arcpy.FeatureCompare_management(inFC, outFC, "OBJECTID", "SCHEMA_ONLY", "","","","","","","CONTINUE_COMPARE")
 print arcpy.GetMessages()

When I run this code, one line of the output message reads:

  Table: Tables have different number of fields (Base: 48, Test: 47).

...if there has been a change, or

  Field: Field properties are the same.

...if the schemas match.

Is there a way I can somehow use this information to set a flag or boolean in my code?

I also considered using a SearchCursor, but I am not so sure that is the right idea.

I am currently able to send an email confirming that the script was successful or if there was an error, but I would prefer to send an email before an error has already resulted in a faulty map service.


texts:
Has anybody had any luck getting the tools in the QGIS2 Toolbox to run? So far I have tried SAGA and TauDEM but both throw up errors. I eventually got TauDEM to produce something from the D8 tool by feeding it empty rasters to fill but the results were complete nonsense. Outputting to temporary files did not help at all.
It's probably yet another path issue, GIS on OSX seems to be one long series of path issues.


texts: Accumulated value from many rasters python/arcgis10

texts:
I need to calculate an accumulated value of temperature (as an additional raster) from many rasters.

Is there any option to do this in ArcGIS 10 or any hints for a python script?


texts: Transformation from local projection to WGS84

texts:
I'd like to reproject point that is in local coordinate system to the WGS84, using gdal. But exception "No PROJ.4 translation for source SRS, coordinate transformation initialization has failed" appears. Unfortunately, I've no idea how to proceed with the issue. Here is my code on c#

            Gdal.AllRegister();
            Dataset dataset = Gdal.Open(@"D:\TifExample\raster.tif", Access.GA_ReadOnly);
            String WKTFromTif = dataset.GetProjectionRef();

            Double[] gt = new double[6];
            dataset.GetGeoTransform(gt);
            Int32 Rows = dataset.RasterYSize;
            Int32 Cols = dataset.RasterXSize;
            Double upperLeftX = gt[0];
            Double upperLeftY = gt[3];
            Double lowerRightX = gt[0] + Cols * gt[1] + Rows * gt[2];
            Double lowerRightY = gt[3] + Cols * gt[4] + Rows * gt[5];

            Double[] upperLeftPoint = { upperLeftX, upperLeftY };
            Double[] lowerRightPoint = { lowerRightX, lowerRightY };

            SpatialReference currentSpatialReference = new SpatialReference(WKTFromTif);

            String EPSG4326WKT = "GEOGCS[\"WGS84 datum, Latitude-Longitude; Degrees\", DATUM[\"WGS_1984\", SPHEROID[\"World Geodetic System of 1984, GEM 10C\",6378137,298.257223563, AUTHORITY[\"EPSG\",\"7030\"]], AUTHORITY[\"EPSG\",\"6326\"]], PRIMEM[\"Greenwich\",0], UNIT[\"degree\",0.0174532925199433], AUTHORITY[\"EPSG\",\"4326\"]]";
            SpatialReference newSpatialReference = new SpatialReference(EPSG4326WKT);

            CoordinateTransformation coordinateTransform = new CoordinateTransformation(currentSpatialReference, newSpatialReference);

            coordinateTransform.TransformPoint(upperLeftPoint);
            coordinateTransform.TransformPoint(lowerRightPoint);

P.S. WKTFromTif has next value for my raster :"LOCAL_CS[\"WGS 84 / Pseudo-Mercator\",GEOGCS[\"WGS 84\",DATUM[\"unknown\",SPHEROID[\"unretrievable - using WGS84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],AUTHORITY[\"EPSG\",\"3857\"],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]"

Does anyone knows how to solve the problem?

Here is the data from gdalinfo

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (4042799.739385183900000,6382762.252179093700000)


texts: setting the style of each leaflet feature

texts:
I am using leaflet and geojson of some regions to show the related calculations about those regions on map by color.

I calculated each region value by a function which is called calculation(). The output of this function is a simple float. Here is what I have done:

I use onEachFeature for loading each feature.
Then I use Ajax to do some calculations. (I use Jquery here)
I save the calculation results in an array.
then, I use .done to set the style of each feature (which is a region) to set its color.

Here is my code:

var region=L.geoJSON(qom, {onEachFeature:
function (feature, layer) {
    $.ajax({    //send the ajax with following parameter
        type:"POST",
        timeout: 2000,
        cache:false,
        url:    "../tool/calculation.php", //the php we are sending info too and has got
        data: { },    // multiple data sent using ajax
        success: function (data) {    //now we have data
            $.each(data, function(calculation_related_key,calculation_related_content) {
                region_array_calculation.push([the_region_id, calculation_related_content['result']]);
            });
        }
    }).done(function (){
        L.geoJSON(qom, {style : set_feature_style(feature, region_array_calculation)}).addTo(mymap);
        });
}
});

Here is setting style code:

function set_feature_style (feature, calculation_array) {
for (i=0;i<calculation_array.length;i++) {
    if (calculation_array[i][0]==feature.properties.region_id) {
        return {
            weight: 0.5,
            opacity: 0.3,
            color: 'Black',
            dashArray: '3',
            fillOpacity: 0.2,
            fillColor:get_feature_color(calculation_array[i][1])
        }
    }
}
}

Here I set the color:

function get_feature_color(input) {
var x= input;
switch (true)   {
case (x>=0 && x<=0.5):
    return 'blue';
    break;
case (x>0.5 && x<=1):
    return 'green';
}
}

But unfortunately the result I get is all in blue.

However the calculation for some of the regions is more than 0.5.


texts: Exponential function variogram sill values

texts:
I am using gstat package for deriving a variogram using exponential model.
My R code goes like this:

g <- gstat(id= "HV",formula=HV~X+Y,data=vari)
vari.var <- variogram(g)
vari.var
plot(vari.var,col="black",pch=16,cex=1.3,xlab="Distance",ylab="Semivariance",main="Omnidirectional variogram NDVI")
model.vari <- fit.variogram(vari.var,vgm("Exp",range=range,nugget=0),fit.method=2)
plot(vari.var,model=model.vari,col="black",pch=16,cex=1.3,main="Omnidirectional variogram NDVI").

I am not sure whether it is correct or not.

Also I am not sure how to calculate sill values.


texts: TypeError: 'int' object does not support item assignment

texts:
I am trying python for the first time. I have a raster layer of Type Float32, and another raster output of Type Integer created from this raster layer outside the python code. I want to scan every column wise, and pick up the Float raster minimum value location within a 5 neighbourhood of Integer raster value of 1 in every column and assign the value 1 at this minimum value location and all other pixels as Null.

I am getting the Type Error : 'int' object does not support item assignment

I am giving the python code below:

new_array = [ [0 for i in range(band.XSize)] for j in range(band.YSize)]

for i, item in enumerate(myarray1):
    for j, element in enumerate(item):
        if (myarray2[i][j] == 1):
            if myarray1[i][j] < myarray1[i-2][j]:
                new_array[i-2][j] = -999
                if myarray1[i][j] < myarray1[i-1][j]:
                    new_array[i-1][j] = -999
                    if myarray1[i][j] < myarray1[i+1][j]:
                        if myarray1[i][j] < myarray1[i+2][j]:
                            new_array[i][j] = 1
                        else:
                            new_array[i+2][j] = 1
                            new_array[i][j] = -999
                    else:
                        if myarray1[i+1][j] < myarray1[i+2][j]:
                            new_array[i+1][j] = 1
                            new_array[i][j] = -999
                        else:
                            new_array[i+2][j] = 1
                            new_array[i][j] = -999
                else:
                    if myarray1[i-1][j] < myarray1[i+1][j]:
                        if myarray1[i-1][j] < myarray1[i+2][j]:
                            new_array[i-1][j] = 1
                            new_array[i][j] = -999
                        else:
                            new_array[i+2][j] = 1
                            new_array[i][j] = -999
                    else:
                        if myarray1[i+1][j] < myarray1[i+2][j]:
                            new_array[i+1][j] = 1
                            new_array[i][j] = -999
                        else:
                            new_array[i+2][j] = 1
                            new_array[i][j] = -999
            else:
                if myarray1[i-2][j] < myarray1[i-1][j]:
                    if myarray1[i-2][j] < myarray1[i+1][j]:
                        if myarray1[i-2][j] < myarray1[i+2][j]:
                            #new_array[i-2][j] = 1
                            new_array[i][j] = -999
                        else:
                            new_array[i+2][j] = 1
                            new_array[i][j] = -999
                    else:
                        if myarray1[i+1][j] < myarray1[i+2][j]:
                            new_array[i+1][j] = 1
                            new_array[i][j] = -999
                        else:
                            new_array[i+2][j] = 1
                            new_array[i][j] = -999
                else:
                    if myarray1[i-1][j] < myarray1[i+1][j]:
                        if myarray1[i-1][j] < myarray1[i+2][j]:
                            new_array[i-1] = 1
                            #new_array[i][j] = -999
                        else:
                            new_array[i+2][j] = 1
                            new_array[i][j] = -999
                    else:
                        if myarray1[i+1][j] < myarray1[i+2][j]:
                            new_array[i+1] = 1
                            new_array[i][j] = -999
                        else:
                            new_array[i+2][j] = 1
                            new_array[i][j] = -999
        else:
            new_array[i][j] = -999


texts: How to monitor map service connections ArcGIS Service?

texts: I would like to see the historical data on the number of users that have connected to a map or feature service in ArcGIS Server. Is this available in the logs?  If not, is there a way to do execute this task?  I am using ArcGIS Server 10.


texts: In Google Earth Engine error occur "List (Error) Element.get: Parameter 'object' is required."

texts:
I got error as "List (Error) Element.get: Parameter 'object' is required" In my Code. Code link is : https://code.earthengine.google.com/1d7b5f17c84ab5b5a2ecc960b38d771e. For simplicity purpose I have reduced the code and same error occurs.
var table = ee.FeatureCollection("users/rakeshsahu1609/subdiv_himalaya_v120319"),
L8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA"),
SRTM = ee.Image("USGS/SRTMGL1_003"),
table2 = ee.FeatureCollection("users/rakeshsahu1609/14_rgi60_SouthAsiaWest");

var ROI8 = L8.filterBounds(table).filter(ee.Filter.calendarRange(2019,2019,'year'))
                            .filter(ee.Filter.calendarRange(7,10,'month'))
                            .filterMetadata('CLOUD_COVER', 'less_than', 10);

var ROI = ROI8

var distinctRows = ee.List(ROI.distinct(['WRS_ROW']).aggregate_array('WRS_ROW'));
var distinctPaths = ee.List(ROI.distinct(['WRS_PATH']).aggregate_array('WRS_PATH'));

var imagePerPath = distinctPaths.map(function(path){
var imagePerRow = distinctRows.map(function(row){
var images = ROI.filter(ee.Filter.and(ee.Filter.eq('WRS_ROW', row), ee.Filter.eq('WRS_PATH', path)));
var Cloud_pre_images = images.map(
     function(img)
   {
      var pixelAreas1 = img.multiply(ee.Image.pixelArea().divide(1000000));
        var Areas1 = pixelAreas1.reduceRegion(
                   {
                     reducer: ee.Reducer.sum(),
                     geometry: img.geometry().buffer(100),
                     scale: 30,
                     maxPixels: 1E13

                    }).get('B1');

        var cloude_per = ee.Number(Areas1)

        return (ee.Image.constant(cloude_per))
                         .copyProperties(pixelAreas1, pixelAreas1.propertyNames())
                        .set("Cloud_per", cloude_per);
        });
      var img_sort =Cloud_pre_images.sort('Cloud_per').first();

     return images.filter(ee.Filter.eq('system:index', img_sort.get('system:index')))
   });
  return imagePerRow;
 });

 var leastCloudies = ee.ImageCollection.fromImages(imagePerPath.flatten()).mosaic();
 Map.addLayer(leastCloudies.clip(table))
 print(imagePerPath);


texts: Kriging observations using additional knowledge from geophysical model

texts:
Here is a simplified explanation of the problem I am trying to solve : I have (1) a 100x100 grid with air pollutant concentrations in a city obtained from an air quality geophysical model, and (2) a set of observations of these air pollutant concentrations at 10 stations within the city (i.e. in 10 different cells of my domain).
The map coming from my air quality geophysical model shows a realistic (physically-consistent) spatial distribution of the air pollutant concentrations but is prone to biases/errors (like any other geophysical model). So I want to correct this spatial distribution using the observations.
I read various posts but I am a beginner in geo-statistics and I am still quite confused. I understood that kriging methods could be used to estimate an air pollution map solely based on the information provided by the 10 observations but in my case, I want to use the extra information provided by the air quality model (since I want the output to respect more or less the spatial distribution predicted by the air quality model, showing roads, parks...).
How can I do that? From what I read (but poorly understood...), it seems that universal kriging might be the solution but I am still quite confused on where the information from the air quality model might be taken into account. Should this extra information be included as what we call the external drift?
NOTA : I think that the Python function UniversalKriging from the PyKrig package might be the solution for me (https://geostat-framework.readthedocs.io/projects/pykrige/en/stable/generated/pykrige.uk.UniversalKriging.html) but I am still not sure how to properly use this universal kriging method and this function specifically. In any case, I am open to other recommendations of packages.


texts: ArcGIS REST map server query doesn't return geometry as JSON?

texts:
I use Python code to query an ArcGIS map server. The result link (as json) only returns the Out Fields not the geometry. When I get the result link as kmz it seems to have the geometry, as I'm able to convert into the following layer:

Which appropriately looks like a subset of the service layer:

I've used the same query method successfully on other map servers (e.g. tiger for US Census), any idea what I'm doing wrong?


texts: Adding scripts within CartoDB infowindow?

texts:
i am trying to implement a simple "zoom to this point" link within a CartoDB infowindow. To do this i figure i need:

To be able to place a script within an infowindow (mustache) template, so that the zoom link can be dynamically set
Access to the lat/lng value that was clicked

Its not clear to me whether the JS should be placed within the mustache template (nested  tags seem a bit messy), or elsewhere - something like an "infowindow.complete" event would be ideal, but not sure if that exists.

Currently i have

</script>

My script does run, but i dont know how to access the lat/lng point. Is the information available within content.data? How can i debug what is in {{content}} ?

This post seems to provide an answer (ie use the 'featureClick' event, which has the lat/lng passed in), however the code doesn't fire for me - presumably because the content in the infowindow is not completed yet (hence my suggestion above)


texts: Is there an existing way to do linear referencing in ArcGIS JavaScript API 2.1?

texts:
First of all, let me specify that this is for the JavaScript API 2.1 or later only, not ADF/Google/Virtual Earth/Silverlight/Flex or any of the other extensions.

We've been creating geometry from a road geometry service to display travel routes on our map. We also have a layer of known routes. We're working on trying to speed things up by using the known routes in a JSAPI map, and calling the road geometry service only when we need to update or add routes.

Using the road geometry service, we can draw the route so that when a vehicle is shown on the route, the road already traveled can be red, and the remaining route as green. Has anybody tried to render lines with different symbols on different parts of the line?


texts: How to connect FME user attributes from reader to writer when there are transformers in between

texts:
I have the following schema (figure) in which a reader takes an excel file and writes something into a database.
In between I have added several transformers.
For some reason I can not pass the user attributes I need from the reader to the writer (red arrows in the writer).
How can I solve this?


texts:
I know there are related questions like this one, this question, and this one all over this site. I have tried the many possible solutions I have seen at the point.

The problem: I am trying to make a basic overlay of 5 years of CA tree mortality data using QGIS 3.6. The data can be found here. The downloads are geodatabase format. Navigate to each year's page and scroll down to the data section for the corresponding gdb. Based on email conversation with folks at US Forest Service, this could possibly be a mismatched geometry and projection file issue, but I don't know how to confirm that.

OTF projection is default in this version of QGIS and works well when adding years 2014-2016 (projected in EPSG:3310 - NAD83 / California Albers) to an OSM basemap. However, when adding years 2017 and 2018, QGIS asks for a projection because there seems to be no associated or an incorrect .prj (or the gdb equivalent). Despite setting these layers to the same EPSG 3310 but they end up near Alaska.

Things I have tried:

I emailed folks at the US Forest Service and they suggested trying EPSG:102003 - USA_Contiguous_Albers_Equal_Area_Conic. This brings the layer closer but still incorrect projection.
Exporting the layers as shapefiles with the CRS set to 3310 (and EPSG: 102003 CRS) and re-adding to QGIS. This does not work.
Disabling OTF projection, setting the project CRS to 3310 and manually setting each layer CRS - 17/18 still end up near Alaska.

Another contact at USFS was able to successfully overlay 17/18 with the prior years data using ArcGIS, so maybe there is something that the QGIS software is missing?

I don't know enough about the inner workings of the software.

Maybe there is some implicit change that can be made to the data to correct the projection issue?

Does anyone have other suggestions?


texts: Getting TMS tiles without clamping to Byte range

texts: How I can get TMS tiles (Z/X/Y.png) without getting them clamped in data range (e.g. GeoTIFF UInt16 -> Gdal2Tiles -> TMS UInt16 instead of GeoTIFF UInt16 -> Gdal2Tiles -> TMS Byte)?


texts:
How can I use OSRM to given a destination and 5 origins. Then OSRM should return the fastest origin, the distance between this origin and the destination and the duration of the route.

Can I do it with Table service? And how can I do it?


texts: Displaying QGIS plugin icons without resources.py file?

texts:
I am working on a new QGIS plugin and I am trying to comply with the requirement below. The full requirement is located here.

  source repository is in good shape:

  no generated file is left in repository (ui_*.py, resources_rc.py, gen. help files…)

I have a resources.py file.

If I remove resources.py file, no icon will be shown. So how can I achieve the visibility of resources without a resources.py file?


texts: Get the correct transform from bounding box with Rasterio

texts:
I have image as numpy array with no projection or geographical data and I have information of the bouding box of this numpy array (the image suppose to match perfectly the bbox).
The bbox is sentinelhub.geometry.BBox type.
I have the following informtion about the bbox:
bbox_size
>>>(850, 1279) #it has 1279 rows and 850 columns, not to be confused

bbox_coords_wgs84
>>>[-56.413932589, -12.93287647, -56.335907465, -12.81704823]  #lower left, upper right

The problem is that whenever  I'm trying to save my numpy array, I believe there is problem with the transform as I get it very disorted.
I'm using the following for the transform:
transform=rasterio.transform.from_bounds(*bbox, width= bbox_size[1], height= bbox_size[0])
transform
>>>Affine(6.1004788115715275e-05, 0.0, -56.413932589,
       0.0, -0.00013626851764705978, -12.81704823)

and the result look like this : (the red crosses are where there real lower left and upper right suppose to be:

Why is this happens? I believe there is something very basic that i'm missing here


texts: Can not cast class javax.measure.unit.Unit (tec.uom.se.unit.BaseUnit and javax.measure.unit.Unit)

texts:
I am having an issue with the below exception being thrown, not sure how to solve it.
Measure distance = .....
Unit<Length> unit = distance.getUnit();

CoordinateReferenceSystem origCRS  = CRS.decode("EPSG:3785")
unit = (Unit<Length>) origCRS.getCoordinateSystem().getAxis(0).getUnit();

exception:
class tec.uom.se.unit.BaseUnit cannot be cast to class javax.measure.unit.Unit (tec.uom.se.unit.BaseUnit and javax.measure.unit.Unit are in unnamed module of loader 'app')
    at com.vanitysoft.canada.boundaries.domain.dao.impl.CAPostalDAOImpl.createCircleRadis(CAPostalDAOImpl.java:131)

the return type of .getUnit() is javax.measure.Unit so I'm not understanding why this being thrown(I don't understand the stack trace).
I'm using GeoTools 22.4


texts: How to export attribute tables to Excel?

texts:
Possible Duplicate:
  How do I export an attribute table in ArcGIS version 9.3.1 to Excel?

We are running 10.0 and are trying to simply export an attribute table to Excel other than exporting to dbf and opening straight away in excel.

The problem is our GIS Manager does not trust everything that is exported as far as formatting etc..changed info from database to excel.

In other words..when we view the attribute table they want the same identical view in excel!


texts: Using GeoServer geowebcache watermark setting

texts:
I have successfully set watermark to my WMS layer, but in my OpenLayers 3 project, I use 'Tiled': true WMS setting and my GeoWebCache doesn't display any watermark.

I have same problem in the layer preview mode of GeoWebCache.

Does watermarking work only with "true" (=singletile) WMS?

I would like to use watermark for example for center of every rendered tile (or center of every fifth tile,..).

Is there an option to render watermark with GWC layer cache?


texts: gdal_grid - fails to find layer and field

texts:
I'm using gdal_grid for the first time and get this issue:

gdal_grid -zfield field_3 -l layer0 OLEX-LAT-190952014.vrt OLEX-LAT-190952014.tiff
ERROR 1: Failed to find layer 'layer0' on datasource 'OLEX-LAT-190952014.short.csv'.
Failed to find field field_3 on layer layer0, skipping.

Here is what my vrt file looks like:

<OGRVRTDataSource>
<OGRVRTLayer name="layer0">
<SrcDataSource>OLEX-LAT-190952014.short.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>EPSG:25832</LayerSRS>
<GeometryField separator=" " encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
</OGRVRTLayer>
</OGRVRTDataSource>

And my data file:

field_1,field_2,field_3
6932193.2570,377160.7010,-33.1300
6932194.6430,377162.4380,-32.6600

The first error mesage indicates that gdal_grid looks for a layer in the csv file - but how can a csv file contain a layer?
The next error message indicates that gdal_grid has read the csv file and can't find field_3 although it is clearly there?

Any ideas what I'm doing wrong?

GDAL version 1.10.1


texts:
I am trying to use Azure Functions as an API that returns a set of boxes(Lat Lng for each corner) along a route, final scope is to identify if a point is located near the route based on an offset distance.
For browser, I have found out that RouteBoxer is doing a great job at this however it can not work in Azure Functions Node.JS
There is also a package leaflet-headless however it uses some dependencies that are no longer available/are deprecated.
I am left with trying on Python. Would you recommend a library that achieves this?


texts: Using raster calculator inside graphical modeler

texts:
I have an issue, which I am struggling with since few hours - I've built a model which aims to make multicriterial analysis based on few vector layers and DEM.

At the step which is problematic for me I have 5 reclassified rasters and they are fine (I checked them), but after the whole process I get following error: "Error executing algorithm Raster calculator 'NoneType' object has no attribute 'crs' See log for more details". So I'm pretty sure, that my problem is calling those raster layers in not a proper way and unfortunately I cannot find any help in the Internet for this issue.

Below screens of the model.

There is also another issue - when I use Reclassify values (simple) from SAGA tools, when I prepare the lookup table and close a dialogue window it resets. Any ideas how to fight that?


texts: Which File to load from National Hydrography Dataset?

texts:
I'm using QGIS 2.0 on MacBook with OS X 10.8.5.
I downloaded a National Hydrography Database file. Unzipped, it is titled NHD11032_NM.gdb.
I try to add the layer as a vector file, and all the sub-files (nomenclature/terminology?) are ghosted, so I try adding the layer as a raster. All the files are visible and available.

Which of the dozens of files do I add to bring the file into my QGIS project?


texts: Creating flood area based on connected pixels of DEM in ArcGIS Desktop?

texts:
I am trying to show the area of flooding that occurs when bodies of water increase by different increments, like 0.5, 1, 1.5 meters, etc.. for all of Tanzania. I have a DEM with all of the bodies of water clipped out and the inverse -the dem of just the bodies of water. My thought is to locate all of the connected pixels to the body of water locations above 0.5, 1, and 1.5 meters (above the particular water source).

I tried experimenting with the cost distance tool in ArcGIS but am not sure if it is the correct tool for the job.

I prefer to use python and/or ArcGIS.

Is there a good way to show what areas would be flooded nearby bodies of water? Is the approach I am considering valid?


texts: Switch back to edit tool after saving edits?

texts: I'm using ArcMap 10.2.2 and after making some edits with the create feature tool, I would save those edits. The cursor would return from the cross hair icon back to the edit tool that lets you select features. Now after saving the cursor remains in create feature mode. Is there an option to switch it back to the edit tool cursor?


texts: QGIS Select by location QGIS 3.10.3 no precision parameter

texts: I have two layers, one polygon of output areas in the UK and a point layer of retail centers. I want to use select by location to find out how many output areas are within 1 km from a retail center. The select by location tool does not have the precision parameter to enter the 1km distance. Does anyone have any suggestion whether there's a plugin that does this? I've done this on ArcGIS but I want to know how it is processed in QGIS.


texts: How to iterate through various layers and export each with an additional word to its file name

texts:
I'm batch dissolving some layers in QGIS and want to save each dissolved shapefile as <layername>_dissolved.shp. I'm using a for-loop but it just shows raise QgsProcessingException(msg) _core.QgsProcessingException: There were errors executing the algorithm. error.
Here's my code:
import processing
layers = QgsProject.instance().mapLayers().values()
for i in layers:
    processing.run("native:dissolve",
                   {'INPUT': i,
                    'OUTPUT': r'C:\Users\imper\Desktop\\' + str(i) + '_dissolved' + '.shp'})
    iface.addVectorLayer(str(i+_dissolved), '', 'ogr')

How should I revise my 'OUTPUT' file path?


texts: Visualization of 3D shapefiles with Mapserver and Postgis

texts: I'm wondering is there a way to visualize 3d data with mapserver? 3D .shp files are in postgis database, and I already connected mapserver to database, but with simple layer in mapfile.


texts:
Can anyone suggest a good, Java-based raster library?  By good, I mean a library that can read common raster formats (ex. tif, img...etc) and save them to other image formats (ex. png).

I say Java-based because I am investigating a simple GIS app for Android.  I am aware that GDAL would seem like the obvious choice, but it seems like it will difficult to deploy GDAL with an Android application.  If anyone knows otherwise, please enlighten me.


texts: CartoDB tooltips don't appear when zoom value is not an integer

texts:
Building a custom CartoDB map here using the createLayer process. With the vast difference in visual appearance and scale between zoom levels of 4 and 5, I am aiming to default to a zoom level of 4.5. This works perfectly fine as far as the map view goes. However, when I am at this zoom level, none of the tooltips show up. They work at any other integer zoom level.

var options = {
    center: [13, -2],
    zoom: 4.5,
    minZoom: 4.5,
    maxZoom: 7
}

var layerSource = {
    user_name: '[username]',
    type: 'cartodb',
    sublayers: [
    {
        sql: "SELECT * FROM ptv_ad_views_by_dma_2015q4",
        cartocss: generateCSS("ads_0600"),
        interactivity: "cartodb_id,dma_name,ads_0600"
    },
    {
        sql: "SELECT * FROM ptv_ad_views_by_dma_2015q4",
        cartocss: generateCSS("ads_0700"),
        interactivity: "cartodb_id,dma_name,ads_0700"
    },
    ... (a dozen other sublayers)
}

var map_object = new L.Map("map", options);
var num_sublayers, map_layer;

cartodb.createLayer(map_object, layerSource)
    .addTo(map_object)
    .done(function(layer){

        num_sublayers = layer.getSubLayerCount();
        map_layer = layer;
        var sublayer;

        for (var i = 0; i < num_sublayers; i++) {

            sublayer = layer.getSubLayer(i);
            sublayer.setInteractivity('cartodb_id,dma_name,'+layerValues[i]);
            var testTooltip = layer.leafletMap.viz.addOverlay({
                type: 'tooltip',
                layer: sublayer,
                template: '

{{dma_name}}

{{'+layerValues[i]+'}} ads

',
                width: 200,
                position: 'bottom|right'
            });

        }

        hideAllSublayers();
        map_layer.getSubLayer(defaultLayer).show();
        populateTitle(defaultLayer);
        currentLayer = defaultLayer;

    })
    .error(function(err){
        console.log("An error occurred: " + err);
    });

Note that for some reason, the standard append line I've seen in support docs does not work and is not necessary for my map, in case you're wondering why I don't have it in there:

$('#map').append(testTooltip.render().el);

Any ideas of why these tooltips don't show up in non-integer zoom levels?


texts: Plotting high latitude bounding box using cartopy

texts:
I have used the sentinelhub python api to download a Sentinel-1 scene corresponding to the north coast of greenland where a bounding box defined by lon-lat is highly distorted. I have tried a few things to plot it but none seem to match the cartopy land mask that I'm using.
# Code to get lon/lat coords

# This is adapted from this guide
# https://forum.sentinel-hub.com/t/get-output-coordinates-using-sentinelhubpy/2804/4

bbox = BBox([-40, 80, -5, 84], crs=CRS.WGS84)
bbox_size = bbox_to_dimensions(bbox, resolution=500)

bb_utm = sentinelhub.geo_utils.to_utm_bbox(bbox)
transf = bb_utm.get_transform_vector(resx=500, resy=500)

pix_row = np.arange(0, bbox_size[1])
rows = np.array([pix_row] * bbox_size[0]).transpose()

pix_col = np.arange(0, bbox_size[0])
cols = np.array([pix_col] * bbox_size[1])

# Convert the pixel positions to UTM
utmx, utmy = geo_utils.pixel_to_utm(rows, cols, transf)

# Convert your UTM pixel positions to WGS84 (EPSG:4326)
lon_degrees, lat_degrees = geo_utils.to_wgs84(utmx, utmy, bb_utm.crs)

I have then tried several different ways of plotting this.
Method 1
North Polar Stereo with lon/lat coords
fig = plt.figure(figsize=(10,10))

ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=-30,true_scale_latitude=82))

# This is my bounding box
# [-40, 80, -5, 84]

ax.set_extent([-45, 0, 89, 79], ccrs.PlateCarree())

ax.pcolormesh(lon_degrees,
               lat_degrees,
               dtp[:-1,:-1],
                vmin = 0,
                vmax = 1,
                transform=ccrs.PlateCarree(),
                cmap='PRGn',
                alpha=1)

ax.scatter([-40,-5],[80,84],transform=ccrs.PlateCarree(),color='r',zorder=2)
ax.scatter([-40,-5],[84,80],transform=ccrs.PlateCarree(),color='r',zorder=2)

ax.add_feature(cartopy.feature.LAND, edgecolor='black',alpha=0.5,zorder=1)

The four corners of the bounding box are not rectangular (which is fine), but for some reason the four corners of the scene are.
Method 2
Mercator plot
Here I've replaced the plt.axes call with:
ax = plt.axes(projection=ccrs.Mercator())

Now the scene is non-rectangular, but the corners of the bounding box are rectangular

Method 3
Do it all in UTM
ax = plt.axes(projection=ccrs.UTM(28))

# This is my bounding box
# [-40, 80, -5, 84]

ax.set_extent([-45, 0, 89, 79], ccrs.PlateCarree())

ax.pcolormesh(utmx,
               utmy,
               dtp[:-1,:-1],
                vmin = 0,
                vmax = 1,
                transform=ccrs.UTM(28),
                cmap='PRGn',
                alpha=1)

Bounding box coordinates are now non-rectangular again, and scene is rectangular.
Bizarrely this is the closest fit between image and land mask, but the top right and bottom left corners of the bounding box no longer correspond to the corners of the scene!


texts:
none of these two unit tests runs fine but they both were supposed to.

public class CompareCrses extends TransformationTest {

  @Test
  public void comparingCrses() throws Exception {
    CoordinateReferenceSystem expectedCrs = CRS.decode("EPSG:32632");
    CoordinateReferenceSystem actualCrs = getEPSG32632ShapeValidFeatures().getSchema()
    .getCoordinateReferenceSystem();
    Assert.assertEquals(expectedCrs, actualCrs);
    // this one won't compile
//  Assert.assertTrue(expectedCrs.equals(actualCrs, false));
  }

  @Test
  public void comparingCrsesByTheirNames() throws Exception {
    CoordinateReferenceSystem expectedCrs = CRS.decode("EPSG:32632");
    CoordinateReferenceSystem actualCrs = getEPSG32632ShapeValidFeatures().getSchema()
    .getCoordinateReferenceSystem();
    Assert.assertEquals(expectedCrs.getName(), actualCrs.getName());
  }
}

I read how to compare crses but I did not found a solution. Variables differs but I don't get why. What am I missing? I am using geotools version 17.2. You can find the shape I use in the tests here.


texts:
I am writing my thesis about WMS and WFS services. I would like to compare these two services.

Can you suggest factors on which I can do this comparison?  (In my mind I have: functionality, fastest data transmission.)  Maybe you can also advise me about any good literature where I can find appropriate information?


texts: How do you change the value of a Boolean column in a QGIS 3 Attribute Table?

texts:
I am editing feature in a PostgreSQL table in QGIS 3 and while I can edit the values in most columns, I cannot edit the Boolean value.

I can enter the cell in edit mode, but keystrokes don't work.  I've tried to type true, 1, ctrl-mousewheel, double-clicking.

Aside from the obvious workaround of a SQL Update statement, is there a trick to changing a Boolean value within QGIS?

QGIS version
3.4.2-Madeira QGIS code revision 22034aa070 Compiled against Qt
5.11.2 Running against Qt
5.11.2 Compiled against GDAL/OGR
2.3.2 Running against GDAL/OGR
2.3.2 Compiled against GEOS
3.7.0-CAPI-1.11.0 Running against GEOS
3.7.0-CAPI-1.11.0 673b9939 PostgreSQL Client Version
9.2.4 SpatiaLite Version
4.3.0 QWT Version
6.1.3 QScintilla2 Version
2.10.8 Compiled against PROJ 520 Running against PROJ


texts: Is it possible to have run-time dynamic table updates to visualization using javascript interacting with postgis and geoserver

texts:
I have an idea of/have seen implementations, of dynamic map retrieval/updating using wms on openlayers by using updates to a postgis database which was hooked into geoserver and having a refresh policy. I have also seen examples done using cURL with php and python. Is this possible using JavaScript without using openlayers? The reason I don't want to use openlayers is because my "base map" is a google earth plugin and I've heard they don't really get along well.

In any case how would I update a postgis database which is connected to geoserver so that I can essentially move a point's coordinates and display this accordingly? And I know that using geowebcache is how geoserver updates the databases it is connected to, but I have heard issues regarding this before as such what was the process for making sure geoserver updates? Preferably I'd like to do this in JavaScript without some kind of wrapper.


texts: Determining how many cells are visible from Viewshed using QGIS?

texts: Using the Viewshed plugin for QGIS (2.16.3) how can I establish how many cells are visible or not from the output Raster?


texts: Multiply two lists of raster from two subfolder

texts:
I have two subfolders, which contains rasters for each country but which express different variables.  I am trying to figure out a way to Loop a multiplication of each raster into a folder with its pair of the other subfolder. Looks like:.

AFG_13      AFG_25.tif_int     for Afghanistan
AGO_13      AGO_25.tif_int     for Angola
etc...

where the first column rapresents the total population and the second is like a mask (infact I need only the pixel with values 25 converted in the pixel value 1)

In other words I need to multiply AFG_13*con(AFG_25.tif_int==25,1) in Python..I tried already to solve the problem in different ways, but still nothing...some advice?


texts: Pop up for features inside a GeoJson File

texts:
I uploaded a GeoJson file into GMaps API that has many polygon features in it, and now I'd like to create a pop up for each polygon. The content will be obtained from the attribute table of the GeoJson file from a field called "NAME".

Do you have a hint how to do that?