texts: Is product purchaser in distribution area?

texts:
I am going to be using PostGIS 2 to create a lookup service. Given a point (lat/lon), I need to either find if a product can be distributed there, or b) a list of products or distributors.

I am assuming I'd create a table with some distributor-related columns, some product-related columns, and a geographical column for stating a geographical distribution area for that distributor-product pair.

There will be hundreds of distributors and thousands of products.

What is the most peformant way to accomplish this? Should I use a geography datatype? Would geohashing be a better approach?
What would the queries look like?
If you have implemented something similar, what gotchas have you experienced?


texts: Location names autocomplete service with custom Nominatim installation

texts:
I'm looking to host an autocomplete service. I have installed Nominatim and imported the planet file with Nominatim. It took about 3 days to dump data into database and an additional 7 days to run indexes on an Amazon EC2 extra large instance without any database tuning.

The problem with Nominatim, is that it is designed to parse a complete address and not isolated letters. I would like to modify the code to accept what we could call partial queries, i. e. if the map is showing France, and the user types the letter 'p', the first result is "Paris, France".

https://github.com/twain47/Nominatim/blob/master/lib/Geocode.php#L343

I'm looking for advice from someone who has attempted this before in order to avoid diving into 1500 lines of code and hitting a snag.

Geocoding on my installation is not performing well even without type ahead, and I don't know if it will be possible to host a decent service.

I found this library https://github.com/komoot/photon which seems to work fast enough, but it doesn't seem to include a bounding box as a parameter to the search. I would like this service to work as good as possible with existing Open Source tools and without a third party hosting service.

Is typeahead achievable using available libraries and data, with a small effort, let's say 80 hours of coding/investigation?


texts: Programmatically change geometry column in OpenLayers WFS layer

texts:
Is there a way to programmatically change the geometry column on which a WFS layer is based and displayed?

WFS (from GeoServer) is correctly returning me two geometry columns, I want to change "on the fly" the geometry column OpenLayers uses to display items, obviously having OpenLayers displaying items (points, in this case) at the new coordinates.

I already tried to change the protocol.geometryName field, even forcing reload to refresh on the layer, but without success.


texts: Merging adjacent polygons using ArcGIS Desktop?

texts:
I want to use ArcGIS Desktop to merge adjacent polygons with each other.

There is no common field between the polygons and hence the Dissolve tool is not the ideal solution. The polygons that I want to merge together are either adjacent or are within a certain distance of each other.

I want to avoid using ArcObjects. However, Python scripts are welcome.

The Dissolve tool merges adjacent polygons when the attribute on which to dissolve is left out (and create multipart is disabled). However, it does not have settings for spatial tolerance. Furthermore, dissolve does not work on polygons that have just a vertex in common.

I have tried the approach with the integrate tool as mentioned in one of the answers. The problem is that the shape of the final feature that is getting created by integrating several feature is not like the shape that would have been formed if the boundaries between the features were dissolved.


texts:
I want to make a specific map in 1:1500 and 1:2000 scale on A2 paper but the label height should be 2.5 mm in composer.

Is there a correlations between label size in "map units" or "points" to produce that specific size in composer.

All this time I just doing Trial and Error method.


texts: Exporting Raster Attribute Table using Open Source software

texts:
I have 120 or so geotiff files that I'd like to build and export the raster attribute table to something a little more friendly (e.g. csv).

I'm curious if there is an open source way to do this. The files are rather large 36001 X 36001 pixels so brute force solutions may be out. Ideally, I'm looking to do this in gdal commandline/OSGEO4w or Qgis. Grass in a pinch.

Some additions in response to comments:

I am somewhat familiar with python. If anyone has a solution in that vein, I'd love to hear it. Same thing goes with R and windows command line.

-The table would look something like this:

  Value | # of Cells
  0     | 100000
  1     | 3214
  2     | 25125
  ...
  98    | 2214213

Basically, I'm looking for a non-ArcGIS way to recreate the ArcGIS raster attribute table. The reason I want CSV is because I have some other mathematical operations I'd like to run, and a csv or some sort of text file works best for my workflow. Its a bit anti-dramatic, but all I really need at this stage is the number of cells for each value in the raster.


texts: ArcGIS Online: embedded map - zoom to feature

texts:
I created a webmap for embedded use: CLICK.
I want to zoom to a certain point with url parameter.

I followed this for the settings: ESRI documentation. Towards the bottom of the page it says following:

  feature—Allows searching for a specific feature in a searchable feature layer. Find Locations by Layer needs to be set up on the map for a searchable field. The syntax is feature=layerID;searchfield;searchvalue, for example, feature=CentralIndianaCities_634;PLACEFIPS;05860.
  The following URL references the central Indiana cities example above and can be used to discover the layer ID of the desired layer in your map. Notice id:CentralIndianaCities_634 under operationalLayers. Replace the webmap ID with the ID of your web map to find your layer ID.

I set up the searchable field in the settings (sorry, it's in Dutch). It just says, in the layer "speeltoestellen", the field "topdesk" must contain the parameter

but for some reason the map doesn't zoom to the selected point. What am I missing to get this to work?


texts: QGIS "freeze(?)" on Windows 10 when using color sampling tool?

texts:
On the Windows 10 computers in our internal IT infrasctructure QGIS freezes when trying to use color sample tool. We use 2.18.21 and I never observed this behaviour before on clients outside our infrastructure. No color is sampled when clicking somewhere with the pipette cursor and the whole windows becomes almost inoperable, apart from killing qgis task in the task manager.

I just observed the following behaviour:

I Select 'sample color'. Pipette cursor apears and can be moved around. But clicking somwhere does not chose a color to by sampled nor focus is given back to 'Select color' dialogue. When doing a Ctrl-Alt-Del, locking the workstation and sign in to the workstation again, the color is selected and focus given back to 'Select color' dialogue.

I just installed QGIS 2.18.24 LTR and QGIS 3.2.1 in c:\tmp with OSGeo4w installer and observed the same behaviour with 2.18.24 but not with 3.2.1, where color sampling works as expected.

What might cause this behaviour?


texts: Snapping Polygon vertices to Point features using ArcGIS?

texts:
I have a set of polygons drawn using an inaccurate GPS PDA. these polygons are sharing nodes.

After post-processing GPS data to get accurate coordinates i have a set of all the nodes, then i have to update the vertices of already drawn polygons. the relation between old and new vertices is based on the nearest point.

How could i modify polygon vertices in this case, using ArcGIS Engine , ArcObjects or just ArcGIS Desktop?


texts:
Main question: How do I manually edit raster layers on QGIS 3.0?

Some context: I georeferenced Google Earth high resolution downloads on to QGIS, and was looking to do image classification for a study on urban areas. Everything I read suggested using ORFEO Toolbox, but my QGIS 2.18, which was supposed to come pre-equipped with ORFEO never seemed to work for me and QGIS 3.0 is incompatible with ORFEO. Then I tried the Semi-automatic Classification Plug-in (SCP) that is also recommended, but it doesn't seem to like my GE screen grabs because they only have 3 bands (everything I read says it should be 4 or more bands). Then I thought I might try just editing the pixels manually so that one color meant built-up, another vegetation, etc., which would take forever, but at least would be more straight forward- but I can't figure out how to edit pixels on QGIS (I read that the Serval Plugin would work, but it is incompatible with QGIS 3.0 and before going back to QGIS 2.18, I wanted to see if there were any other recommendations). Ultimately, I'd like to create a vector layer with an attribute table that has the %age of land-type per polygon (land-type percentages coming from the raster that I create from the georeferenced GE image).

I am using a MacBook Air.


texts:
I want to use default icon provided by leaflet instead of using something else. But unfortunately i am not seeing anything, not even error.

var myIcon = L.icon({
            iconUrl: 'my-icon.png',
            iconSize: [38, 95],
            iconAnchor: [22, 94],
            popupAnchor: [-3, -76],
            shadowUrl: 'my-icon-shadow.png',
            shadowSize: [68, 95],
            shadowAnchor: [22, 94]
        });
        L.marker([50,10], {icon: myIcon}).addTo(this.map);


texts: QGIS - Lookup elevation or slope from geoTIFF from a series of lat/lon values

texts:
I have a geoTIFF file (pretty large) covering a good chunk of Spain around Almeria. I got it from the EU-DEM v1.1 data (25m resolution), freely available in which I have a route in lat/long and I want to look up the elevation values at each point or the surface-slope between points, but elevation might be easier.

So far,

loaded the geoTiff in.
figured out how to add a new layer that can be pointed.
installed the point sampling tool and I found a GDAL function called slope that made a slope layer.

But I have no idea how I get QGIS to load in my sample lat/long points? and if that would even work as I see the coordinate system in the bar across the bottom seems to show pixel X, Y.


texts: How to get street information from latitude and longitude

texts: How can I get street details from particular longitude and latitude? Is that particular street highway/city street ? How can I get this particular information?


texts:
New to QGIS, I'm facing the Problem that my lat, lon information has seemed to be lost creating polygons.

I created layers importing csv files being gps points. From these points I used the function points to line and then line to polygon which I saved as shapefile. It seems that now in my polygon there is no information about the latitude longitude anymore.

What can I do to retrieve the Information or how can I keep it during the process?


texts: Disabling access to ArcGIS Server Manager externally?

texts:
The background of my server

I installed the ArcGIS server 10.3.1 (No Web Adaptor installed)

I have Windows 2012 Server x64 with SQL Server 2012 R2.

I'm also using IIS 8 as a web server on the same machine.

I am wondering if I could disable the external access to the link for the ArcGIS Server Manager.

For example, assuming the link is http://gisserver.domain.com:6080/arcgis/manager/. Could I know if there is a way to disable/block the external users to access to the link?


texts: How to debug ogr2ogr ERROR 1: Did not get at least 2 values or invalid number of set of coordinates?

texts:
I'm trying to convert a government data contour GML (http://geo.data.linz.gv.at/katalog/geodata/hoehenschichten/2011/1m-Hoehenschichten.gml) to Shapefile.

It seems like it was created by FME:

<?xml version="1.0" encoding="UTF-8"?>
<gml:FeatureCollection xmlns:gml="http://www.opengis.net/gml" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:fme="http://www.safe.com/gml/fme" xsi:schemaLocation="http://www.safe.com/gml/fme 1m-Contour.xsd">
<gml:boundedBy>
<gml:Envelope srsName="AUT-GK31-5" srsDimension="3">
<gml:lowerCorner>67741.3127211605 341935.72972973 245</gml:lowerCorner>
<gml:upperCorner>79959.5666279546 360502.420825918 602</gml:upperCorner>
</gml:Envelope>
</gml:boundedBy>
<gml:featureMember>
<fme:Contour gml:id="ida7f8be19-1a86-42be-886c-2238b778e00d">
<fme:_elevation>434</fme:_elevation>
<fme:_length>2987.20354766336</fme:_length>
<gml:curveProperty>
<gml:LineString srsName="AUT-GK31-5" srsDimension="3">
<gml:posList>68485.5237522595 355869.11751307 434 68487 355867.757575758 434 68488.6860465116 355866 434 68492 355862.369426752 434 68493.6929133858 355861 434 68497 355858.324840764 434 68498.921052631$
</gml:LineString>
</gml:curveProperty>
</fme:Contour>
</gml:featureMember>

Ogrinfo (OGR 1.9) recognizes multi line strings:

> ogrinfo 1m-Hoehenschichten.gml
Had to open data source read-only.
INFO: Open of `1m-Hoehenschichten.gml'
  using driver `GML' successful.
1: Contour (Multi Line String)

But

 ogr2ogr -f "ESRI Shapefile" -nlt MULTILINESTRING25D  -skipfailures  out 1m-Hoehenschichten.gml Contour

fails with

ERROR 1: Did not get at least 2 values or invalid number of
set of coordinates <gml:posList>68485.5237522595 355869.11751307 434 68487 355867.757575758 434 68488.6860465116 355866 434 68492 355862.369426752 434 68493.6929133858 355861 434 68497 355858.324840764 434 68498.9210526316 355856 434 68502 355851.9375 434 68502.5818965517 355851 434 68503.5340909091 355849.465909091 434 68505.6853448276 355846 434 68507 355843.896551724 434 68509.2459893048 355841 434 68512 355836.846774194 434 ....

My interpretation so far is that OGR does not recognize that the file contains 3D geometries.

How can I find out if the file is broken or whether I'm dealing with an OGR bug/missing feature?


texts: postGIS, length, area, density of line network

texts:
I am trying to figure out the length and density of all major roads and minor roads within a study area to return something like:

|STUDYAREA|a_road_length|a_road_density|other_road_length|other_road_density

There are two polygons within the study area - an inner focus and an outer buffer.  I want the length and density grouped by these two areas.
Major roads are only one type (a_roads) but other roads are stored in two different relations

|a_road|name|geom

|b_road|name|geom

|minor_road|name|geom

So far I can get the length of each separately using ST_Length, the area of the polygons using ST_Area and the density of a_road but not all the elements for all the roads together in one result.

Perhaps need to use subquery or

create VIEW

but unsure how to put all these elements together.

In the example above red line length and density for inner and outer area then combine green and yellow length and density for inner and outer area to return as "other roads".


texts:
I am beginner to map plotting. I am currently trying to create some custom map-plots on historical Europe on the datawrapper website.
I have obtained a GeoJSON file from this link (https://github.com/aourednik/historical-basemaps/blob/master/world_1920.geojson)
However, its for the whole world and I need only for Europe in 1920s.
So, is there some way I can generate a custom GeoJSON
(I am inexperienced in tools like ArcGIS)


texts: How do I create/find .cpg files after creating a shapefile in QGIS?

texts:
I've been asked to map an area of land using QGIS, so I created a new vector layer polygon and saved it to my PC without problems.

They've asked for the following files: .shp, .shx, .dbf, .prj and .cpg.

When I go to the folder containing my saved polygon files, all files except .cpg are there.

How do I create a .cpg file for my polygon?


texts: How to convert curved features into geojson

texts:
I'm using a MapGuide system which provides some CURVEPOLYGON(s).

Using python shapely's loads() raises a WKTReadingError.

Tried to find a tool that can approximate curved shapes into a polygon (as a set of points) but no luck.

I suspect ogr2ogr to be able to do this (not sure if I'm stupid) but wasn't able to do so...

Has anyone been successful when dealing with such shapes?


texts: How to join linestrings without reversing directions in PostGIS

texts:
I have a collection of linestrings, which represent one-way streets. I linemerge these together when they have matching attributes, to get longer linestrings for better rendering. But I want to avoid reversing the direction of any of the linestrings during the join, since the direction of one-way streets is important.

In this first example, there are three streets which can be merged together, since they are all 'nose to tail':

select st_astext(st_linemerge(st_geomfromtext(
  'MULTILINESTRING((0 0, 1 1), (1 1, 2 2), (2 2, 3 3))'
)));
          st_astext
-----------------------------
 LINESTRING(0 0,1 1,2 2,3 3)
(1 row)

This is fine. But in the second example, the middle linestring points in the opposite direction:

select st_astext(st_linemerge(st_geomfromtext(
  'MULTILINESTRING((0 0, 1 1), (2 2, 1 1), (2 2, 3 3))'
)));
          st_astext
-----------------------------
 LINESTRING(0 0,1 1,2 2,3 3)

The direction of the middle linestring has been reversed during the merge, from (2 2 -> 1 1) to (1 1 -> 2 2).

Is there a method to merge linestrings in PostGIS that is similar to st_linemerge, but that avoids reversing the direction of linestrings?


texts:
In the context of QGIS/PyQGIS, I would like to hide the unselected features after performing, for example, a simple selection.

feature = layer.getFeatures().next()
layer.setSelectedFeatures([feature.id()])

After running this commands a feature becomes highlighted. Besides, I'd like to hide the unselected features. How can I do that?

I believe this is possible because the Query Builder performs a filter that hides the undesired features.


texts: How to use QgsCoordinateTransformContext to export vector layer to another CRS?

texts:
I am having troubles to export a vector layer to another Coordinate Reference System (CRS) using WriteAsVectorFormatV2 and QgsCoordinateTransformContext.
I looked in many posts related to QgsCoordinateTransformContext but did not find any example on how to use it for my needs.

Here is what I am trying to do:

##########
# Test 1 #
##########
# This is a vector layer of point geometry and CRS set to EPSG:4326
input = iface.activeLayer()

context = QgsCoordinateTransformContext()
ref_crs = QgsCoordinateReferenceSystem("EPSG:4326")
dest_crs = QgsCoordinateReferenceSystem("EPSG:4258")
# According to the documentation, if the third element is void, it should use the default PROJ conversion string
context.addCoordinateOperation(ref_crs, dest_crs, "")

s_out = r"test_out_1.gpkg"
o_save_options = QgsVectorFileWriter.SaveVectorOptions()
o_save_options.layerName = lyr_airports.name()
tpl_err = QgsVectorFileWriter.writeAsVectorFormatV2(lyr_airports, s_out, context, o_save_options)
# The layer is written and I can load it in QGIS but it does not have information about CRS

##########
# Test 2 #
##########
# This is a vector layer of point geometry and CRS set to EPSG:4326
input = iface.activeLayer()

context = QgsCoordinateTransformContext()
ref_crs = QgsCoordinateReferenceSystem("EPSG:4326")
dest_crs = QgsCoordinateReferenceSystem("EPSG:4258")
# I obtained the PROJ string with the following command: projinfo -o PROJ -s EPSG:4326 -t EPSG:4258
context.addCoordinateOperation(ref_crs, dest_crs, "+proj=noop")

s_out = r"test_out_2.gpkg"
o_save_options = QgsVectorFileWriter.SaveVectorOptions()
o_save_options.layerName = lyr_airports.name()
tpl_err = QgsVectorFileWriter.writeAsVectorFormatV2(lyr_airports, s_out, context, o_save_options)
# The layer is written and I can load it in QGIS but it does not have information about CRS

In both test cases, I get the output layer and I can load it back in QGIS, where it is displayed without any information about the CRS.

Can anyone advise on what I am doing wrong? Thanks in advance for your help.


texts: The Info Windows content change for Google Maps Engine cannot be adopted

texts:
When I edit the Info Window, I modified

 View Details

as

 View Details

However, in the map preview, the new attribute target="_blank is not added into the Info Window.

Did I miss anything?


texts: How can i set layers name from map to combobox?

texts: I'd like to list map layers from map into combobox. I tried to set itemsource of combobox to: ItemsSource = "{Binding ElementName=Map, Path=Layers}. It does not list the layer names. It just lists "Esri.ArcGIS.Client.FeatureLayer". How can I list layers with their names in my combobox?


texts: Adding a FeatureLayer ARCGIS JavaScript API

texts:
I am currently exploring ArcGIS JavaScript API. I have my own map services and feature services up on my local server. I would want to use the FeatuerLayer to select and highlight feature on the map. But when I add my feature service I get the following error.

TypeError: Unable to draw graphic (geometry:null, symbol:null): _14 is undefined

I am not able to figure out what this could be, and I am not sure if the feature layer has been added to the map or not. How could I check if it has been added? I tried using a dojo.connect method with the FeatureLayer with an onClick event. But this doesn't seem to work. Have placed the code below

Code:

esri.config.defaults.io.proxyUrl = "proxy.ashx";
dojo.require("dijit.layout.BorderContainer");
dojo.require("dijit.layout.ContentPane");
dojo.require("esri.map");
dojo.require("esri.layers.FeatureLayer");
dojo.require("esri.tasks.query");

var map;
var selectionToolbar;
var infoTemplate = new esri.InfoTemplate();
infoTemplate.setTitle("${ROADNAME}");
infoTemplate.setContent("ROAD NAME: ${LINK_ID}
" + "CAT${CAT}");

function init() {
  try {
    var initExtent = new esri.geometry.Extent({
      "xmin": 103.55,
      "ymin": 1.13,
      "xmax": 104.16,
      "ymax": 1.56,
      "spatialReference": {
        "wkid": 4326
      }
    });

    map = new esri.Map("map", {
      extent: esri.geometry.geographicToWebMercator(initExtent)
    });

    var basemap = new esri.layers.ArcGISTiledMapServiceLayer("http://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer");

    //my map service layer
    var basemap_sing = new esri.layers.ArcGISTiledMapServiceLayer("http://karthikpc:8399/arcgis/rest/services/Carriage_Mercantor/MapServer", {
      displayLevels: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
    });
    // my feature layer which i am not able to use
    var featureLayer = new esri.layers.FeatureLayer("http://karthikpc:8399/arcgis/rest/services/Carriage_Mercantor/FeatureServer", {
      mode: esri.layers.FeatureLayer.MODE_ONDEMAND,
      infoTemplate: infoTemplate,
      outFields: ["*"]
    });

    map.addLayer(basemap);
    map.addLayer(base_sing);
    map.addLayer(featureLayer);

    dojo.connect(map, "onLoad", initSelectToolbar);
    dojo.connect(featureLayer."onClick", helloworld);
  } catch (e) {
    alert('An error has occurred: ' + e.message);
  }

}

function helloworld() {
  alert("hello world");
}

dojo.addOnLoad(init);


texts:
I have a SpatialLinesData frame with 1657 different routes.  Each route is of a different length and has a differing number of nodes or waypoints (lat, lon).

For example, routes@line[[87]] is made up of 7 waypoints.  I am wondering if there is a simple method to extract the lat, lon coordinates, or [,1] and[,2] of every route and convert them into a list of 1657.

When I try this with one routes@line[[87]], for example, I get an error:

    > as.list(routes@lines[[87]])
    Error in as.list.default(routes@lines[[87]]) :
      no method for coercing this S4 class to a vector

Does a method exist for converting line objects to lists?

    > str(routes)
    Formal class 'SpatialLinesDataFrame' [package "sp"] with 4 slots
      ..@ data       :'data.frame': 1657 obs. of  3 variables:
      .. ..$ start_time: Factor w/ 1631 levels "2016/09/09 00:00:02",..: 1 2 3 4 5 6 7 8 9 10 ...
      .. ..$ duration  : int [1:1657] 786 1248 738 786 651 660 616 889 408 475 ...
      .. ..$ difftime  :Class 'difftime'  atomic [1:1657] 2 4 19 67 92 119 170 202 206 213 ...
      .. .. .. ..- attr(*, "units")= chr "secs"
      ..@ lines      :List of 1657

    > routes@lines[[87]]
    An object of class "Lines"
    Slot "Lines":
    [[1]]
    An object of class "Line"
    Slot "coords":
              [,1]     [,2]
    [1,] -79.92976 43.25867
    [2,] -79.92984 43.25865
    [3,] -79.92986 43.25868
    [4,] -79.92989 43.25870
    [5,] -79.92983 43.25871
    [6,] -79.92983 43.25871
    [7,] -79.92983 43.25871

    Slot "ID":
    [1] "87"


texts: QGIS: split polygon features with different polygon layer

texts:
I have two polygon layers. One containing vegetation features and one containing the borders of different monitoring areas. I want to calculate the total size of different vegetation types per monitoring area, however, the vegetation features often go over the borders of two areas (see example below, with red monitoring area border).

I want to split the vegetation features that cross borders so that they fit to the borders of the monitoring areas. How can that be done? I found the split line by line tool in QGIS but couldn´t find the same for polygon. Splitting manually is not an option as I have many thousand features.


texts: What are the units of the area of DIVA-GIS shapefiles

texts:
I've been checking several other post and even sent messages to DIVA-GIS, but still cannot figure out the answer by myself.

I want to know the area of a shapefile downloaded from DIVA-GIS: http://www.diva-gis.org/gdata and sowmehow figure out the units of it.

So for instance, I got the shapefile for the US and extracted only the state of Alabama using R:

usa.shp <- readShapePoly("countries/USA_adm1.shp")
ala <- usa.shp[usa.shp$NAME_1=="Alabama",]

Then I estimate the area using gArea() function from rgeos:

gArea(ala)
[1] 12.867

And this is also the information provided to me in the shapefile:

ala$Shape_Area
[1] 12.867

But when I check about the area of Alabama in Wikipedia this is the information I get:

52,419 sq mi
135,765 km2

I haven't worked with geospatial data before and I would want to stranform this area to km2 if that's possible.


texts:
Is there any guideline or procedure to integrate Solr in GeoNetwork for metadata searching?
I have tried search in GeoNetwork website, it only mention Apache lucene.


texts:
I am programming with ArcObjects and trying to output a raster to a folder path that contains spaces. Before creating the raster I am trying to validate that the destination path is valid. The following call to IsWorkspace() returns false when the path has spaces. Is this by design or a problem? How should one test that a folder is a valid file-based raster workspace?

Dim aType As Type = Type.GetTypeFromProgID("esriDataSourcesRaster.RasterWorkspaceFactory")
Dim obj As System.Object = Activator.CreateInstance(aType)
Dim pWSF As ESRI.ArcGIS.Geodatabase.IWorkspaceFactory = obj
return pWSF.IsWorkspace("C:\MyGISData\My Output Folder")


texts: Features with geomety issues or unexpected CRS in Geopackage tables

texts:
In my project I've been asked to make some data validations to GeoPackage tables previous to start processing the data. I'm doing this with FME Dekstop.
Among other properties I should validate that features in each table have the expected coordinate system and geometry type. I have some experience working with shapefiles and I've come across features with null geometry and other geometry problems or a different coordinate system than the rest of features, but I have few experience working with GeoPackage and don't know whether this same issues may occur with GeoPackage tables.
If that's not possible it would simplify the validation and reporting/logging process. Assuming all features have same properties I would just check the defined geometry type and coordinate system for the table.
Using QGIS I've been trying to create features in a GeoPackage table with different properties that the rest to no avail and I can't find an example or issue anyone may have had with that. So I need your experience to reply to this question and/or give me some tips to intentionally generating incorrect features in a GeoPackage table using QGIS.


texts: Configurable Map Viewer (CMV) - average the values of the points

texts:
I have been able to select points and find the sum of the points. What I really need is to average the points.
 Here is my (acquired from someone else) code for finding the sum.

define([
'dojo/_base/declare',
'dijit/_WidgetBase',
'dijit/_TemplatedMixin',
'dijit/_WidgetsInTemplateMixin',
'esri/toolbars/draw',
'esri/InfoTemplate',
'esri/layers/FeatureLayer',
'esri/symbols/SimpleMarkerSymbol',
'esri/symbols/SimpleLineSymbol',
'esri/Color',
'esri/tasks/query',
'dijit/form/Button',
'dojo/_base/lang',
'dojo/on',
'dojo/_base/array',
'dojo/dom',
'dojo/text!./Select/templates/Select.html',
'xstyle/css!./Select/css/Select.css'
], function (declare, _WidgetBase, _TemplatedMixin, _WidgetsInTemplateMixin, Draw, InfoTemplate, FeatureLayer, SimpleMarkerSymbol, SimpleLineSymbol, Color, Query, Button, lang, on, arrayUtil, dom, SelectTemplate, css) {
    return declare([_WidgetBase, _TemplatedMixin, _WidgetsInTemplateMixin], {
        widgetsInTemplate: true,
        templateString: SelectTemplate,
        selectTools: null,

        postCreate: function () {
            this.selectTools = new Draw(this.map);

            this.selectQuery = new Query();

            on(this.selectTools, "DrawEnd", lang.hitch(this, function (geometry) {
                this.selectTools.deactivate();
                this.selectQuery.geometry = geometry;

                this.featureLayer.selectFeatures(this.selectQuery,
                FeatureLayer.SELECTION_NEW);
            }));
this.fieldsSelectionSymbol =
             new SimpleMarkerSymbol(SimpleMarkerSymbol.STYLE_CIRCLE, 10,
new SimpleLineSymbol(SimpleLineSymbol.STYLE_SOLID,
new Color([255,0,0]), 1),
new Color([0,255,0,0.25]));

            var infoTemplate = new InfoTemplate(this.title, this.content);

            this.featureLayer = new FeatureLayer(this.featureLayerURL,
            {
                mode: FeatureLayer.MODE_ONDEMAND,
                infoTemplate: infoTemplate,
                outFields: ["*"]
            });

            this.featureLayer.setDefinitionExpression(this.defExpress);
            this.featureLayer.setSelectionSymbol(this.fieldsSelectionSymbol);
            this.featureLayer.on("selection-complete", lang.hitch(this, 'sumGasProduction'));
            var clearMsg = this.clearSelectMsg;
            this.featureLayer.on("selection-clear", function () {
                dom.byId('messages').innerHTML = "" + clearMsg + "";
            });

            this.map.addLayer(this.featureLayer);

        },

        select: function () {
            this.selectTools.activate(Draw.EXTENT);
        },
        clearSelection: function () {
            this.featureLayer.clearSelection();
        },

        sumGasProduction: function (event) {
                        var productionSum = 0;
                        var sumField = this.fieldToSum;

                        arrayUtil.forEach(event.features, function (feature) {
                            productionSum += feature.attributes[sumField];

                        });
                        dom.byId('messages').innerHTML = "" + this.selectSumMsgPrefix + " " + productionSum + " " + this.selectSumUnits + ". ";
                    }
    });
});

Secondly my clear selection button is not working.

Here is my code for clear selection:

    <button data-dojo-type="dijit/form/Button"  data-dojo-attach-event="onClick:clearSelection" >
        Clear your Selection
    </button>

I have no idea how to get the average. I was trying to use what you said in another answer but I could not figure it out. I have found an example that does find the average

https://developers.arcgis.com/javascript/jssamples/query_bypoly.html

      function average(fset, att) {
    var features = fset.features;
    var sum = 0;
    var featuresLength = features.length;
    for (var x = 0; x < featuresLength; x++) {
      sum = sum + features[x].attributes[att];
    }
    return Math.round(sum/featuresLength);
  }


texts:
I created three functions the first open the raster file and test if it's correct or not. the second calculate a parameter and the third create my output raster file. the issue is that the output file is not being created. I'm stuck at this level please help me if you can.

    # -*- coding: cp1252 -*-
from osgeo import gdal
import numpy as np
import sys ,os
from osgeo.gdalconst import *

output_dir='C:\wamp\www\Données CRTS\irrisat_dessimination-smiej'
etp = os.path.join(output_dir,"L8_L8_ETpot_24_30m_2016_223.tif")
rain = os.path.join(output_dir,"GLDAS_NOAH025SUBP_3H.A2016223.sumday.rain.tif")
eff = os.path.join(output_dir,"eff_v11.tif")

class calculate_iwr :
    path=''
    def openRaster(self,raster):
        raster = gdal.Open(path)
        band=None
        if raster is None :
            print ("Erreur : Impossible d'ouvrir le raster: " )
        try :
            columns = raster.RasterXSize
            lines = raster.RasterYSize
            band = raster.GetRasterBand(1).ReadAsArray(0,0,columns,lines)
        except:
            print ("Erreur : Impossible d'extraire la bande")
        return raster, band, columns, lines

    def parameters(self,parameter) :
        bands = []
        bands.append(self.band)
        if parameter == "IWR" :
            iwr=np.divide(np.subtract(bands[0]-bands[1]),bands[2])
        return iwr

    driver=gdal.GetDriverByName('GTIff')
    driver.Register()

    #std=np.std(NDVI)
    #mean=np.mean(NDVI)
    #new_image=driver.Create(dst_f, colonnes,lignes, 1, gdal.GDT_Byte)
    #for canal in cannaux :
        #i=i+1
        #new_image.GetRasterBand(i).WriteArray(canal)
    #new_image.GetRasterBand(1).WriteArray(codeim)
    #new_image=None
    #histogramme (NDVI,std,mean)

    def createImage(new_image):

        new_image=driver.Create(path, self.columns,self.lines, 1, gdal.GDT_Byte)
        new_image.SetProjection(projectionfrom)
        new_image.SetGeoTransform(geotransform)
        new_image.GetRasterBand(1).WriteRaster( 0, 0, columns, lines, output_dir )
        new_image=None
        return columns, lines

def main() :

    openRaster(etp)
    openRaster(rain)
    openRaster(eff)
    parameters("IWR")
    createImage(os.path.join(output_dir,"iwr.tif"))


texts: Understanding LSOA and MSOA boundaries (UK)?

texts:
I'm currently attempting to do a statistics project summarising data across UK administrative boundaries, in particular Lower Layer Super Output Areas (LSOAs) and Middle Layer Super Output Areas (MSOAs) as defined by the Office of National Statistics (ONS).

However, I'm struggling to understand how LSOAs and MSOAs relate to each other.

Looking at the boundaries map available at neighbourhood statistics gives the impression that LSOAs are a breakdown of MSOAs however I can't find any documentation confirming this or a list mapping LSOAs to MSOAs.

I was hoping someone might be able to help explain how they relate to each other or point me in the way of some relevant documentation.

For reference I managed to find a basic description here, but it doesn't fully explain the relationship.


texts: How can I export a list of an mxd's content (layers names) into a table or a pdf?

texts:
I am using Arcgis 10.2.

I have created a set of mxd from which I need a list of the layers and their characteristics (Name, sources, symbologie - maybe as in legend, if there is an Definiton Query, and other representational aspects). I would like this list to be in a table or a Report or something simple to read/understand.

I need this in order to inform the colleagues what data they have to work, how and why in order to simplify the integration in the mxd templates.

This seems like a Report Job but I couldn't figure out why it was not doing what I needed. This sounds such a Basic Task that there must be a tool somewhere?

Have any of you seen some tools doing that type of export?

I hope that made sense. Will definitely clarify as much as needed.

Many thanks

EDIT: See Alex' reply


texts:
I am trying to run Shapely on a new machine, and tried both installing it from source, using the installer on its webpage, and using pip install shapely. Every reference I try to make from the library except from shapely import * throws the error WindowsError: [Error 126] The specified module could not be found. Has anyone had a similar experience when trying to leverage Shapely, and figure out how to get it to function correctly?

Python Error Thrown:
WindowsError: [Error 126] The specified module could not be found

Environment:
Windows 7 64-bit
Python 2.7


texts: How to create DEM from Lidar point cloud data?

texts:
I'm using ArcMap and I want to create DEM(s) over an area to have different outputs depending on which settings/filters I put in. I have access to laser point cloud data as my input files.
I've worked with the LAS dataset tools but I'm not really sure if I'm doing everything correct or not. Because the points aren't regular, maybe some interpolation method is needed first? The problem with that is, I don't have any input point features.

This link describes fairly good the steps I've made: http://www.edc.uri.edu/blog/using-las-datasets-create-functional-outputs


texts: Iterating through pair of points to create lines

texts:
I have a table with ~12000 rows like:

NOMBRE     LAT     LONG      FECHA                INT
AL011900    15  -42.1   1900-08-27T00:00:00Z    18.0054
AL011900    15.2    -43.4   1900-08-27T06:00:00Z    18.0054
AL011900    15.3    -44.7   1900-08-27T12:00:00Z    18.0054
AL011900    15.4    -45.6   1900-08-27T18:00:00Z    18.0054
AL021900    19  -59.3   1900-09-13T12:00:00Z    33.4386
AL021900    19.5    -60 1900-09-13T18:00:00Z    36.0108
AL021900    20  -60.6   1900-09-14T00:00:00Z    38.583
AL041905    36.3    -48.6   1905-10-10T18:00:00Z    46.2996
AL041905    37.9    -47.9   1905-10-11T00:00:00Z    43.7274
AL041905    39.6    -47.1   1905-10-11T06:00:00Z    41.1552
AL041905    41  -46 1905-10-11T12:00:00Z    41.1552

Fecha is "date time" and the table is already sorted using it.

I will like to create lines between sequential points per each NOMBRE, keeping in the attribute of each line the INT values of the first point of each line. Please note that there are several INT per NOMBRE and that's the reason why I need to create lines for each consecutive pair of points.

That will be: the line from point 1 to point 2 will have the INT attribute of point 1. The line from point 2 to point 3 will have the INT value of point 2. And so on for each of the NOMBRE.


texts: QGIS 3 on Mac OS - Where is the configuration directory?

texts:
I've got a fresh clean install of QGIS 3 running on Mac OS and I can't find where the plugin directory should be.

I don't have a directory at ~/.qgis3.

I've search my drive with locate qgis and locate plugins | grep qgis and can't see it anywhere.

I actually have two QGIS 3 installs, one using Homebrew and the other using Kyngchaos, and when I start either, they list the same installed plugins, so there must be a directory somewhere!

So, where is the QGIS 3 config directory on Mac OS hiding?

Edit: Nyall helped me find it! It's at ~/Library/Application\ Support/QGIS/QGIS3/profiles/default/ (the default component is a named profile, which is a new feature).


texts: postgis st_relate polygon inside another and sharing a boundary

texts:
I am trying to select polygons inside an other one WITH a common boundary.
It seems that st_relate is the way but I can't find the right matrix, even if I don't see where the fault is in :
st_relate(e.geom, d.geom, '2FF11F212')

[For example I want to select the bottom green polygon only


texts: I have a NDVI that I created from a Landsat 5 image. Is there a way to extract all the pixels above a certain threshold and save it as a shapefile? I am trying to create a vegetation mask of sorts to filter out areas in a city that have vegetation, leaving only the structures.


texts: How to delete parts of a multipolygon inside an sf-object in R?

texts:
I have spatial data (about some glaciers) as an sf-object. My shp-file can be downloaded here:
https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT

Some of the glaciers are separated into different polygons. In the end, I want to resolve the inner boundaries and delete all polygons who are smaller than a certain cutoff (lets say 0,05km²).

So far I managed to merge adjacent polygons using st_union

glaciers <- readOGR("C:/***", "Glacier_Clip_Polygon")
glaciers_sf <- st_as_sf(glaciers)
glaciers_combined <- st_union(glaciers_sf)

Merging them results in a Multipolygon with all inner boundaries resolved. So far so good. It now looks like this:

> glaciers_comb
Geometry set for 1 feature
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 81.38991 ymin: 42.47015 xmax: 81.54986 ymax: 42.57237
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
MULTIPOLYGON (((81.46754 42.47274, 81.46791 42....

But now I don't know how to address the polygons separately or how to delete the ones with an area smaller than my cutoff.
Can anyone help me?


texts: Saving Processing result to GeoJSON string instead of QgsVectorLayer or GeoJSON file

texts:
I run processing.run("native:buffer") and from newly created vector I get all features and export them to string using QgsJsonExporter().
My problem is that newly created vector gives new IDs to all features and they aren't equal to attribute ID's values. And after exporting them to GeoJSON using QgsJsonExporter() I get GeoJSON features with identifiers.
However I can solve this by declaring processing.run("native:buffer") algorithm OUTPUT to .geosjon file and opening it afterwards.
Is there a way that I can get GeoJSON string without IDs without having to save OUTPUT to .geojson and reading it to string afterwards.
vl = QgsVectorLayer(sys.argv[1],"mygeojson","ogr")

def runProcessingNativeBuffer(inputGeoJson, distance):
    processing.algorithmHelp("native:buffer")
    return processing.run("native:buffer",
                         {'INPUT': inputGeoJson,
                          'DISTANCE': distance, 'SEGMENTS': 5,
                          'END_CAP_STYLE': 0, 'JOIN_STYLE': 0,
                          'MITER_LIMIT': 2, 'DISSOLVE': False,
                          'OUTPUT': "TEMPORARY_OUTPUT"})

ats = runProcessingNativeBuffer(vl, sys.argv[2])

vl1=ats['OUTPUT']

features = vl1.getFeatures()
exporter = QgsJsonExporter()
print("GEOJSON",exporter.exportFeatures(features))

IDs in properties are fine I don't need additional IDs in "features"


texts: Using Data Driven Pages - like function for full map layer symbolization?

texts: I have a hexagon layer, and in each hexagon is a count of how many trees are in the hexagon. There are several species within this data. I am trying to create a definition query that will run through each species type the same say Data Driven Pages runs through different polygons. The primary difference is that the map extent will remain the same, but the data being displayed is being changed. In the end, I would like a separate map for each species. Is there a simple way to do this? Do I need to use arcpy or is there a feature within ArcGIS that can handle this?


texts: QGIS Graphical Modeler Add Specific files

texts:
When I create a model in the Graphical Modeler Tool I need to add layers using ‘Add File’.

When I run the script, I have to add the same files each time (Which can lead to adding numerous files). Since none of the layers change, how do I add links to specific files?

The image below shows only two layers which need to be added each time (Listed Buildings and then Building Centre Points). In my detailed script, this happens countless times.


texts: Multiple Field Update on comma separated values

texts:
I am renaming data within columns from names to numbered codes to prepare for importing into a database. This is easy enough if I have single entries within the field, but I'm struggling to update multiple comma separated values.
For single entries, the below example works fine in field calculator:

case when "Enviro" = 'A'
then '118'
when "Enviro" = 'B'
then '119'
when "Enviro" ILIKE '%C%'
then '120'
when "Enviro" ILIKE '%D%'
then '121'
when "Enviro" ILIKE '%E%'
then '122'
when "Enviro" ILIKE '%F%'
then '123'
end

However if the "Enviro" field contains both A and B for example (A,B) then the code for 'B' will automatically overwrite any reference to 'A' in the final output. Ideally, I'd want this particular example to read 118,119 in the final version. Any ideas?


texts: Cascading order of basins (next downstream basin) in QGIS, R, or Whiteboxtools

texts:
I am conducting a hydrologic analysis of DEMs. Suppose that I have the following basins in an area as shown below:.
 I need to get the cascading order of each of the basins. In other words, I need to know the next downstream basin for each one of those basins. Is there a function that can do this in QGIS, Whiteboxtools, or R.
 I am aware of the trace downstream function but it gives the flow line to the outlet, not the next downstream basin.


texts: Converting Mapserver default URL to a more user friendly version

texts:
I use the default MapServer URL to view the map in browser.

http://localhost/cgi-bin/mapserv?map=/home/mwh/osm-demo/basemaps/osm-google.map&mode=browse&template=openlayers&layers=all

How can I convert it to more general style such as http://foo.com ?

I just need to show my map in a specific URL like in maps.google.com or www.openstreetmap.com and nothing more. For example I want to hide my directories showing such as /home/mwh/osm-demo/basemaps/ in the original URL. How can I achieve this?


texts: Syntax of conditional statements in QGIS Field Calculator(s)

texts:
I have a shapefile which I need to modify. There is one column "CLASS" with the attributes: "A", "B" and "C". I need to change "A" into "1", "B" into "2" and "C" into "3". I tried it with case when function. With one case it's not a problem but how does it work with 3 in row?

Is it also possible to build a model with the FieldCalculator function?

I guess I figured it out! Is this right?

CASE
  WHEN "VSt_K" IS 'E' THEN '5'
  WHEN "VSt_K" IS 'A' THEN '1'
  WHEN "VSt_K" IS 'B' THEN '2'
  WHEN "VSt_K" IS 'C' THEN '3'
END

But I can't use this in the modeler FieldCalculator function...


texts:
I have what is hopefully a fairly simple question.

I am trying to run an OD cost matrix for every minute of a day for research I am performing. All of the settings will be the same for each run, except the time will change. I will be automating the process in Python. I know how to create a cost matrix with the settings, add the locations, and solve it all using Arcpy. However, I am hoping to be able to just change the time setting for each run without having to create a whole new cost matrix and add the locations again (similar to what you can do manually in ArcMap).

I have not been able to find a way to do this in Python (if there even is a way).


texts: Georeferencing postcodes into bubble map with Cartodb

texts:
I have recently discovered Cartodb and I am trying to create a bubble/thematic map of postcodes. I have a csv document with one column of postcodes from around the UK. I have uploaded the dataset and when georeferencing the data i choose from the postal code tab and enter United Kingdom as a country.

The issue is that only a handful of the post codes got georeferenced while the rest remain null. The postcodes are valid so not sure why it's not working.

Is there something I am doing wrong?

EDIT
here is an screenshot of the data, you can see a few have their long/lats while others don't


texts: in ArcGIS js API 4.2, ArcGISDynamicMapServiceLayer replaced with MapImageLayer but AGS services are forced to use the previous version

texts:
I'm trying ArcGIS js api 4.2 to generate some services using ArcGIS desktop 10.4.1 and ArcGIS server 10.4.1. but it act so weird for me. i can't see my sercives, all i can see is in this Screen shot that shows just the name of service.

I've done a bunch of search over the internet but there is no solution for this problem apparently. i can't change the ArcGISDynamicMapServiceLayer.js (which is deprecated in js api 4.2) with "/mixin/MapImageLayer.js" (the righ path in js api 4.2). How can i change the wrong path?


texts: Toggle unchecked button is not closing the Socket. QGIS

texts:
I have a QGIS python plugin which works fine when toggle button is checked but when i unchecked it then its not closing the socket why?

here is my code

def show_markers(self):
    s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
    host = socket.gethostbyname('192.168.225.12')
    s.connect((host, port))
    scale=0
    while True:

        if self.iconAction.isChecked():
            print ('Checked')
            m = QgsVertexMarker(self.iface.mapCanvas())
            data = s.recv(SIZE)
            data1 = s.recv(SIZE)
            c = data.decode()
            d = data1.decode()
            x = float(c)
            y = float(d)
            print("printing X :", x)
            print("printing Y :", y)
            rect = QgsRectangle(float(x)-scale, float(y)-scale, float(x)+scale, float(y)+scale)
            me = self.iface.mapCanvas()
            me.setExtent(rect)
            me.refresh()
            m.setCenter(QgsPointXY(x, y))
            m.setColor(QColor(255, 0, 0))
            m.setIconSize(7)
            m.setIconType(QgsVertexMarker.ICON_X)  # or ICON_CROSS, ICON_X
            m.setPenWidth(3)

        else:
            print('Unchecked going to close socket')
            s.close()

A sender computer send the lat/long values & this program marks the lat/long values on map canvas it happens when my toggle button is checked that's good. But when i make it unchecked by calling socket close() method why its not closing the socket (even it prints else part 'Unchecked going to close socket' but s.close() not closing the socket).


texts: Does OSM Road segments contain intersections - reg

texts:
I was searching for OSM documentation online, but I couldn't find one. This is my question,

In OSM(raw data format), does the road segment has intersections included in them. If so, are the intersections are identifiable?. Also, are there any detailed information available about OSM (WiKi is not complete)?.


texts: Creating multiple feature classes from multiple text files using ArcPy

texts:
I have a folder with multiple text files. I am writing a code that reads all of the text files, and creates a point feature class (i.e. a shapefile) for each one.  The text files are all tab delineated with the same set up "Id, name, type, lat, long" I have posted the code I have below.  A shapefile is created for each text file, and the fields are created for the attribute table, however, the points are not created, and the rows are not being inserted into the attribute table.  I had it working for an individual text file but I would like to do it for multiple text files.

#importing arcpy and os
import arcpy, os
from arcpy import env

arcpy.env.overwriteOutput = True

Location = r"C:\temp"

# Set environment settings
env.workspace = Location

for file in os.listdir(Location):
    if file.endswith(".txt"):
        try:
            fcName = os.path.splitext(file)
            outName = fcName[0]  + '.shp'
            outFolder = r'c:/temp/Data'

            #Creating the feature class
            datasetCreate = arcpy.CreateFeatureclass_management(outFolder, outName, "Point", "", "", "")
            cursor = arcpy.InsertCursor(datasetCreate)

            #Creating the fields for the attribute table
            arcpy.AddField_management(datasetCreate, "Name", "TEXT", "", "", 40)
            arcpy.AddField_management(datasetCreate, "Type", "TEXT", "", "", 20)
            arcpy.AddField_management(datasetCreate, "Lat", "DOUBLE", "", "", 20)
            arcpy.AddField_management(datasetCreate, "Long", "DOUBLE", "", "", 20)

            #Open the text file and read it
            textfile = open(file, "r")

            #Create a variable for the headerline
            headerline = textfile.readline()

            #Get rid of the headerline, and split the text file into columns
            valueList = headerline.strip().split("\t")
            print valueList
            Lat = valueList.index("Lat")
            Long = valueList.index("Long")
            Name = valueList.index("Name")
            Type = valueList.index("Type")
            ID = valueList.index("ID")

            #Insert a cursor using arcpy
            for point in textfile.readlines():
                segmentedPoint = point.split("\t")
                latvalue = segmentedPoint[Lat]
                longvalue = segmentedPoint[Long]
                namevalue = segmentedPoint[Name]
                typevalue = segmentedPoint[Type]
                idvalue = segmentedPoint[ID]

                #Create a point object and set the lat and long fields to the x and y coordinates
                vertex = arcpy.CreateObject("Point")
                vertex.X = latvalue
                vertex.Y = longvalue

                #Use the cursor to insert each new row for the attribute table
                feature = cursor.newRow()
                feature.shape = vertex
                feature.setValue("ID", idvalue)
                feature.setValue("Name", namevalue)
                feature.setValue("Type", typevalue)
                feature.setValue("Lat", latvalue)
                feature.setValue("Long", longvalue)
                cursor.insertRow(feature)

        # If an error occurred print the message to the screen
        except:
            # If an error occurred print the message to the screen
            print arcpy.GetMessages()

            #delete the cursor
            del cursor


texts: Landsat 8 Level 2 products. Converting Pan & Thermal bands to spectral reflectance

texts:
We have downloaded around 20 Landsat Level 2 Surface Reflectance scenes for our AoI. The processed bands available in the downloaded scenes are, as you people know - Aerosol,B,G,R,NIR,SWIR,Cirrus.

The unprocessed bands - Pan and 2 Thermal bands, we would like to convert those in to surface reflectance as well as we are working on band ratioing, band combinations & PCA.

I've downloaded LaSRC - which USGS uses to convert DN -> surface reflectance - from their github repository https://github.com/USGS-EROS/espa-surface-reflectance. But couldn't able to run it as it requires few dependencies it seems (our's is windows 10 system. But has QGIS, GRASS and few other OSGEO softwares installed).

My questions are,

1) Can't we convert Thermal bands (which will be in the units of kelvin) to surface reflectance ? If we can't, then how to do band rationing and combinations with thermal bands ? (Tried ENVI's Thermal Correction tool. But not working)

2) Converting Pan band to surface reflectance is straight forward i suppose. Though we are going to use Pan for visual enhancement, do we still need to convert those to surface reflectance ?


texts: OpenLayers - How to create 2 vector layers; 1 draggable and the other is not?

texts: Is there a way to make one vector layer (or a certain kind of feature draggable), while the other types of feature or layer is not?


texts:
I am overlaying some WGS84 (EPSG:4326) coordinates over a few DEM rasters with EPSG:102710 (NAD_1983_STATEPLANE_NEW_HAMPSHIRE_FIPS_2800_FEET)

When they are these two different layer CRS the coordinates overlay perfectly where they are supposed to be, however, if I change the coordinates to have the same EPSG as the DEM rasters, then the coordinates are way off. What is the explanation for this?

The rasters were from this site: http://www.granit.unh.edu/data/search. "Digital Elevation Model - Quad 54" for example.

Unfortunately, I cannot provide the coordinates.


texts: Save categories of graduated symbols, but not their values in QGIS 3

texts: I have about 100 vector layers of different months and years showing EVI anomalies based on county boundaries in Germany. When I use graduated symbology (equal count, 5 classes) for each layer I get a nice map of Germany showing the counties in different colours depending on the EVI anomaly. Now the problem is, that the 5 classes are created depending on the max and min value of each specific vector file and since each month has different max and min values of anomalies, I get different classifications. So the colours are actually misleading when I compare the vector layers with each other. How can I get the same classification for each layer? I have tried to save the style of one layer and applied it to another layer, but then the second layer has the same values as the style layer, i.e. I end up with two identically looking layers - which is not what I want.


texts: Arcmap reads column of integers in .csv as strings instead of numbers

texts:
When reading a column of integers (land cover classes) from a CSV file, ArcMap (10.5) automatically assigns a field type of 'string', rather than a number format (e.g. 'Long'). This makes it so you cannot use that field for symbology or mathematical operations.
Strangely, ArcMap reads the column of numbers from an original CSV as integers, but not from the version modified in Excel.

While I could add a field and then use field calculator to convert the string values to numeric, which is the only solution I could find online, this shouldn't be necessary. There must be a way to get numeric fields in CSVs to be interpreted as numeric, as most CSVs work fine.

Edit: This question differs from the suggested duplicate in that it relates to CSV files and not excel files. Also the CSV file is not being imported with Table to Table, or Excel to Table, but is being used directly as a table in arcpy. The workaround also differs from those other 'solutions', which are also just workarounds. None of the workarounds, including my own, explain the fundamental underlying reason why ArcMap only sometimes decides to read number columns as text. Calling them solutions is like a doctor saying that a painkiller will cure the disease. I don't want to have to convert the file unnecessarily.

Where is the real solution to why ArcMap only sometimes interprets a number column as text?


texts: Difference model between two raster (DEM) how to show values bigger than and smaller than 0,00

texts:
I got to DEM layers over the same area with positive and negative values in them. So I ran a raster calculator raster1-raster2.
Now the raster got values with 0,000659 (This should be counted as 0,00, I just want 2 decimals) and so on.
I would just like to display values that are bigger or smaller than 0,00 either by symbology or by doing a new calcuation and extracting all values that are not equal to 0,00.

Using Qgis


texts: How to get extent of filtered wms layer and How to zoom it with openlayers 4?

texts:
I am looking for zoom to filtered wms layer, but it is not working. I am using Openlayers 4 with geoserver. How to get the extent of filtered wms layer. I tried in API but it is not able applicable. Here is the my code

<!DOCTYPE html>
<html>
  <head>
    <title>Integrated Land Records Management System</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.6.5/css/ol.css" type="text/css">
    <!-- The line below is only needed for old environments like Internet Explorer and Android 4.x -->
    <script src="https://cdn.polyfill.io/v2/polyfill.min.js?features=requestAnimationFrame,Element.prototype.classList,URL"></script>
    <script src="https://openlayers.org/en/v4.6.5/build/ol.js"></script>
    <style>
        .print_button {
            position: absolute;
            z-index: 100;
            right: 20px;
            top: 20px;
        }
        .bhucode {
            margin-left: 30px;
        }
    </style>
  </head>
  <body>
  Select Bhucode
  <select id="bhucode" class="bhucode" >
    <option value="">Select Bhucode</option>
  </select>
    <div id="map" class="map">
    </div>

    <script>
      var wmsSource = new ol.source.TileWMS({
        url: 'http://localhost:9090/geoserver/ilfs/wms',
        params: {'LAYERS': 'ilfs:dharani_samplevillage',
            'TILED': true, 'format': 'image/png',
            'VERSION': '1.1.1', 'SRS':'EPSG:4326',
                "TRANSPARENT": true,
                "ISBASELAYER": true,
                "RATIO": 1,
                "TILED": true,
                "CQL_FILTER": "bhucode=11223344444555 and kid=243"
        },
        serverType: 'geoserver',
        transition: 0
        //crossOrigin: 'anonymous'
      });

      var wmsLayer = new ol.layer.Tile({
          source: wmsSource,
          opacity: 0.7,
          visible: true
      });

      var view = new ol.View({
          projection: 'EPSG:4326',
          center: [0, 0],
          zoom: 18
      });

      var map = new ol.Map({
          layers: [wmsLayer],
          target: 'map',
          view: view,
          renderer: 'canvas'
      });

      console.log("Layer Extent : "+wmsLayer.getSource().getExtent());

      var extent = map.getView().calculateExtent(map.getSize());
      console.log("Extent - "+ extent);

      map.getView().fit(wmsLayer.getExtent(), map.getSize());

      var centerMap = map.getView().getCenter();
      console.log("Center : "+ centerMap);
      map.getView().setCenter(centerMap);

      var srs = JSON.stringify(view.getProjection().getCode());
        //ol.proj.Projection
      console.log("View SRS : "+srs);
      console.log("WMSSource SRS : "+ JSON.stringify(wmsSource.params));
  </script>
  </body>
</html>

Please help me with the above code to get the extent, and how to apply the extent to below.
map.getView().fitExtent(extent, map.getSize());


texts: How to get route between Two lat/long using ArcGIS Map in Android?

texts: What I'm trying to do is have two edit boxes for source and destination. Based on those two locations I need to draw the route on an ArcGIS map in Android.


texts: multipolygon: switch coordinate arrays with javascript/jquery


texts: How to maintain OpenStreetMap labels/details (from QuickMapServices) when exporting from QGIS?

texts:
My question is about OpenStreetMap place labels in QGIS, and how to keep them at a similar size when exported to as they appear in the software. I realise variations of this question have been asked before and it used to be an issue with the OpenLayers plugin but I have yet to find a workable solution.

I am having issues with losing the scale of labels of my backdrop mapping when exporting my QGIS print template to a PDF. I need the labels to maintain their size seen in both the map view (screenshot 1) and print template (screenshot 2). However when I export, the labels shrink to a unusable size (screenshot 3).

Does anyone know a work around to this issue?

I have found you can lower the DPI in the print template export setting which produces larger labels, however this also lowers the resolution of the map to a near unusable level.

Someone had similar problems: Why do labels in exports using HERE tiles appear tiny?

However I have tried this solution (screenshot 4) and unfortunately I don't see any change in the label size to the output? (screenshot 5)

I am not sure if I am putting in this new parameter the URL correctly or not?


texts: Do I have to transfer gwc-layers from my GEOSERVER_DATA_DIR directory to a production server

texts:
I am working with Geoserver and I know that it is possible to transfer the contents to the GEOSERVER_DATA_DIR between servers. I have already moved my GEOWEBCACHE_CACHE_DIR to a different location so that the tiles are not transferred. I am using the embedded geowebcache in Geosserver version 2.2.5.

I however have a number of xml files in GEOSERVER_DATA_DIR/gwc-layers:

LayerGroupInfoImpl--570ae188_124761b8d78_-7fab.xml
LayerGroupInfoImpl--570ae188_124761b8d78_-7fac.xml
LayerGroupInfoImpl--570ae188_124761b8d78_-7fad.xml
LayerInfoImpl-18e1ce0d_13ef93f68ce_-7ffd.xml
LayerInfoImpl-1c3a926c_13ef90929b2_-7a75.xml
LayerInfoImpl-2fbcd275_13e632965e7_-7feb.xml
LayerInfoImpl--4a9eec43_132d48aac79_-7fff.xml
LayerInfoImpl--570ae188_124761b8d78_-7fae.xml
LayerInfoImpl--570ae188_124761b8d78_-7fb0.xml
LayerInfoImpl--570ae188_124761b8d78_-7fb4.xml
LayerInfoImpl--570ae188_124761b8d78_-7fb6.xml
LayerInfoImpl--570ae188_124761b8d78_-7fb8.xml

So my question is - are these files generated automatically if I tranfer the contents of GEOSERVER_DATA_DIR or do I have to include them in the transfer?


texts: Applying filter on the geometry for a given feature?

texts:
I have multiple layers, with multiple features, where I am interested in extracting information from a given feature within a rectangle. I am able to achieve this using the following filter,

currLayer.SetSpatialFilterRect(minx, miny, maxx, maxy)

This works very nicely when the feature has a single geometry associated with it. However, some of the features contain multiple geometries, which is evident when calling

feat.DumpReadable()

where some features contain POLYGON entries with 2+ geometries. And since the filter intersects one (or more) of the geometries, all are returned.

I was therfore wondering, is it possible to apply the filter to a given features geometries? Instead of to the POLYGON entry of the feature itself.

Additional information: Below is a link to an image showing the issues. The green bounding box is the applied spatial filter, where two objects are correctly selected and indicated by X. And since the box intersects two of the blue areas, they are also included. However, all the defined geometries within the blue areas are neither intersected or included within the green box. But due to being part of the overall feature, they are not filtered.

ar = []

currLayer = ds.GetLayerByName("layername")
currLayer.SetSpatialFilterRect(minx, miny, maxx, maxy)
currLayer.GetFeatureCount()
for feat in currLayer:
    multi_geom = feat.GetGeometryRef()
    for iPnt in range(multi_geom.GetGeometryCount()):
        pnt = multi_geom.GetGeometryRef(iPnt)
        ar.append(pnt.GetPoints())


texts: Accessing raster attribute table in QGIS?

texts:
I have some Python code that makes Raster Attribute Tables for GeoTiffs using GDAL.

In ArcGIS I can see these attribute tables fine and they behave how I like (as classifications) but in QGIS the whole concept of an attribute table seems to be missing.

Is there any way to access a RAT in QGIS?


texts: Creating threshold for NDVI using Google Earth Engine?

texts:
I want to create a threshold for the NDVIs in this collection.  I would like to say something along the lines of:

If pixel is less than 0, assign pixel a value of 0; otherwise assign the pixel a value of 1.

This should turn each raster NDVI into a raster that only contains ones and zeroes.

//NDVI functions for all satellites
var mergedAllFunction = function(mergedCollection) {
  var ndvi = mergedCollection.normalizedDifference(['NIR', 'red']).rename('NDVI');
  return mergedCollection.addBands(ndvi);
};

// find all data and filter them by date and add NDVI
var mergedCollectionNDVI = mergedCollection
    .filterDate('1988-01-01', '2018-10-01')
    .filterBounds(studyarea)  // filter to area-of-interest;
    .map(mergedAllFunction).select('NDVI');
print ('mergedCollectionNDVI', mergedCollectionNDVI);


texts: Create field with length of lines intersecting with polygons

texts:
In QGIS I`m trying to find an expression (or anything really) that could could do the following: given a layer of polygon and a layer of line that intersects it, I wish to be able to populate or create a field in polygon with the length of the the line where it intersects with the polygon, the reference for the length being always the start point. But this value can only be integer and only store even values.
This image illustrate what I'm trying to achieve.


texts:
I'm working with the shapefiles for the US census 1999 msa files which can be downloaded here:

https://www2.census.gov/geo/tiger/PREVGENZ/ma/ma99/cm_sa_99_shp.zip

Unfortunately, there is no CRS attached with these shapefiles. After reading other questions, it seems that guess and check for the CRS is the best strategy. EPSG 4326 seems to work well, but after loading the shapefile with the R sf package, there seems to be some slight distortion in the shapefile after using sf::st_set_crs with 4326. My question is should I worry about this distortion and is there a better CRS (or method to find the CRS)?


texts: Polygon shapefile field calculator problem

texts: My project is set in WGS84 only, and I have checked 'enable on the fly'. I have created a polygon shapefile and my problem is that when i use field calculator in the attribute table on my polygon shapefile (settings are absolutely identical as above) I use the geometry: $area to create a column for area_m2 but all I get is 0, Where the other column I have made has values for perimeter. Any ideas?


texts: OpenLayers: Draw feature of defined size (in meters)

texts:
I'm trying to draw a rectangle polygon on the map having a defined size in meters.

Background is: I want to show the user a rectangle of a defined paper size (e.g. DIN A4) on the map to choose a region to print for a specific scale (e.g. 1:25000).

My math to calculate the size in meters that would be visible on the paper is:

widthInMeters = paperWidthInCm * (scale/100)

Given the width of A4 landscape (29.7cm) I'd get 7425 meters (29.7 * (25000/100)).
Height would be 7425 * 0.7070 = 5249m

Now I want to draw a rectangle of this size (7425m x 5249m) on the map.

My first try was to just set the extent accordingly because the coordinates units (Web Mercator EPSG:3857) are in meters. But that, I think, would only work directly on the equator. Further up north or south this does not work at all because of projection distortion.

Given a defined width and height in meters and a center (or top left) coordinate how can I create an appropriate geometry covering that area?


texts: Can't load PyQGIS on QGIS3 due to problem with API version of sip module

texts:
Again having issues with building QGIS 3 on Manjaro Linux, from AUR. After solving my previous problem (QGIS 3.0 can't load plugin processing, no module named 'osgeo', on Manjaro Linux), I now get the following error message:
Couldn't load PyQGIS.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
  File "/usr/share/qgis/python/qgis/__init__.py", line 72, in
    from qgis.core import QgsFeature, QgsGeometry
  File "/usr/share/qgis/python/qgis/core/__init__.py", line 34, in
    from qgis._core import *
RuntimeError: the sip module implements API v11.0 to v11.3 but the within renderer)

.. versionadded:: 2.6 module requires API v12.4

Python version:
3.6.4 (default, Jan  5 2018, 02:35:40)
[GCC 7.2.1 20171224]

QGIS version:
3.0.1-Girona 'Girona', exported

Python path:
['/usr/share/qgis/python', '/home/thiago/.local/share/QGIS/QGIS3/profiles/default/python', '/home/thiago/.local/share/QGIS/QGIS3/profiles/default/python/plugins', '/usr/share/qgis/python/plugins', '/home/thiago/anaconda3/lib/python3.6/site-packages', '/usr/lib/python36.zip', '/usr/lib/python3.6', '/usr/lib/python3.6/lib-dynload', '/usr/lib/python3.6/site-packages']

I have recompiled gdal, python-gdal, have used pip to install the latest version of sip, but I still get the same error, after several rebuilds. This error means QGIS starts without python support, making it pretty much useless for any analysis.


texts:
I have a Geoeye 1 satellite photo that's a 4-band multispectral B,G,R,NIR, 11 bit raster.

Ultimately I want to publish this raster as a web service (using Mapserver).
I opened up the raster in Qgis just to look at it and I noticed that it took an unusually long time to load it, which does not bode well for an quick and responsive web service.

So I concluded that it would be good if the raster could be transformed into a "normal" RGB raster.

Is it possible to process a RGB+NIR raster into a "standard" RGB raster? How can this be done?

Here's what gdalinfo reports on the bands of the raster http://pastebin.com/1bQ9Lh16
As you can see there's a fourth band with no color intepretation, yet if I would just remove this band I end up with a quite blue image. So somehow it seems to contribute to the visual appearance nevertheless.


texts: Generating SLD files from Autodesk's layer

texts:
I have Autodesk layer files (with .layer suffix).

I'd like to create SLD files based on these layer's information (e.g. the colors). I tried to generate the SLD XML file by using xslt and parsing the layer files.

Does anyone have a better solution for generating SLD files?


texts: Making feature layer based on attributes (2 columns) in ArcPy?

texts:
I'm new to Python, but know this tool can solve my problem. I have a large dataset and I want to extract layers based on the attributes of two columns.

For example, my dataset has columns A-N. In this exercise, columns D and C are of relevance. These columns both have variables ranging from 1-9. I want to loop through this data set and extract every combination based on these two columns (D=1 and C=1, D=2 and C=1, D=1 and C=2, etc.).

I assume I need to use a combination of the 'select by attribute' tool and 'make a feature layer' tool, but am having trouble putting together the two.

Can you provide a basic script structure for this solution?


texts: QGIS error when accessing data from Database via eVis and ODBC

texts:
Using the "eVis Database Connection" I can access an ODBC data source which contains records with Lat and Long fields. The connection is made OK, and the "Database File Selection" window opens where I can define the name of the new layer and the x and y coordinate fields. Having done this, clicking OK gives the error message:

Layer is not valid : The layer file:///C:/Users/ABCDE/AppData/Local/Temp/QGIS3.HJlEUi?
delimiter=%09&delimiterType=regexp&xField=WP_Lat&yField=WP_Long is not a valid  layer
and can not be added to the map. Reason:

Looking at the QGIS Log Messages, under the "Delimited Text" tab I can see:

2020-01-15T18:39:48  WARNING  Errors in file C:/Users/ABCDE/AppData/Local/Temp/QGIS3.HJlEUi
2020-01-15T18:39:48  WARNING  X field is not defined in delimited text file
2020-01-15T18:39:48  WARNING  Y field is not defined in delimited text file

I have looked at the temporary file QGIS3.HJlEUi it appears to be a delimited text file in exactly the correct format:

WP_Lat  WP_Long

-15.765315  40.515876

-15.779443  40.534107

-15.816472  40.628185

-15.816369  40.611462

-15.816465  40.628178

-15.816469  40.628175

Can any one see what might be going wrong?

QGIS 3.8.3
eVis 1.1.0 (core plugin)
Windows 10


texts: How to undo ArcGIS "Move Edits to Base" mistake?

texts:
I'm trying to help a client who mistakenly checked the "Move Edits to Base" option when versioning a feature class in ArcGIS.  They would now like to fix this mistake.  However, they don't want to lose any of the versioned edits in that feature class, and they cannot reconcile and post everything in it.  There are a couple of years of work in there from several thousand long-lived versions.

Does anyone know of a way to accomplish this?

Thanks!


texts:
I am trying to use file geodatabase (.gdb) data in QGIS.

I have tried opening the dataset in ArcMap but I keep getting an error message so trying in QGIS.  I have got as far as opening and viewing the data however, the data in the attribute tables are just numbers.  These numbers relate to feature type names that would show up in ArcMap via the 'show domain and subtype descriptions' button.

Is there a way of doing this in QGIS?  The data is pretty much useless without this link.


texts: Google Earth Engine: get image value on click

texts:
When clicking on a GEE map, I would like to get the value of an Image's (1-band raster) pixel.

I'm using Map.onClick to extract values from point on a map. According to GEE examples (e.g. https://developers.google.com/earth-engine/tutorial_api_03), xx.reduceRegion(ee.Reducer.mean(), clickpoint, scale) should work for ImageCollection objects but does not seem to work for Image objects.

// Load 1-band image
var dataset = ee.Image('WWF/HydroSHEDS/03VFDEM').select('b1');
print(dataset);
Map.style().set('cursor', 'crosshair');
Map.addLayer(dataset);

// Following is adapted from https://code.earthengine.google.com/?scriptPath=Examples:User+Interface/Ocean+Timeseries+Investigator
// Set a callback function for when the user clicks the map.
var header = ui.Label('Some text', {fontWeight: 'bold', fontSize: '18px'});
var toolPanel = ui.Panel([header], 'flow', {width: '400px'});

Map.onClick(function(coords) {
  // Create or update the location label (the second widget in the panel)
  var location = 'lon: ' + coords.lon.toFixed(4) + ' ' +
                 'lat: ' + coords.lat.toFixed(4);
  var click_point = ee.Geometry.Point(coords.lon, coords.lat);
  var demValue = dataset.reduceRegion(ee.Reducer.first(), click_point, 90);
  var demText = 'Habitat suitability: ' + demValue;
// var demText = 'Habitat suitability: ' + demValue.get('b1');
  toolPanel.widgets().set(1, ui.Label(location));
  toolPanel.widgets().set(2, ui.Label(demText));

  // Add a red dot to the map where the user clicked.
  Map.layers().set(1, ui.Map.Layer(click_point, {color: 'FF0000'}));

});

// Add the panel to the ui.root.
ui.root.add(toolPanel);

Using demValue directly or demValue.get('b1') doesn't change anything.

The resulting panel gives (-86.6746,35.3036 being a random point):

  Some text
  lon: -86.6746 lat: 35.3036  Habitat suitability:
  ee.Dictionary({ "type": "Invocation", "arguments": { "image": {
  "type": "Invocation", "arguments": { "input": { "type": "Invocation",
  "arguments": { "id": "WWF/HydroSHEDS/03VFDEM" }, "functionName":
  "Image.load" }, "bandSelectors": [ "b1" ] }, "functionName":
  "Image.select" }, "reducer": { "type": "Invocation", "arguments": {},
  "functionName": "Reducer.first" }, "geometry": { "type": "Point",
  "coordinates": [ -86.67460937499999, 35.30364788916715 ] }, "scale":
  90 }, "functionName": "Image.reduceRegion" })

GEE script: https://code.earthengine.google.com/1204e25cbbdb782d9af4c7944d047bf2


texts: How can I use map service ( published on Arcgis Server 10.4.1) from another computer?

texts: With ArcGIS Server 10.2, the map service url is in a format like http://{ip}+ port/arcgis/rest/services/yourservice, in Arcgis Server 10.4.1, the map service url replace the ip with local domain, for now I can only load the service on my local machine using Arcgis Javascript. How can I make the map service public to all the local network users, like a different computer on the same network?


texts: Cumulative cost upstream a raster stream network with GRASS

texts:
I have a river network in raster and locations of many dams on the river.

In the image above, the BLUE point represent one random cell. Dam 1 is downstream to the cell.

I am looking to calculate the cumulative installed capacity of all dam upstream each cell. Of course, the installed capacity data are collected.

In the example, the sum of Dam2 and Dam3 capacity should be assigned to the blue cell. And I want the result in raster.

Does anyone know how to solve the problem using GRASS?


texts: Align descriptions after labels in ArcGIS legend

texts: I have produced a legend which includes both a label and description for each symbology. I would like the descriptions to align despite the label lengths differing between layers. In this example, "example 1" is part of a dataset which has much longer layer names, hence a wider gap. I have tried wrapping text but can't see a means of separating the label and description. I have also tried offsetting the description but this takes quite a lot of manual time to repeat and it seems like there must be a quicker and more efficient way of doing it.


texts: Can you change the location of .qgis from the home directory to another folder on Windows in QGIS1.8?

texts: Is it possible to change the location of the .qgis folder to another location ( for example on a network drive H:\ ) in Windows?


texts: How to change the style of Openlayers 3 zoom box?

texts:
I want to zoom to an area selected by drawing a box on the map following this example:

http://openlayers.org/en/latest/examples/box-selection.html

Working fine. Now I want to change the line color and background of the box. Some examples have this code (edited to show all my working CSS changes but ol-dragbox):

<style>

.ol-overviewmap {
    bottom: 30px !important;
}

.ol-control button {
    font-size : 18px;
}

.ol-zoomslider {
    top : 5.5em;
}

.ol-dragbox {
    background-color: rgba(255,255,255,0.4);
    border-color: rgba(100,150,0,1);
}

.ol-scale-line {
    background: rgba(0, 60, 136, 0.5);
    padding: 2px;
}

</style>

... it is changing all classes but not ol-dragbox.

I'm curious because in this example I can't see anything adding this kind of behaviour ( zoom to box ) and nor in my map even it is working well:

My map:

map = new ol.Map({
    layers: [ landLayer, routesLayer ],
    target: 'map',
    renderer: 'canvas',
    loadTilesWhileAnimating: true,
    loadTilesWhileInteracting: true,
    controls: ol.control.defaults().extend([
       new ol.control.ScaleLine(),
       new ol.control.ZoomSlider(),
       new ol.control.MousePosition({
           undefinedHTML: '',
           projection: 'EPSG:4326',
           coordinateFormat: function(coordinate) {
               return ol.coordinate.format(coordinate, '{y}, {x}', 4);
           }
       })

    ]),
    view: theView
});

I not added any control or interaction to create this zoom-to-box feature.
Inspecting the OpenLayers CSS file ( v3.15 ) I can't find any style to .ol-dragbox. Where this blue box came from and how can I change it?


texts: Extracting raster statistics *by pixel category* on QGIS or GRASS GIS?

texts:
i know i can churn out raster statistics with r.report in GRASS GIS and QGIS, but what i want is the statistics of manually-defined categories of pixel ranges.

what i have is a clipped raster (of a single band) with a range of values that has not been classified. i only know how to categorize the pixels by going to the style tab, and what i want is a count of pixels in each category that i have defined.

to elaborate further, i about ~80,000 pixels and a maximum pixel value of ~25. i want to categorize my pixels like this:

(1) no value - pixel value range 0.0 - 1.5
(2) low value - pixel value range 1.6 - 6.5
(3) medium value - pixel value range 6.6 - 15.0
(4) high value - pixel value range 15.1 - 25.0

now i want a count of the number of pixels in each category. how do i produce that statistical report? i tried the "zonal raster statistics" tool but there's something amiss here. (See attached image for an idea of what i get with a simple r.report of my raster)

I am an amateur in GIS!


texts: Forbidden error with WPS process chaining

texts:
We have published some geo-processing map services (it’s called WPS in OGC terms) in ArcGIS Server. So I have looked into GeoServer to see if we can do the same. I’ve found following issues. Not sure if it’s related to “HTTP” protocol.

When tested “Process Chaining”, got 403 error. To reproduce the issue, follow the steps below:
-Go to the Demo WPS Request Builder interface, and select gs:Clip from the first combo-box Select the Vector Layer sf:streams as input.
-For the clip input select SUBPROCESS and then click on Define/edit. Another WPS Builder dialog window should appear, but got 403 error

When tested WPS “Execute process”, got “308 Permanent Redirect” response. To reproduce this issue, follow the steps below:

-Go to the Demo WPS Request Builder interface, and select JTS:intersection from the first combo-box.
-The interface asks for two mandatory parameters, a and b, which must be geometries. Choose application/WKT for all input and output parameters and insert the following inputs for a and be:
        a --> POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))

        b --> POLYGON((5 5, 5 15, 15 15, 15 5, 5 5))

-Click on Execute Process after having selected application/WKT as type of the generated result, we expect to get back a polygon representing the intersection between a and b which should be:
POLYGON ((5 10, 10 10, 10 5, 5 5, 5 10))
Instead, 308 Permannent Redirect response is returned.


texts: MakeFeatureLayer throwing error 00840: The value is not a feature layer

texts:
I'm working through tutorials from the book "GIS Tutorial for Python Scripting," by David W. Allen, and I have this code pretty much exactly the way the book wants it (I altered one variable name from LayerSStations to LayerStation because the double-S was giving me problems seeing it) and every time I run it it throws the error "00840: The value is not a feature layer."

Here's my code, up to the point it throws the error:

import arcpy

#create a map document object for the map document
thisMap = arcpy.mapping.MapDocument(r"C:\Users\Andy\Documents\GIS DataBase\EsriPress\GISTPython\Maps\Tutorial 3-3.mxd")

#Create data frame object
myDF = arcpy.mapping.ListDataFrames(thisMap)[0]

#Create list object of the map elements
myMapElements = arcpy.mapping.ListLayoutElements(thisMap)

#Get a list of all layers in the table of contents
#Use a for loop to find layers Parcels and SamplingStations
#assign them to new variables for easier reference
myLayers = arcpy.mapping.ListLayers(myDF)

for lyr in myLayers:
    if lyr.name == "Parcels":
        print "Parcels have been found."
        layerParcels = lyr
    elif lyr.name == "SamplingStations":
        print "Sampling Stations have been found."
        layerStation = lyr

#create feature layer to hold only Operational sampling stations
arcpy.MakeFeatureLayer_management(layerStation,"ActiveStations_lyr",'"Status=\'Operational\'')

It's failing on that last command, the MakeFeatureLayer_management line. Any idea what I'm doing wrong?

I've been over and over the variable names looking for wrong things, I have done Rubber Duck debugging, and I can't nail down the problem.

EDIT: Added more of the original code to the snippet.


texts: Adding curved features to map with ArcGIS API for JavaScript?

texts:
I've tried without much success to add a curved line feature to the ArcGIS JavaScript map.

I created a polygon with a curved edge in ArcMap:
curved line http://imageshack.com/a/img538/9702/p6mj9B.jpg

However, when I load this same feature in the web API, it becomes segmented:
segmented curve http://imageshack.com/a/img673/9196/rPnJ4l.png

The problem started when I tried to create a curved feature from the JavaScript API itself, which only seems to allow full circle features (not circular arcs).  This makes me think that curves are still unsupported in the JavaScript API.
Do you know any way around this? Do you have any other ideas on alternatives I can use to create and correctly view curved features?


texts: Touristic poi-data (api or dump) with popularity measure?

texts:
I'm looking for freely available touristic (hotel, restaurants, bars, attractions, etc.) poi-data that have some measure of popularity.

With popularity I mean, explicit popularity like on Tripadvisor (i.e: #10 out of 1500 attractions in Paris), but no api/ dump is available for TA. Or otherwise more soft-measures like the nr of checkins (Foursquare for example)

Anyone knows of any data-set (dump or api) which makes this type of info available? First use would be to display markers on a map, where only the most popular are shown. (as showing all 5000 available bars in say Paris is not very practical)

See list below of some sources I'm checking out

Factual Global PLaces (dataset from former SimpleGeo): likely no popularity, not tried.
openstreetmap (not tried yet, perhaps my next check)
foursquare (has an api with implicit popularity (nr of checkins)) but hard/impossible to query the api to get complete coverage (besides it's against their ToS
Facebook Graph api looking for places (not tried)


texts:
I am looking for help digitizing or adding info from survey plans into QGIS.

In this example I'm trying to add the lease boundaries highlighted in yellow and the pipeline right of ways to the right of it.

Right now I just eyeballed a feature (well centre) so I used the advanced digitizing tools to draw from it based on the lengths and angles noted in the survey drawing, but I was wondering what the proper way to do this would be.

Part of the problem is I do not know how to tell what angle the lease is sitting on.

Maybe QGIS is not the best open source tool for this task.

Survey Drawing


texts: I have a table, this table contains attributes written in GML and I want to transfer these GML to geometry in QGIS

texts: Some attribute written as a .GML file in PostgreSQL and I want to transfer these GML to polygon, so that the table becomes a vector table.


texts: Creating Protocolbuffer Binary Format (PBF) file using QGIS?

texts:
I have a set of geospatial data, lines (streets) in geojson and shp file format which which i want to convert to PBF (Protocolbuffer Binary Format).

I am presently working with QGIS.

Please does anyone have an idea of how i can go about this and what tool(s) can be used?


texts: Updating existing GeoServer ImageMosaic with GeoTIFF using REST API

texts:
While trying to update an existing ImageMosaic in GeoServer using the REST API, following the documentation (provided here http://docs.geoserver.org/stable/en/user/rest/examples/curl.html#uploading-and-modifying-a-image-mosaic)

using the following command

curl -u admin:geoserver -XPOST -H "Content-type: text/plain" -d "file:///path/to/file/img.tiff" "http://localhost:8080/geoserver/rest/workspaces/opengeo/coveragestores/klicfinal/external.imagemosaic"

I get an empty response, while my ImageMosaic remains un-updated. Does anyone have a clue as to what could be causing this?

EDIT: I've already tried using the recalculate parameter. Furthermore, after closer examination, the shapefile of the imagemosaic has been updated and already shows the location and the reference to the new geotiff file. However, the newly added image does not show up in the WMS layer preview, even after disabling caching.

Even the GML that you can get using http://localhost:8080/geoserver/rest/workspaces/opengeo/coveragestores/klicfinal/coverages/klicfinal/index/granules.xml?
shows the updated reference to the newly uploaded GeoTIFF.

Hostname was NOT found in DNS cache
Trying 127.0.0.1...
Connected to localhost (127.0.0.1) port 8080 (#0)
Server auth using Basic with user 'admin'
POST /geoserver/rest/workspaces/opengeo/coveragestores/klicfinal/external.imagemosaic HTTP/1.1
Authorization: Basic YWRtaW46Z2Vvc2VydmVy
User-Agent: curl/7.35.0
Host: localhost:8080
Accept: */*
Content-type: text/plain
Content-Length: 83

upload completely sent off: 83 out of 83 bytes
HTTP/1.1 202 Accepted
Date: Tue, 11 Apr 2017 14:33:17 GMT
Server Noelios-Restlet-Engine/1.0..8 is not blacklisted
Server: Noelios-Restlet-Engine/1.0..8
Transfer-Encoding: chunked

Connection #0 to host localhost left intact


texts: How to obtain geometric centroid for DataFrame where each row is a unique vertex?

texts:
I have looked at several related questions, but have not found a working solution beginning with a DataFrame of vertices for multiple polygons. There are existing solutions that address shapefiles where each polygon is contained in a single row.

My issue is the output of gCentroid is clearly weighted toward uneven boundaries. It seems to be generating the mean of the polygon vertices rather than the expected centroid.

library(ggplot2); library(maps); library(sp); library(rgeos)

county_df <- map_data('county')  #mappings of counties by state
county_df <- subset(county_df, region %in% c("illinois", "indiana", 'michigan', 'minnesota', 'wisconsin'))   #subset just for GL
county_df$county <- county_df$subregion
county_df_plot <- county_df

coordinates(county_df) <- c("long", "lat")

# Get centroids
ctrs <- lapply(unique(county_df$group),
    function(x) gCentroid(SpatialPoints(county_df[county_df$group==x,])))
ctrsout <- setNames( ctrs , unique(county_df$group ) )

# Create data frame
centroid_df <- do.call(rbind, lapply(ctrsout, data.frame, stringsAsFactors=FALSE))
uniq_logical <- !duplicated(county_df$group)
centroid_df$state <- county_df$region[uniq_logical]
centroid_df$county <- county_df$county[uniq_logical]
names(centroid_df) <- c("longitude", "latitude",'state', 'county')

ggplot(county_df_plot, aes(long, lat)) +
  geom_polygon(aes(group=group), colour='black', fill=NA) +
  geom_point(data=centroid_df, aes(longitude, latitude), size=1)


texts:
I have a satellite image of only one part of some city with a lot of vegetation. I need to find out where is that image taken.

Is there any way to find out where it is taken on Google Maps, Google Earth or something similar?

It is a .jpg image so I don't know how to start search.


texts:
I'm in need to process a large batch of Landsat images with MapTiler and since the version I'm using (1.0 Beta 2) does not support multiple input files, I created a simple Python script that calls gdal2tiles for every file.

The problem is that using, apparently, the same settings, I get different coordinates.

This is the correct output from MapTiler:

 <?xml version="1.0" encoding="utf-8"?>
    <TileMap version="1.0.0" tilemapservice="http://tms.osgeo.org/1.0.0">
      <Title>mar_apr.tif</Title>
      <Abstract></Abstract>
      <SRS>EPSG:900913</SRS>
      <BoundingBox minx="-26.36147186954573" miny="-65.83667775218426" maxx="-21.96598523780907" maxy="-62.31334995428857"/>
      <Origin x="-26.36147186954573" y="-65.83667775218426"/>
      <TileFormat width="256" height="256" mime-type="image/png" extension="png"/>
      <TileSets profile="mercator">
        <TileSet href="6" units-per-pixel="2445.98490468750010" order="6"/>
        <TileSet href="7" units-per-pixel="1222.99245234375010" order="7"/>
        <TileSet href="8" units-per-pixel="611.49622617187504" order="8"/>
        <TileSet href="9" units-per-pixel="305.74811308593752" order="9"/>
        <TileSet href="10" units-per-pixel="152.87405654296876" order="10"/>
        <TileSet href="11" units-per-pixel="76.43702827148438" order="11"/>
        <TileSet href="12" units-per-pixel="38.21851413574219" order="12"/>
        <TileSet href="13" units-per-pixel="19.10925706787110" order="13"/>
      </TileSets>
    </TileMap>

And this is the incorrect output from gdal2Tiles:

<?xml version="1.0" encoding="utf-8"?>
    <TileMap version="1.0.0" tilemapservice="http://tms.osgeo.org/1.0.0">
      <Title>mar_apr_CLIP1.tif</Title>
      <Abstract></Abstract>
      <SRS>EPSG:900913</SRS>
      <BoundingBox minx="-25.32504747468283" miny="1.94797423573108" maxx="-21.33272919903126" maxy="5.10698976386579"/>
      <Origin x="-25.32504747468283" y="1.94797423573108"/>
      <TileFormat width="256" height="256" mime-type="image/png" extension="png"/>
      <TileSets profile="mercator">
        <TileSet href="6" units-per-pixel="2445.98490468750010" order="6"/>
        <TileSet href="7" units-per-pixel="1222.99245234375010" order="7"/>
        <TileSet href="8" units-per-pixel="611.49622617187504" order="8"/>
        <TileSet href="9" units-per-pixel="305.74811308593752" order="9"/>
        <TileSet href="10" units-per-pixel="152.87405654296876" order="10"/>
        <TileSet href="11" units-per-pixel="76.43702827148438" order="11"/>
        <TileSet href="12" units-per-pixel="38.21851413574219" order="12"/>
        <TileSet href="13" units-per-pixel="19.10925706787110" order="13"/>
      </TileSets>
    </TileMap>

Does MapTiler perform any other operations besides running gdal2Tiles?


texts: How to speed up shortest path routing using pgRouting?

texts:
I need to develop the shortest path routing query in PostgreSQL. I have my own map data that is a shapefile. From that shapefile I have made PostgreSQL database. It is a large database with around 60k entries. I have made one edge table from the dataset and then created the topology. Lastly tried with astar and dijkstra algorithm (pgr_astar and pgr_dijkstra). It is giving the shortest path but the pgr_astar or pgr_dijkstra query takes at least 3 to 4 sec to provide the result. But I need it faster, may be in msec.

There are lots of applications with routing in market and in those cases routing results come within fraction of seconds. What are the tricks to make it faster?

I have checked a few tricks to speed up. There are highway hierarchies, busy node hierarchies. But I can't use these to my application as most of my roads are forest roads. And my start and/or destination point is a point somewhere in between forests; not a city centre of a big city (so, not a busy node).

I tried to generate clusters and then to pre-compute routes between clusters. But thus the database is growing significantly. So, I am not sure whether it will increase the speed or not.

Can you give any more easy idea to solve this problem?


texts: Distance between the rotated label and the line in map composer

texts:
When trying to print rotated map in composer manager labels are shifted:

In map canvas everything is OK!

I need to save again my layer in project projection for avoid this offset.

It was tested in version 2.6 there were no offset:

here is my loaded files: the line and point vector layers with their styles, project file.

The "trick" with adding the map object in print composer, rotate the map canvas, then rotate the map object also in the composer is very unconvinient for me, because I use atlas generation and creating new composer is unnecessary.
When I choose the map that I want to print in composer atlas list menu, I have already defined rotation angle value (calculated from projection formula),
and even if I manually set it to 0 and then my angle again, there is still offset.

Is it possible to get rid of this distance without creation new composer?

Some ideas?)


texts:
I have a csv file Author, Name, Geometry.

The Geometry column in the csv has rows like this

{"type": "Point", "coordinates": [11.6285463, 42.4193742]}

How do I import the json as geometry?

I can import the csv to an existing table using the following:

copy locations FROM path/to/file/pleiades-locations.csv' DELIMITERS',' CSV HEADER;

I have the geometry column set as json, but no geospatial data is available when I connect with QGIS, do I have to do some conversion?


texts: Identifying port numbers for ArcGIS Online Basemap?

texts:
we have web map application on DMZ server, it was created by "ArcGIS API for JavaScript". our client needs to block all port numbers that access on Internet except URL of (arcgisonline rest services) for reasons related to Security Policy.

So is there any idea how to find port numbers or IP of arcgisonline_basemap