texts: Conditional dissolve of adjacent polygons

texts:
I have a shapefile with many non-overlapping polygons, touching each other (see sample screenshot below):

All these polygons have a field called 'avg_ht' in the attribute table, which is a 32-bit float (average height of trees).
Some of these polygons have to be merged with others based on the following criteria. Consider polygon1 and polygon2. They should be merged if:

They are adjacent to each other
abs(polygon1.avg_ht - polygon2.avg_ht) <= 1.0

Condition #2 above basically checks if the difference between the average heights of the candidate polygons is less than 1.0 units.

What is the best way to do this? I have license for the eliminate tool but am not sure if it is the best way.


texts: Modifying values in raster on basis of value of overlapping cells in another raster?

texts:
In ArcGIS 10.1 (advanced license) I have a slope raster, and a raster representing a road network. In the road raster, each pixel is given value 0.5, and the areas outside roads are set to NoData. Now, for some purposes I have to accomplish, I need to half the slope value in those cells that fall within roads. I cannot multiply SLOPE by ROAD rasters because I would end up with NoData in those areas outside the roads.

How can I keep the slope value unchanged for the cells outside the roads?


texts: Python - Calculate image extent based on center coordinate, number of pixels and rotation angles

texts:
I have images with a center geo-coordinate and rotation angles in x- and y-direction. Now I want to georeference them by setting up a GeoTransform object. For non-rotated images, this is easy, using
ul_x = cent_x + ((cols / 2.) * image_resolution)
ul_y = cent_y + ((rows / 2.) * image_resolution)
geotrans = (ul_x, image_resolution, 0, ul_y, 0, -image_resolution

But how can I calculate the actual upper left corner coordinate in case of a rotated image?
I guess there is a way to somehow utilize affine.Affine in a similar way as given in the accapted answer of Retrieve pixel value with geographic coordinate as input with gdal, but I don't know how.


texts: ArcGIS Pro Topology - Fix multiple errors in one step through Inspector

texts: Does someone know if it is possible to select multiple "error-features" in ArcGIS Pro Topology error inspector and fix all of them with the same solution (e.g. create new feature)?  Currently in ArcGIS Pro 2.2.3 it does not seem possible.


texts: How can I make LAS data output from an RPAS survey processed in Agisoft Metashape, to work in lidR with catalog_apply?

texts:
I am creating a Normalized Canopy Height Model using a ground DTM from lidar and a set of surface elevation LAS files from an RPAS (drone) processed and output with Agisoft Metashape (these will be the vegetation). I want to filter/decimate the RPAS LAS files (stored in a LAS Catalogue) and then subtract the lidar DTM in a function called with catalog_apply in lidR.
It works for some chunks (all with warnings: WARNING: data type 22 of attribute 1 ('normals') is deprecated), but most are empty. Finally it bails with an error cannot allocate vector of size xxxMB.
I do not know if its a code error or a data error. I'll provide the code in case someone can see an error. Otherwise is there a way to "fix" the las data? (I already tried las2las to change from 1.2 to 1.4 for the LAS format).
library(lidR)

RPAS_LASCatalogue <- readLAScatalog("D:/Projects2020/p2020_1021_RPAS_EFI_test/RGB/RPAS_LAS", select = "xyzc")
LIDAR_LASCatalogue <- readLAScatalog("D:/Projects2020/p2020_1021_RPAS_EFI_test/lidar_LAS/original", select = "xyzc", filter = "-keep_first -drop_z_below 0 -drop_overlap", progress = TRUE)

plot(RPAS_LASCatalogue)
plot(LIDAR_LASCatalogue)

LIDAR_DTMGrid <- grid_terrain(LIDAR_LASCatalogue, 1, knnidw(), use_class = c(2L, 9L))

RPAS_LASFilterFunction = function(chunk)
{
  RPAS_LASChunk = readLAS(chunk)
  if (is.empty(RPAS_LASChunk))
    return(NULL)

  RPAS_LASChunk <- lasfilter(RPAS_LASChunk, Classification == 3L | Classification == 4L | Classification == 5L)
  RPAS_LASChunk <- filter_duplicates(RPAS_LASChunk)
  RPAS_LASChunk <- decimate_points(RPAS_LASChunk, homogenize(density = 2, res = 1, use_pulse = FALSE))
  if (is.empty(RPAS_LASChunk))
    return(NULL)

  RPAS_LASFile_Normalized <- RPAS_LASChunk - LIDAR_DTMGrid
  RPAS_LASFile_Normalized <- filter_poi(RPAS_LASFile_Normalized, Z > 1)
  return(RPAS_LASFile_Normalized)
}

opt_chunk_size(RPAS_LASCatalogue) <- 50
opt_output_files(RPAS_LASCatalogue) <- "D:/Projects2020/p2020_1021_RPAS_EFI_test/RGB/RPAS_LAS/RPAS_LASFile_Normalized_{ID}"

catalog_apply(RPAS_LASCatalogue, RPAS_LASFilterFunction)


texts: set anchor for a PictureMarkerSymbol in arcgis android sdk

texts:
I am trying to create a PictureMarkerSymbolin arcgis android, however I found that the marker use the center of the drawable as the anchor point, is it possible to set the anchor?

update:

This is the drawable(with 40x56 in pixel):

This is  my code for adding marker:

protected void onCreate(Bundle savedInstanceState) {
    final MapView mapView = (MapView) findViewById(R.id.map_imp_esri);
    mapView.addLayer(new ArcGISTiledMapServiceLayer("http://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer"));
    Point p = new Point(-9664114.480484284,3962469.970217699); //Birmingham
    mapView.setResolution(mapView.getMinResolution());

    GraphicsLayer graphicsLayer = new GraphicsLayer(mapView.getSpatialReference(), new Envelope(-180, -90, 180, 90));
    mapView.addLayer(graphicsLayer);

    //add marker
    Drawable d = getResources().getDrawable(R.drawable.marker_a_normal);
    PictureMarkerSymbol sym = new PictureMarkerSymbol(d);
    sym.setOffsetX(-convertToDp(d.getIntrinsicWidth()) / 2);
    sym.setOffsetY(-convertToDp(d.getIntrinsicHeight()) / 2);
    Graphic g = new Graphic(p, sym);
    graphicsLayer.addGraphic(g);
}

private int convertToDp(double input) {
    final float scale = getResources().getDisplayMetrics().density;
    return (int) (input * scale + 0.5f);
}

But this is what I got::

The four screenshot are captured during zooming in the map.

Note, at the image 1, the marker event does not at the united states.

Do I miss anything?


texts:
I have create a layer that is the represents the overlay areas based on some condition. The overlay is represented in purple:

I want to know the area of the overlay - which is a separate layer that has been created using the geoprocessing tool intersect. To calculate the area of the overlay I have done the following:

select the layer overlay in the layer panel -> right click on open attribute table -> activate toggle editing mode -> open field calculator -> write query as following:

The problem is that no calculations are made. In the attribute table there is only NULL for every cell entry.

Am i missing something? How can I calculate the area of overlay?


texts: QGIS use item id in print composer for relations

texts:
i'm trying to use the item id in the print composer of QGIS to build relations.
for example:
i have a label with the text "test" in it and it has the ID=label1
-> now i have a second label in my print composer and want to show the text of label 1 in this label automatically (and it should change if i change label 1)

is something like this possible?


texts:
I'm trying to divide a polygon in smaller polygons from a line, but I think I can't use the st_split function. What I need is to create small polygons inside a big one, using linestring grid.

I've tried some ways, but I can't get the result. What I've tried:

Divide one polygon from a LineString using st_split()

From a boundary polygon.

And Linestring table:

I'd need the following polygons:

Problem: I can't Split a polygon from several lines, neither a polygon from a multilinestring.

The other method I'm trying is to create a polygons from the lines with st_polygonize() The SQL I was trying is:

SELECT
   g.path[1] as gid,
   g.geom::geometry(polygon, 22033) as geom
FROM
   (SELECT
     (ST_Dump(ST_Polygonize(geom))).*
   FROM linestable
) as g;

Extracted from Creating polygons from line segments using PostgreSQL and PostGIS

Problem: I can only get one polygon (the boundary).

Can someone tell me which would be the best way to get the polygons from the linestring, or if am I missing something?

Note: Tables are in the same SRID, and geometries are snapped into a grid. In QGIS I can run the polygonize process from lines to polygon perfectly.

As John's demand, here is the linestring table.
https://drive.google.com/file/d/0B603y_m735jfS014S0EyVnpMUEU/view?usp=sharing


texts:
I'm writing an ArcGIS 10.1 .pyt toolbox to accept a feature class as input (not a shapefile).   Referencing ESRI's online guide to  'Geoprocessing Data Types of Parameters and Environments' [PDF warning], I can set the datatype to either 'DEFeatureClass' or 'IFeatureClass'.  What's the difference between these two datatypes?  Which should I use and why?

Example code from the getParameterInfo method:

treePolys = arcpy.Parameter(
    displayName = "Input tree polygons feature class:",
    name = "tpolys",
    datatype = "DEFeatureClass",
    parameterType = "Required",
    direction = "Input")


texts: I have a problem trying to convert a .shp file which I successfully save as from QGIS. Im in MapInfo trying to use universal translator to convert that shp file into TAB file but it keeps giving error as translation failed. Is there any solution or work around to my problem?


texts: Bing image services stopped working in Silverlight

texts: I have been using the Bing image and street services for several years with great success.  However today the services have gone down.  Did I miss the memo about them going away?  Or are they experiencing technical difficulties?  I have switched my app to use ESRI data instead for the time being.


texts:
I am using a Mac OS High Sierra, version 10.13.6.

While installing the more recent version of QGIS, I ran into a popup that said "Failure, The GDAL 2.3 framework is required."

Naturally, I opened up terminal, and entered tried to bash install GDAL, but it doesn't appear to be working.

sudo pip3 install GDAL

gives me

Command "python setup.py egg_info" failed with error code 1 in /private/tmp/pip-install-g7dpscpt/GDAL/

Does anyone know why this is?


texts: OpenLayers 2.12 and http basic authentication woes

texts:
I am trying to use OpenLayers 2.12 to display WMS layers from a server which has HTTP Basic Authentication enabled.

I have tried to handle the authentication by putting the username and password in the URL parameter in my JavaScript code. Example layer creation:

myLayer = new OpenLayers.Layer.WMS('background',
        'https://username:password@ws.nls.fi/rasteriaineistot/image?',
        {
            layers: 'background',
            bbox: '-380188,6249943,1347312,8226943'
        },
        {
            displayInLayerSwitcher: true,
            isBaseLayer: false,
            projection: 'EPSG:3067',
            visibility: true
        });

Of course this is not secure since the credentials are stored in JavaScript code and does not work in all browsers. Internet Explorer 8 gives security error pointing to OpenLayers.js and refuses to display the map at all. Firefox 13 pops up some authentication dialogs which I can cancel(the map displays correctly after that). In Chrome 23 the authentication seems to work flawlessly.

Can you confirm that it is not possible to handle the HTTP basic authentication in a cross-browser manner by encoding it in the URL and giving it to OpenLayers like in the example?

Can you suggest alternative ways to handle the HTTP basic authentication so that it works transparently to the user (no authentication popups displayed)? Perhaps use some kind of proxy server to work around this.


texts: How do I color XY coordinates with different colors in ArcGIS for Desktop?

texts:
I have added XY coordinates on map layers in ArcMap.

These coordinates now show up as points.

However, I want to be able to color each coordinate point a different color.

Any tips/explanation?


texts: Geoserver - Layered or Georeferenced PDF

texts: Has anyone had any experience with outputting maps from GeoServer into a Layered and / or Georeferenced PDF?


texts: Generate Street List From Centerline With "From" and "To" Street Name Attributes?

texts:
I'm trying to sort out an efficient way with which to generate a street list within my municipality that includes "from street" and "to street" street name attributes sorted by municipal ward (I've completed this portion of the ask).
Data I have:

I have joined the Ward data to the Centerline successfully. I'm just not sure about my next step.
As per the screenshots, the intersection data has a connecting ID to the Centerline.
How I translate that relationship into correctly being able to pull street names in a "from" and "to" context is where I am foggy. Within the point data, the INTERSEC_1 field has the cross streets, but I'm not sure if that alone is 100% the solution to procure what I need via a table join or some other measure.


texts: Renaming a variable input to replace spaces with underscores to be used in an output's name

texts:
This may be easier than i am making it. Using Model Builder, I have a variable input "Name", which is set as a model parameter, used for a 'Select by Attribute' early in the model. I would like to then use this variable as an input into the final feature class output's name. My issue is that some of these "Names" have spaces in them (E.G. Name = "Watershed 1"). Is there a way, using Model Builder, that i can replace the spaces with underscores/remove the spaces before inserting the variable into the output's name?

Resulting in (....spatialjoin_dissolve_Watershed_1)

Since this model will be run individually for multiple areas i need the name of the watershed to be included, without having to have the user enter in the watershed name with underscores as a parameter.

Thanks for any help!


texts: How to determine when the map has finished rendering

texts:
I have written a few functions to compare QGIS performance to MapInfo. Most of them involve doing some analysis then dumping the results on the map, as such the process would be considered 'finished' when the map has finished refreshing/rendering. I believe in this context refreshing/rendering are pretty much the same thing?

I've had a look around the documentation and have found references to (all under QgsMapCanvas class) renderStarting(), renderComplete(QPainter), mapCanvasRefreshed(). These are all 'signals', but I'm not sure how to use them (shame there are no examples in the documentation). Could anyone give me an example? Or is there a better approach? The important thing is to get an end-time when the map has finished rendering


texts: Delete row with specific text in a specific field at a specific time

texts:
I have a feature class that resides in an Enterprise geodatabase, I need to delete a row that has specific text in a specific field at a specific.
Data is being submitted daily to this feature class via Survey123, a "dummy" point is created daily for other purposes.  I would like to delete that "dummy" point daily at a specific time using ArcPy.  I am novice with ArcPy.
import arcpy

myFC = ("target dataset")

commentfield = 'Comments'
sql = """{0} LIKE '%NameOfTextToSearchFor%'""".format(arcpy.AddFieldDelimiters(myFC,commentfield))

with arcpy.da.UpdateCursor(myFC, commentfield, sql) as cursor:
    for row in cursor:
        cursor.deleteRow(row)

I couldn't figure out how to indent the last two lines


texts:
I am wishing to clip a set of polylines (black lines in image below) to the outer boundary of a polygon. Any voids within the polygon should be ignored. My ideal output are the dashed yellow lines. The initial lines may or may not be straight. The image is a simplified example, in reality the polygon is much more complex and there are hundreds of lines. I don't think a convex hull would work (but I might be wrong). I'm open to solutions in arcgis, qgis, arcpy, shapely, etc. Coding would preferably be in python by I'm open to other options if necessary. Arcgis would also be preferable to make it easier for my coworkers to share the tool but is not a requirement.

The best i can think of right now is to intersect an individual line with the polygon creating a set of points at all boundary intersections. Sort the points by distance to the start of the line. The furthest and closest (FAC) points will be the outer boundary of the polygon. Then use the FAC points to select the proper vertices from the original line and create the yellow dashed line from the appropriate points. It should work but seems more complicated than necessary.

A few additional thoughts:

The lines are linear "enough" that a simple distance calculation between points should work, linear referencing shouldn't be necessary.
This would be easy in arcpy if there was a tool to split a line at a point but I cannot find one.

Thoughts anyone?


texts: View on update event listener in the ArcGIS API for JavaScript 4.x

texts:
In the ArcGIS API for JavaScript version 3.x, there is a an event listener for ArcGISDynamicMapServiceLayer. In version 4.x of the API, I found that the MapImageLayer is almost the same as the ArcGISDynamicMapServiceLayer. To listen for update-starts and update-end in version 3, I would do add in an event listener like this:

mapImageLayer.on("update-start", function() {
   // ...
})

How would I achieve the same result in version 4? I read their API documentation but did not find anything and am wondering if this is even possible.


texts: I want to make a shapefile from an excel spreadsheet, after clicking on Add Data button and leading to the folder containing the excel file, it doesn't show any excel spreadsheet and the drop down list of "show of type" just has one option which is "Datasets, Layers and Results" .


texts:
In the gdalbuildvrt documentation, I read

  This utility is somehow equivalent to the gdal_vrtmerge.py utility [...]

So what is the difference (if any) between gdalbuildvrt and gdal_vrtmerge.py?

What does "somehow equivalent" mean here?

Is one an alias for the other?
Do they have the same command line interface?
(I.e., do they accept the same parameters and options,
and do so with the same semantics?)

(It seems so.)
Did one replace or deprecate the other?
Should one be preferred over the other?

And where can I find more information about gdal_vrtmerge.py? It doesn't seem to be included in my GDAL installation, but it's in the GDAL sources at gdal/swig/python/samples/gdal_vrtmerge.py (download | view on trac) on trunk.


texts: Setting PDF page numbers/names when created by Export Map enabled with Data Driven Pages?

texts:
how to create PDF with Grid Number in  Data Driven Pages i tried to export PDF but it is export with page no like 1,2,3..... but i want to export with Grid Number is there any way to export like that could you please tell me.

http://resources.arcgis.com/en/help/main/10.2/index.html#/Exporting_Data_Driven_Pages/00s900000038000000/


texts:
My main purpose is creating a set of users with specific roles.

I'm running GeoServer 2.13.

On the first step I wanted to create a role with access to one workspace - "UNIVERSITY"

So I've created:

A new user "test_student" in Default XML user/group service.
A new role "STUDENT" without parents.
A rule for the "STUDENT" role with a permission to write all layers and groups in workspace "UNIVERSITY".

Then I logged in the GeoServer account under the new user (test_student) and got nothing:

But I expected to get all layers in specified workspace. By the way when I set "ADMIN" as a parent for my role "STUDENT" I got all data without any restrictions.

How to solve this problem?


texts: Running ArcGIS Split By Attributes tool gives error?

texts:
While I was doing an exercise for one of my subjects, it was asked to use the "split by attribute" feature on the ModelBuilder in order to get separate shapefiles. The issue is that as soon as I open the tool, wether its inside or outside the ModelBuilder, new or existing project, it shows up an error in the "Input table" section. If I hover my mouse on the error, it wont show any message

If I try to run it, it won't output me anything


texts: Selecting top 10 records in attribute table using ArcMap

texts:
I am using ModelBuilder in ArcMap 10.2 to try and find a way to select the first 10 records sorted by Category and then by Distance. For example, my table has 100 banks and 100 churches, which is the "Category" field, and I have the distance from my site to each record stored in the "Distance" field. Essentially I want to be able to select the 10 closest banks and the 10 closest churches and export the selected records, using a tool in ModelBuilder.

Do you have any suggestions?


texts: display geojson file with folium (using geopandas to generate file)

texts:
I can not display properly a geojson file in folium using the GeoJson function. My file is made with geopandas, has some tags as fields and also a geometry field. In each row, in the field of geometry, it has a LineString (shapely object) in epsg:4326 with similar coords to the one of the center of the map of folium. The code to save my dataframe that later I read:

df.to_file('/path_to_json_file/filename.json',driver='GeoJSON')

What I tried:

import folium
import geopandas as gpd

map = folium.map.LegacyMap(location=[43.263667, -2.923202],zoom_start=14)
file = '/path_to_json_file/filename.json'
df1 = gpd.read_file(file)

#to have more info about the dataframe
print(df1.loc[0,'geometry'].coords[0:3])
>>>[(43.2614429399418, -2.927346548046099), (43.261432339487435, -2.927343485313923), (43.26142161221715, -2.927341407690571)]
print(df1.crs)
>>>{'init': 'epsg:4326'}

# 2 possibilities I tried:
folium.GeoJson(json.load(open(file)),name='geojson').add_to(map)
folium.GeoJson(file,name='geojson').add_to(map)

folium.LayerControl().add_to(map)
map.save("folium_map.html")

The thing is if I try only to display a row using folium's PolyLine instead, it works fine:

folium.PolyLine(df1.loc[0,'geometry'].coords[:]).add_to(map)

Maybe I have to declare the GeoJson in other way? I check the official documentation, but it does not appear how to pass the dataframe exactly. Checking the description of the function it seems it accepts using the path of the file like I do. The code executes without error but when I display the html it does not appear any line of the dataframe. My version of folium is v0.5.0


texts:
I'm trying to study how to find the distance between two points if their latitudes and longitudes are given.

I have a question about the proof of part of the haversine formula given at Math Forum. It says that the length of chord AD, two points at the same latitude, lat1, on a sphere of radius 1, is

2*sin(dlon/2)*cos(lat1)

but I couldn't get how they obtained it. Could you help me?


texts: USGS Earth Explorer: Bulk Download from metadata export

texts: I want to create a large bulk download order for Landsat Level-3 FSCA data. My search returns ~6700 results, but as far as I can tell I can only add one page of results (10) at a time. Is there a way to add all results to a bulk download order? I exported the metadata to csv ...


texts: Is it possible to add dynamic controls to a dockable window?

texts:
I am writing a small application that uses the dockable window. What will happen is that for each layer in the TOC a label and checkboxes will appear within the dockable window. However what do i need to do to add the control onto the Dock Window.
As an example i have

                IDataset dataset = enumDataset.Next();
            while (dataset != null)
            {
                MessageBox.Show(dataset.Name);
                //create a label object
                Label label = new Label();
                label.Name = "lbl_" + dataset.BrowseName;
                label.Text = dataset.Name;
                label.Top = top;
                label.Height = height;

But there is no

this.Controls.Add(label);

Does anyone know what the code should be to add the control?

Thanks


texts: Can't understand Return value in gaiaLineLocatePoint_r

texts:
I am trying to develop a C application that finds the nearest point on a linestring from another point, in spatialite. I found
the two API's gaiaLineLocatePoint_r and gaiaLineLocatePoint , but I am not able to understand their output. How to get the coordinates from them?

explanation of the API is here


texts:
I found map on ArcGIS Server I would like to use on my Android device:
https://experience.arcgis.com/experience/67697fc9eed74133aecb4b72dd348d12
What application I should use in order to use this map for navigation? Usually I'm using Locus Map software, but I can't find how to add such map.


texts: Lizmap - QGIS server - PostgreSQL error layer not valid

texts: I've installed Lizmap in Windows 10 but a problem is that I haven't the edition button and can't display layers from PostgreSQL. I have a "service non disponible" with "layer not valide" error in QGIS Server log when layer is added into QGIS project from PostgreSQL. Connection with PostgreSQL is working and I have authz, count... tables created in PostgreSQL after installation, but always error persist. What is the key? I installed Apache 2.41+ fcgi module, qgis-ltr 3.4 and qgis-ltr server 3.4, PostgreSQL 11 with PostGIS. Also with QGIS Server getmap url to QGIS project I have a same error "layer not valid" in Chrome.


texts: Changing Button's Enable-State on ArcMap10-Start

texts:
How can the enable-state of a AddIn-Button be influenced on start of ArcMap? My button should be disabled on startup of ArcMap and than turned enabled when the Edit-session is started.
In ArcGIS 9.3 this was quite simple:

public override bool Enabled
{
get
{
     editor = (ESRI.ArcGIS.Editor.IEditor)m_application.FindExtensionByName("ESRI Object Editor");
     if (editor.EditState == ESRI.ArcGIS.Editor.esriEditState.esriStateEditing)
     {
         return true;
     }
     return false;
}

Any suggestions?


texts: Diagnosing cause of slow performance in MXD

texts:
Every now and then I'll have an MXD that suddenly begins to not only draw at agonisingly slow speeds, but also become sluggish in general. As a routine solution, I generally run MXD doctor, remove unnecessary layers, clear results, and remove layers one by one til I identify the problem.

However, I am wondering if there is a way to quickly return diagnostics that would identify problematic layers in the first instance. The in-built 'Analyze Map...' function has not generally been helpful in these cases. Is there a way to programmatically return the drawing time of each layer - or any other creative solutions to return helpful info?


texts: pgRouting: Is it possible to do multiple (concurrent) Dijkstra/A*/A_shooting_star queries at a time?

texts:
Disclaimer: the author of the question is a novice in PostgreSQL and pgRouting.

As I understand doing a shortest path query in pgRouting is a complicated SELECT underneath.

And as I know unlike any UPDATE queries, SELECT queries in PostgreSQL can be done simultaneously, because they don't interfere with the the database integrity anyhow.

If the above is correct, the question is:

Can I query pgRouting for shortest paths, for example, each 1 millisecond and get the same timing results as in case with truly sequential shortest path queries? Is case 2 possible?

Usual case (1):

[Time elapsed: 0 seconds]
Query for shortest path A
Get result for the shortest path A
[Time elapsed: 1.5 seconds]
Query for shortest path B
Get result for the shortest path B
[Time elapsed: 3 seconds]

A case with simultaneous queries (2):

[Time elapsed: 0 seconds]
Query for shortest path A
[Time elapsed: 0.001 seconds]
Query for shortest path B
Get result for the shortest path A
[Time elapsed: 1.5 seconds]
Get result for the shortest path B
[Time elapsed: 1.5001 seconds]


texts: Can the tools in QGIS be used in python code?

texts: Can the tools in QGIS, e.g. Vector>>>Geometry Tool>>>Extract Vertices be used in Python code? For example, to run through a shapefile to get all the points for  for polygons within the file?


texts: ArcPy: Select_Analysis Invalid Characters

texts:
The following script produces the ERROR 000210: Cannot create output... .shp and ERROR 000354: The name containts invalic characters Failed to execute (Select.).

The line with arcpy.Select_analysis is causing the problem

The result should be shapefiles with the SAME NAME as the Input, but in the Output-Workspace-Folder!

Is there a problem with the path or would you have another idea?

import arcpy, os, math
from arcpy.sa import*
arcpy.CheckOutExtension("Spatial")

pathworkspace = r"D:\aprettylongpath\InputWithShapefilesInside"
outputworkspace = r"D:\aprettylongpath\OutputStillEmpty"

arcpy.env.workspace = pathworkspace
mylist = arcpy.ListFeatureClasses()
where_clause = '"GRIDCODE" =2'

for z in mylist:
    print "Shapefile..."+z
    inShapefiles = z
    arcpy.Select_analysis(inShapefiles,os.path.join(outputworkspace,z),where_clause)

The result from the print in the 3rd last line is:

Shapefile...1111111-rgbi-10cm_compband.shp


texts: How to add a raster layer to Geoserver via the REST-API?

texts:
I want to automatically (via script) add new TIF-files on the harddrive to Geoserver as a new layer. I use the REST-API for that.

In the first step I added the TIF as a coveragestore, and that step works. So I use this JSON:

{
    "coverageStore":
    {
        "name": "test123",
        "type": "GeoTIFF",
        "enabled": true,
        "workspace":
        {
            "name": "TEST"
        },
        "url": "file:\/\/\/data\/test123\/test123.tif"
    }
}

I add that via wget and that part works:

wget --user=test '--password=test123' --header='Content-Type: application/json' --post-file=coveragestore-test123.json 'http://localhost:8080/geoserver/rest/workspaces/TEST/coveragestores'

So far it s fine. But now I want to publish that as a layer, and that doesn't work. My JSON here looks like this:

{
    "layer":
    {
        "name": "test123",
        "type": "RASTER",
        "defaultStyle":
        {
            "name": "TEST:style123",
            "workspace": "TEST"
        },
        "resource":
        {
            "@class": "coverage",
            "name": "TEST:test123"
        }
    }
}

And I send it with this:

wget --user=test '--password=test123' --header='Content-Type: application/json' --post-file=layer-test123.json 'http://localhost:8080/geoserver/rest/layers'

It comes back with an HTTP-error 405.

What do I do wrong?


texts: I have some lat long coordinates but do not have any other information (datum, CRS etc) . I know where the place on Earth is but I cannot locate it (is a few miles off) using the coordinate numbers I have using for example Google Maps. One coordinate is (9.22618  -70.25967)


texts: How is it possible that Geoserver WMS map works perfectly in GIS software but the same map cached using integrated GeoWebCache using the same coordinate system gives the same map in other location. I tried different data with EPSG:3059, but always the same offset.  Where should I look for a misconfiguration or what could be the problem?


texts: How to get (x,y) coordinates given MinGX, MinGY, MaxGX, MaxGY?

texts: I have a set of data that stores the MinGX, MinGY, MaxGX, MaxGY values of each location. I need to find the (x,y) coordinates of these locations. Is there any way to do this without using/downloading any software? Also, is it possible to plot the map given these 4 values by using any software?


texts: Showing unknown SRS (native SRS) when importing shapefile into GeoServer

texts:
I am facing a problem on an unknown SRS when I am importing a shapefile into GeoServer. I am trying to import a shapefile into GeoServer but when I am publishing the shapefile, it is giving the unknown SRS in Native SRS on GeoServer. I am giving the link of screenshot of issue: ( http://mindcrewgis.com/screensho_error/shape_file.png)

Then after I set the SRS (900913) and publish the layer, I am getting another error on the layer preview when I click on any feature on the layer. It is giving an error like:

  java.lang.RuntimeException: org.geotools.data.DataSourceException: An exception occurred while parsing WKB data org.geotools.data.DataSourceException: An exception occurred while parsing WKB data An exception occurred while parsing WKB data null

And it does not give any feature information.
The screenshot is here: http://mindcrewgis.com/screensho_error/feature_error.png

But the same shapefile works fine on QGIS on my Desktop.
Please give me the solution to the above problem.


texts: Ideas for showing Land Use and Zoning on a Single Map

texts:
I was wondering if any of you had any creative ways to clearly show zoning and landuse together on a single map.  I thought about showing one in a thick outline and the other in a fill with no outline but wondered if any of you had any better ideas.

Does anyone have a better idea?


texts: Error while running batch processing in QGIS 3

texts:
I have built a straightforward model in the Graphical Modeller in QGIS 3.4 to select features by expression and then save the selected features to a new shapefile.

It works when running singly but fails when run as a batch process with the following message:

The same model in QGIS 2.14 works when run both individually and as a batch process.

I have seen similar questions QGIS- Model works on single layer execution, errors out when run as batch but no explanation or indication how to get around this in QGIS 3 (apart from using 2.14 or writing a Python script).

Is this a known issue? Or is there a fix to get it to work in QGIS 3?

I am trying to demonstrate this feature to students - my current plan is to use QGIS 2.14 for this part of the workshop.


texts: Group and aggregate SpatialPolygonsDataFrame by column in R

texts:
I'm working with rgdal, sp, maptools, etc.

I have a SpatialPolygonsDataFrame that has each row/Polygons list with only one Polygon.  It also has an identifier column which I want to use to group by, and other attribute columns.

I would like to get the equivalent of group_by, where the attribute columns are aggregated, but where the corresponding Polygons' Polygon are concatenated in the same list.

Is there a function to do this?


texts: Importing shapefile to file geodatabase feature dataset using ArcObjects and C#?

texts:
I'm trying to import a shapefile into a featuredataset under a geodatabase using FeatureDataConverter.ConvertFeatureClass. All the examples that I found are not concerning the import under an existing featuredataset.
I checked that the spatial reference it's the same and the domain it's enough big to contain the features. When I try to run stop at featureDataConverter.ConvertFeatureClass and exit with out any information.

...
ISpatialReferenceFactory srFactory = new SpatialReferenceEnvironmentClass();
                Int32 EPGS = 25832;
                ISpatialReference spatialReference = (srFactory as ISpatialReferenceFactory4).CreateSpatialReference((int) EPGS);
                IEnvelope envelope = RvrFtrLyr.AreaOfInterest;
                spatialReference.SetDomain(envelope.XMin-1000, envelope.XMax+1000, envelope.YMin-1000, envelope.YMax+1000);
                IFeatureDataset targetFeatureDataset = coreEngine.CreateFeatureDataset(targetWS, "Rivers", spatialReference);

                coreEngine.FeatureLoader(targetFeatureDataset, "Lines", RvrFtrClsSHP);
...

internal void FeatureLoader(IFeatureDataset targetFDS, String name, IFeatureClass sourceFeatureClass)
        {
            //Source feature
            IDataset sourceDataset = (IDataset)sourceFeatureClass;
            IWorkspace sourceWS = (IWorkspace)sourceDataset.Workspace;
            IFeatureClassName sourceFtrClsName = (IFeatureClassName)sourceDataset.FullName;

            //Target feature
            ESRI.ArcGIS.esriSystem.UIDClass clsid = new ESRI.ArcGIS.esriSystem.UIDClass();
            clsid.Value= "esriGeoDatabase.Feature";
            IFeatureClass ftrCl = targetFDS.CreateFeatureClass(name, sourceFeatureClass.Fields, clsid, null, esriFeatureType.esriFTSimple, "Shape", null);

            IDataset dataset = (IDataset)ftrCl;
            IFeatureClassName trgFCN = (IFeatureClassName)dataset.FullName;
            IFeatureDatasetName trgFDN = (IFeatureDatasetName)trgFCN.FeatureDatasetName;

            IFieldChecker fieldChecker = new FieldCheckerClass();
            IFields sourceFields = sourceFeatureClass.Fields;
            IFields targetFields = null;
            IEnumFieldError enumFieldError = null;

            fieldChecker.InputWorkspace = sourceWS;
            fieldChecker.ValidateWorkspace = targetFDS.Workspace;

            fieldChecker.Validate(sourceFields, out enumFieldError, out targetFields);
            if (enumFieldError != null)
            {
            Console.WriteLine("Error in fields validation");
            }

            String shapeFieldName = sourceFeatureClass.ShapeFieldName;
            Int32 shapeFieldIndex = sourceFeatureClass.FindField(shapeFieldName);
            IField shapeField = sourceFields.get_Field(shapeFieldIndex);

            IGeometryDef geometryDef = shapeField.GeometryDef;
            IClone geometryDefClone = (IClone)geometryDef;
            IClone targetGeometryDefClone = geometryDefClone.Clone();
            IGeometryDef targetGeometryDef = (IGeometryDef)targetGeometryDefClone;

            IFeatureDataConverter featureDataConverter = new FeatureDataConverterClass();
            IEnumInvalidObject enumInvalidObject = featureDataConverter.ConvertFeatureClass(sourceFtrClsName, null, trgFDN, trgFCN, targetGeometryDef, targetFields, "", 1000, 0);

            IInvalidObjectInfo invalidObjectInfo = null;
            enumInvalidObject.Reset();
            while ((invalidObjectInfo = enumInvalidObject.Next()) != null)
            {
                // Handle the errors in a way appropriate to the application.
                Console.WriteLine("Errors occurred for the following feature: {0}",
                    invalidObjectInfo.InvalidObjectID);
            }
        }


texts: OpenLayers3 change Layer source URL (or replace features loaded from another URL)

texts:
I have a map with a Vector layer, containing features from a GeoJSON source:

var map = new ol.Map({
    layers: [new ol.layer.Tile({
                 source: new ol.source.OSM()
             }),
             new ol.layer.Vector({
                 source: new ol.source.Vector({
                     url: 'http://example.com:5000/geo/data/zones/1',
                     format: new ol.format.GeoJSON()
                 }),
             })],
    renderer: 'canvas',
    target: 'map',
    view: new ol.View({
        center: [737514.438475665, 5864533.629390752],
        zoom: 13
    })
});

I have multiple URLs that returns s GeoJSON string:

http://example.com:5000/geo/data/zones/1
http://example.com:5000/geo/data/zones/2
http://example.com:5000/geo/data/zones/n

I need to be able to switch the URL of my Layer source (or fetch and display features from another URL).

I have tried to find the 'url' property on the ol.Layer.Vector instance:

l=m.getLayers().getArray()[1]
l.getProperties()

and on the ol.source.Vector instance:

s = l.getSource()
s.getProperties()

but I haven't found anything about the 'url' property.

Could you provide a way to do that ?

is it possible to simply update the source URL (and automagically refresh the layer features) ?
shall I remove existing features, load new features using my own logic and add the loaded ones ?
shall I remove the whole Layer, re-create it, and re-add it ?


texts: Is it possible to connect Maptitude 5 with Oracle Spatial?

texts:
I have a version of Maptitude 5, and I need to connect it to Oracle Spatial.
I have the client installed on the machine, including the connection to the database using Qgis and gvSIG and this works normally.
Using Maptitude:

Open File Menu
Files of type: Oracle Spatial Layer
In the window that I've opened, entered User, Password and Service

This then gave an error that the connection was not performed with these parameters, but these are the parameters used in Qgis and gvSIG.

I know this is an old version of Maptitude, but I wonder if anyone uses or has used this version connected to Oracle Spatial.


texts: How to specify alpha channel selection inside SLD raster style on GeoServer

texts:
Is there a way to specify alpha channel inside raster SLD in GeoServer.
For example, I have a multispectral map with 6 bands(B-G-R-RE-NIR-Alpha), where the sixth band is alpha, and I would like to show it as RGBA using SLD style.
I can order channels in the ChannelSelection tag to show RGB without the Alpha channel, but how can I add an alpha channel here?

<ChannelSelection>
  <RedChannel>
    <SourceChannelName>3</SourceChannelName>
  </RedChannel>
  <GreenChannel>
    <SourceChannelName>2</SourceChannelName>
  </GreenChannel>
  <BlueChannel>
    <SourceChannelName>1</SourceChannelName>
  </BlueChannel>
</ChannelSelection>


texts: How to plot Circle of radius 1km on Google Maps

texts:
I have circle of radius 1km in QGIS (centered at: -122.69690,58.80688 WGS84, North America), when I export it to Google map it isn't 1km in radius (I checked using the measure tool on Google maps) and the shape isn't a perfect circle as seen below.

On Google maps, North and South direction from the center are more that 1km radius while in East and West direction it is less than 1km.

How can this be fixed?

In QGIS

In Google Maps


texts:
I'm attempting to map cell phone towers and their range in QGIS. I'm wondering if it's possible to create buffers which change color intensity according to their distance from a central point. The buffer color is most intense at the center and gradually becomes lighter as one moves away from the center.

I speculate that I must copy the range of each tower to the buffer layer attribute table and then select "Graduated" in the style menu of the buffer. Unfortunately, the buffer layer has only 1 row in its attribute table.

Are there other solutions?


texts: Splicing Large Geometry into Smaller Geometries for Export in Google Earth Engine

texts:
I've been working on a task for a few days that involved iterating geometries and exporting time series data within these geometries. I want to be able to create a 1000 km^2 tile and then produce 100 tiles that are 100km^2 within this larger tile. Within the smaller tiles I want to export NDVI time series data to my Google Drive.

The desired output is 100 .tif file exports for a large area that can then be further analyzed in a different data processor.

This is the code that I am working with ("point" can be defined arbitrarily):

var targetArea = ee.Number(1e12)

function toSquare(buffer) {
  return point.buffer(buffer, 1).bounds()
}

function getBuffer(area) {
  return ee.Number(area).sqrt().divide(2)
}

var buffer = getBuffer(targetArea)
var minBuffer = buffer.multiply(0.99)
var maxBuffer = buffer.multiply(1.01)
var square, area
for (var i = 0; i < 100; i++) {
  var buffer = minBuffer.add(maxBuffer).divide(2)
  square = toSquare(buffer)
  area = square.area(1)
  minBuffer = ee.Number(ee.Algorithms.If(area.lt(targetArea), buffer, minBuffer))
  maxBuffer = ee.Number(ee.Algorithms.If(area.gt(targetArea), buffer, maxBuffer))
}

print(area, 'area')
var perimeter = square.perimeter(.001);
print(perimeter, 'perimeter');
Map.addLayer(square)
Map.centerObject(square)

This is the link to the script pasted above and it creates a roughly 1km^2 tile with the data that I want to be exported:
https://code.earthengine.google.com/c48fb6bf2bbc8bc0511d95ab9c360fdb

The obvious problem is that the export would be too large for Earth Engine to handle at higher quality. Is there a way to break this down into smaller tiles through some sort of loop? I want the tiles to touch, so that no data is left out within the large square.


texts: Remove Yellow Squares Created by the LRS Plug-in

texts:
Does anyone know how I can remove the yellow squares which appear in my QGIS project when using the "Locate" tool within the "LRS" Plug-in?
I have refreshed the Project and also can't deselect them! They are also apparent on the screen if I turn all the layers off.


texts: How to do different thematic analysis to same GeoJSON in Leaflet

texts:
In order to do several Thematic Analysis to my GeoJSON Polygon file, I have to load each time the duplicate version of it, in order to have that color thematic of the specified attribute (This is not good idea).
I want to use only one GeoJSON file & do several thematics.

var mymap = L.map('mapid').setView([35.402484, 8.125763], 13);
var OpenStreetMap_Mapnik = L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 19,
    attribution: '&copy; OpenStreetMap contributors'
}).addTo(mymap);

//Color class for Population
var populationtClass = [
[0, 'rgb(249, 253, 255)'],
[500, 'rgb(195, 225, 243)'],
[1000, 'rgb(180, 219, 243)'],
[2000, 'rgb(124, 188, 227)'],
[5000, 'rgb(89, 160, 203)'],
[9000, 'rgb(30, 123, 180)'],
[19800, 'rgb(8,  86, 133)']
 ];
 //Color Class for Job
 var JobClass = [
[0, 'rgb(255, 255, 255)'],
[20, 'rgb(243, 250, 213)'],
[300, 'rgb(241, 252, 197)'],
[800, 'rgb(234, 252, 164)'],
[1000, 'rgb(224, 247, 132)'],
[2000, 'rgb(214, 250, 68)'],
[3000, 'rgb(180,222,15)']
 ];

 //function Returning the Color
 var coloring = function (thresholds, value) {
                for (var i = 0; i < thresholds.length - 1; i++) {
                    if (value < thresholds[i + 1][0])
                        return thresholds[i][1];
                }
                return thresholds[thresholds.length - 1][1];
};

//Loading Geojson of Zone1 and color it by population value
var Zone1=new L.geoJson(Zone1,{
      style: function (feature) {
              var population = feature.properties.Pop_2015;
                                return {
                 fillColor: coloring(populationtClass, population),
                        fillOpacity: 0.75,
                        weight: 1,
                        color: 'black'
                    };
                },
onEachFeature: function (feature, layer) {
        layer.on({
                        'mousemove': function (e) {
                            e.target.setStyle({
                                weight: 7,
                                color: 'red'

                            });

                          },
                        'mouseout': function (e) {
                            Zone1.resetStyle(e.target);

                        },
                        'click': function (e) {
                            mymap.fitBounds(e.target.getBounds());

                        }});}});

//Loading Geojson of Zone2 and color it by Jobs value
var Zone2=new L.geoJson(Zone2,{
      style: function (feature) {
              var JOBS = feature.properties.Job_2015;
                                return {
                 fillColor: coloring(JobClass, JOBS),
                        fillOpacity: 0.75,
                        weight: 1,
                        color: 'black'
                    };
                },
onEachFeature: function (feature, layer) {
        layer.on({
                        'mousemove': function (e) {
                            e.target.setStyle({
                                weight: 7,
                                color: 'red'

                            });

                          },
                        'mouseout': function (e) {
                            Zone2.resetStyle(e.target);

                        },
                        'click': function (e) {
                            mymap.fitBounds(e.target.getBounds());

                        }});}});

//I dont want to load twice my Geojson files in order to get different thematic
//I want to use Only one Geojson file & do two thematics once for poupulation & other for Job

L.control.layers({
'BaseMap' :OpenStreetMap_Mapnik
    },{
'Population':Zone1,
'Jobs':Zone2
}).addTo(mymap);

Here is my example running


texts: AGOL Dashboard - Select Attribute from List

texts:
I've compiled a database of proposals and geolocated them, while giving each entry / feature a "tag" of 1-3 values, separated by commas. In AGOL's Dashboard's Category Selector, I can only create a drop down menu for the entire string of the values. For example:
I have: "Bridge, Roads, Sings" as a single field. I want to be able to select "Bridge" and have all features that contain "Bridge" be selected. Is that possible?


texts: 3D model in ArcScene only with coordinates possible

texts:
How do I create a 3D model in ArcScene with only a .txt file existing?
I have the coordinates and need to put in in a 3d model.


texts: Spatial Join not adding votes correctly when mapping election precinct data by County using ArcMap

texts:
I was given precinct level election data but need to make a map at the county level. The precinct data does not have any Counties listed in the attribute table, so I took a county layer and did a spatial join with my precinct data. I did Target feature: County and Join feature: Precinct. Join_One_To_One. Selected the "County_Name" feature to join feature. Set my merge rule to "join." The new attribute table for this join shows the "vote totals" per county as completely wrong - it seems like it's only bringing in one precinct per county and not adding all of them together. For example, Los Angeles County shows 700 votes when it should be hundreds of thousands. I tried doing a spatial join the other way (Target:Precinct and Join Feature: County, and the vote totals seemed correct! But then, I'm not sure how to make a County level map when the view is all individual precincts.
Can anyone help me with this?


texts: Apend child to a layertree with nodeType: "gx_layer"

texts:
Im building an application with geoext and opengeo suite SDK.
I´m trying to create a layer tree with some "grouped" layers like the last node here
I saw somewhere that I should set the node Type to "gx_layer" to get that effect, but if I do that then the entire tree won’t load. I only get the root node and when I click on it , the icon changes to the “loading” image,  but then nothing happens


texts: Finding overlapping annotations using ArcGIS Desktop?

texts:
I use about 10 annotations to produce my map. Every single one has been created separately using Tiled Labels to Annotation. I am using ArcGIS Desktop 10.3

Is there a way to find the annotations that overlap each other?


texts:
I'm working in ArcGIS 10.3. I have a shapefile, called chunk6, that looks flattened but can be viewed normally when merged with another similar shapefile (called chunk5).

In its own data frame, chunk6 it looks like this

In the data frame of chunk5, I can't see chunk6 unless it is merged with chunk 5. Below I can only see chunk5 even though chunk6 is also selected

Here I have merged chunk5 and chunk6 which produced the shapefile called chunk6_Merge. In the merged shapefile I can now see both chunk6 and chunk5 in the chunk6_Merge shapefile

Any idea why chunk6 is flattened when viewed by itself?

Note that both shapefiles have a a defined geographic coordinate system (GCS_South_American_1969) and are not projected.


texts: Performing overlapping zonal statistics in ArcGIS?

texts:
I have buffer lines overlapping and within these lines I want to know the type of land use. However, overlapping is not allowed with Zonal statistics.

How do I do this in ArcGIS because I want all the data and it's a big dataset. Also I don't know coding (Python).


texts: Create SQL View in GeoPackage in GeoServer

texts:
In GeoServer 2.15.2 I can "publish" an entire layer from a GeoPackage.  No issues there, works great.

However, If I attempt to "Configure a new SQL view..." using that same layer in the GeoPackage (or any layer from any GeoPackage), GeoServer fails to recognize "geom" as the geometry column.  In the "Edit SQL view" screen, the "type" for "geom" is listed as "BigDecimal".

What I would like to do via the "SQL View" is to publish a layer from a GeoPackage and use a SQL query to scope the results (cannot do this if I just publish the entire GeoPackage layer).

Note:  As I mentioned, GeoServer does indeed recognize "geom" when you "publish" the entire GeoPackage layer . . . it can successfully identify the geometry type and SRID.  This fails when trying to "Configure a new SQL view...".

Update: GeoServer does not prevent me from creating/saving a SQL View layer using a GeoPackage as the "Store".  However, it is useless because the geometry type and SRID are unknown (and the Bounding Boxes are populated with invalid data).  However, if I manually edit the "featuretype.xml" file for the layer and populate it with the correct values, I can at least view it via "Layer Preview".


texts: Is it possible to calculate a value based on values of another row?

texts:
I have a point shape file in qgis. For one of the attributes, I would like to offset the values, i.e. the value for point 1 goes to point 10, the value for point 2 goes to point 11 and so on. I see the $rownum attribute in the field-calculator, but I cannot see any examples on how (if possible at all) I can use that to access values in another row...

Is this possible using the field calculator or do I have to make a python script?

(Yes, I know for the example, rows 1 - 9 will get null- values and the last 9 values will "fall off the end")

(Using qgis 2.0.2)


texts: Create new point feature of a shapefile if a condition is met when creating a new line feature

texts:
I am working with multiple ESRI Shapefiles, both Line and Point types. Every shapefile has a set of fields and some of them share common values (by default) or have values that refer to each other from a logical point of view. I need the shapefiles to be like they are so I cannot alter or join fields on them to achieve the result.
Basically, I am trying to find a way to create a New Point Feature of the shapefile STOCK, compiling all its fields in the process (ID, CABLE, LENGTH, CITY, COORDINATES), if a condition is met when creating a New Line Feature in the shapefile CABLES.
The exact example of pseudo-code is
IF(CABLE.LENGTH > 100, STOCK.CREATE_POINT(ID, CABLE.ID, 10, X.CITY, CABLE.COORDINATES))

Where CABLE.COORDINATES is exactly the coordinate of the last node of the line geometry. I need it to be this so that when I finish creating a line and it is longer than 100, a new point feature with its label and symbol already set in the STOCK shapefile appears at the end of the line and I can then move it to the appropriate final position.
I am using QGIS 10.10 for compatibility reasons of the base-project. I am trying to upgrade it with new automated features to make the job faster but I am new to gis. I will then apply more complicated boundaries based on multiple fields to get the exact STOCK Length I need.


texts: Editing PostGIS topology in a web browser

texts:
I am trying to conceptualize the best process for editing topology from PostGIS in a web browser. Ideally what I would like to accomplish is the loading of a topology (imagine the lower 48 states of the US), adjusting the borders, and updating a PostGIS table with the changes. I am fairly comfortable with the delivery, representation, and uploading of the geometry, but I am a bit stuck on the editing part.

Although it seems like OpenLayers 3 supports topology editing (despite this answer), it does not seem to enforce any rules (for example, it allows the user to drag a boundary to overlap another polygon), nor does it seem like OpenLayers supports exporting vector layers as a GeoJSON, though I'm sure it's possible to write that functionality. Likewise, Leaflet offers the excellent Leaflet.draw plugin along with GeoJSON export functionality, but does not offer the ability to edit a topology.

Unless I'm missing something, it seems like my best bet at this point is to write a tool that allows for the topological editing of polygons, most likely for Leaflet because I am more comfortable with it. I imagine the process would be something like:

Click on a boundary
Find all features that have a geometry at this location
Create a vertex at this location in all found feature geometries
While dragging, update all features that share this vertex
When done dragging, make sure that no segments of the altered geometries intersect other features

Ideally it would function very similarly to the OpenLayers example above, but would also enforce some rules and make exporting easy. Does this sound like a reasonable process? Am I missing anything? Do you have any specific suggestions for any part of the process, i.e. any existing tools that accomplish these pieces of functionality?


texts: LineSymbols not rendering in ArcGIS Runtime for WPF SDK 10.1.1 without AcceleratedDisplay

texts:
I've posted this over on ESRI's ArcGIS Runmtime for WPF forums but I thought I'd repost here to get some more developer exposure. Any input would be greatly appreciated. Original post here: http://forums.arcgis.com/threads/82721-LineSymbol-and-SimpleLineSymbol-not-rendering-without-AcceleratedDisplay

With the AcceleratedDisplay enabled, SimpleLineSymbols and LineSymbols both render correctly. If I turn it off, both fail to render completely. Nothing else has changed such as their geometries or symbologies. Is this a known issue and/or am I doing something wrong? Additionally, SimpleFillSymbols, PictureMarkerSymbols, and TextSymbols seem to work just fine regardless of AcceleratedDisplay settings. That is to say, if I render the symbols as SimpleFillSymbols they will show up (BorderBrush with no FillBrush).

Code Behind:

    private void MakeSimpleSymbol()
{
            //Create a simple line symbol
            SimpleLineSymbol sls = new SimpleLineSymbol();
            sls.Color = new SolidColorBrush(Colors.Blue);
            sls.Style = SimpleLineSymbol.LineStyle.Solid;
            sls.Width = 2;

            //Create a polygon object
            Polygon p = new Polygon();
            p.SpatialReference = new SpatialReference(102100);

            //Populate it with some points
            ESRI.ArcGIS.Client.Geometry.PointCollection pointCollection = new ESRI.ArcGIS.Client.Geometry.PointCollection();
            ObservableCollection<ESRI.ArcGIS.Client.Geometry.PointCollection> rings = new ObservableCollection<ESRI.ArcGIS.Client.Geometry.PointCollection>();
            MapPoint wgsPoint = new MapPoint(0,0);
            MapPoint wgsPoint2 = new MapPoint(50,0);
            MapPoint wgsPoint3 = new MapPoint(0,50);
            pointCollection.Add(wgsPoint);
            pointCollection.Add(wgsPoint2);
            pointCollection.Add(wgsPoint3);
            rings.Add(pointCollection);
            p.Rings = rings;

            //Create a graphic and set its properties
            Graphic graphic = new Graphic();
            graphic.Geometry = p;
            graphic.Symbol = sls;

            //Create a graphics layer and add the graphic to it
            GraphicsLayer gLayer = new GraphicsLayer();
            gLayer.Graphics.Add(graphic);

            //Add the graphics layer to the map
            _map.Layers.Add(gLayer);
}

XAML:

<Window x:Class="ESRIWPFTestApplication.MainWindow"
     xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
        xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
        xmlns:esri="http://schemas.esri.com/arcgis/client/2009"
        Title="MainWindow" Height="350" Width="525">
    <Grid>
        <!-- AcceleratedDisplay is OFF -->
        <esri:Map x:Name="_map" UseAcceleratedDisplay="False">
        </esri:Map>
    </Grid>
</Window>


texts: Adding Lambert coordinate system to QGIS

texts:
I want to define Lambert coordinate system for Iran in QGIS.

What format it should be? wkt or proj string?
How can I get the parameters I need?
Which Lambert is suitable for Iran?


texts: Does ArcMap 10 have a command line switch that is roughly equivalent to Visual Studio's /SafeMode?

texts:
Anytime I had an  ArcMap 9.x related incident open with ESRI, they always had me do the usual steps of uninstalling extensions and tools, renaming the ESRI profile folder to ESRI_old, and renaming the ESRI registry key under HKCU to ESRI_old.

I am unable to find a list of arcmap.exe command line switches, but I feel like a /SafeMode type of switch would really speed up troubleshooting.  The only command line switch I know of is /log.

Does anyone know of a /SafeMode equivalent or have a link to arcmap.exe command line switches?


texts: Refreshing map after changing Feature's style in OpenLayers 3

texts:
I have an OpenLayers 3.2.0 map that features some vector sources (ol.source.Vector) and associated vector layers (ol.layer.Vector)
When Features (ol.Feature) are added to the vector sources, they are given a data property that is set to the javascript object that the feature represents. TypeScript follows...
vectorSource.addFeature(new ol.Feature({
    geometry: /* ... */,
    data: vectorData,
}));

The vector layers then have a style-function that reads the data property and retrieves its style:
vectorLayer = new ol.layer.Vector({
    source: vectorSource,
    renderBuffer: /* ... */,
    style: function (feature: ol.Feature, resolution: any) {
        var data = </* TypeScript Type */>feature.get('data');
        if ((data) && (data.style)) {
            return [data.style];
        }
        else {
            /* return default style */
        }
    }
});

Sometimes, events unrelated to the map cause the styles to change. For example, when an object becomes invalid, its style changes. Clearly, since data.style is entirely within my control, changing it is trivial.
The problem is that the map does not know that the style has changed. If I change an object's style and then zoom the map, forcing it to redraw, I notice that my style-functions run and return the new style and the feature is redrawn. How do I programatically force the map to refresh?
After some searching and experimentation, I have tried:

Calling render() on the ol.Map itself.
Calling dispatchChangeEvent() on the ol.source.Vector
Calling redraw() on the ol.layer.Vector

These were suggested but none of them worked, which isn't surprising since only the first method is even listed in the OpenLayers 3.2.0 API documentation and it is not marked as stable.


texts: Deleting a standalone table from the TOC-- C# ArcObjects

texts:
I'm looking to delete a standalone table that has been added to the arcmap document from an access database.  I only want to delete it from the TOC, not from the file and I can't seem to get it to work.

Here is my current code:

        IMxDocument mxDocument;
        mxDocument = (IMxDocument)ArcMap.Application.Document;

        mxDocument.ActiveView.Refresh();

        ////IWorkspaceFactory workspaceFactory;
        ////workspaceFactory = new AccessWorkspaceFactory();
        IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
        IMap map = mxdoc.FocusMap;

        // Get the table
        IStandaloneTableCollection stTableCollection = map as IStandaloneTableCollection;
        IStandaloneTable standaloneTable = null;
        ITable table = null;

        IFields fields;
        fields = table.Fields;

        ITableCollection tableColl = map as ITableCollection;

        table = (ITable)map.Layer[0];

        int x = map.LayerCount;
        for (int i = 0; (i < x); i++)
        {
            if (map.get_Layer(i).Name == name)
            {
                table = (ITable)map.get_Layer(i);
            }
        }

        mxDocument.FocusMap.DeleteLayer(table);
        mxDocument.UpdateContents();
        mxDocument.ActiveView.Refresh();


texts: Adding values to point by its location in vector layer using QGIS?

texts:
I have two layers. One is the location of houses, the other is the boundaries of areas in my home town.

I'd like to add the name of those areas to the points, so I can select them, group them and display them by the area they are located in. Like "240 of those points are in the area of xyz"

I did try and select them by the intersect tool but had no success with that.


texts: Not updating the row number properly. Attributing the last created point with "1" & vice versa

texts:
We have to capture the sample points in survey maps. While doing so, we are facing a peculiar problem. Our survey points start from left to right <- along the survey lines. Once we have captured all the points, say 86 points and update the ID filed using $rownum-1 (to include zero), what QGIS does is it takes the last created point as the first and attributing it with ID = 0 (which should be 85) and the first created point as ID = 85 (which should be 0) . I've attached an image to explain the same.

1) $ID is not that much effective for our line of work. Coz at times we might miss some points and delete some and to update points with $ID is creating confusion.


texts: Merging polygons in one layer by boundaries of polygons from another layer in QGIS

texts:
So I have two polygon layers, polygons in one layer are greater than in the other one. In some cases it is possible to merge smaller polygons so that they fit within the boundaries of a bigger polygon.
Is there a way to do it in QGIS automatically and attribute the ID number of a bigger polygon to all smaller polygons falling within its boundaries?
Trying spatial join did not result: when using Spatial join predicate "within" Qgis only does operation for polygons that do not intersect with any other bigger polygon (see highlighted in yellow on green), and gives NULL for small polygons that are still within one bigger polygon but intersect another bigger polygon.

Using spatial join "intersect" results in duplicates and triplicates of smaller polygons with ID numbers of ALL polygons they intersect with. So, all the polygons highlighted in yellow on violet background will receive ID number of the polygon they belong to, but then will be duplicated with ID number of a bigger polygon they border with.
I want to get the smaller polygons with the ID number of only one bigger polygon, the one they belong to, so that I will be able to merge them based on that ID number. I do not want to get duplicates and triplicates of the same features with IDs from bordering bigger polygons.


texts: QGIS behaviour when dealing with spatialite layers (splitting/merging features)

texts:
I'm wondering what is the behaviour of qgis when you work on spatialite layers.

For instance what I've discovered so far is that when you split a geometry in two. One remain the original (same id, same attributes) and the other is a duplicata (negative id but still same attributes) that can be modified before the commit.
After you save changes there is an insert for the new feature but is there an update for the original feature because it's geometry changed ?

And furthermore, when you merge features i've noticed that no DELETE is emitted to the database.

I'll be glad if someone points me toward some doc about it.


texts: Split Leaflet layer control into multiple tabs panel

texts:
I have a Leaflet map with several layers (21) grouped into 7 categories.

Because of the large number of layers, I would like to give up on the standard L.control.layers look and turn it into a Bootstrap panel with tabs, based on my categories.

So I currently have this:

And I want to achieve this:

Is is possible to split the Leaflet layer control into a multiple tabs panel?

I have added my code to a JSFiddle. I only used some markers, for simplicity.

-- Edit --

I have tried using Leaflet.StyledLayerControl plugin to group my layers and move them to each tab based on their category.

Result: here.

Problem: After moving the elements, their styling and the click event are lost.


texts: Loading data and symbol them to the map dynamically

texts:
I am using ArcGIS SDK for Android 10.2.2, and we want to add some geographic data(graphic) to the map.

While our data is deployed by a common web service(no geographic web service like wms or wfs) which can be queried by a bbox.

Now we want to add the queried point to the map after the map zoom,pan or something else which make the map update.

But it seems that the only map event we can use is OnZoomListener and OnPanListener, but we found that we the view(no matter center or the resolution) of the map is changed by program like this:

mapView.setResolution(...); // or
mapView.centerAt(...);

These listeners will not be triggered.

So I tried to make the points a kind of Layer, because I think the layer must be sensitive to the map update event, but I have see almost all the api about the layer, I can not found a class to extend.

In fact, the ArcgisFeatureLayer is the most close to our requirement, it seems that it can automatically load new data according to the map view, but we do not deploy our data by arcgis feature service and we do not plan to.

So I wonder if there is an alternative way to make it work?

Update the data model:

for example, we get data from this url:

http://server/ws?bbox=30,30,40,40&xx=xx

Then we will get:

[{id:1,x:33,y:44,name:'xx'},{id:2.....}]


texts: Does GDALWARP have support for OpenCL speedup to orthorectification?

texts:
So I am attempting to orthorectify some images very large images. The process can take several hours so currently I am relying on using 8 threads to perform the speed up. However, the GPU does not appear to be in use at all according to nvidia-smi. I am very confused why this is the case:

Here is the script I am using:

#! /bin/sh

CUR_DIR="$(cd $(dirname ${BASH_SOURCE[0]}); pwd)"
export LD_LIBRARY_PATH="$CUR_DIR/../lib"
$CUR_DIR/gdalwarp -overwrite -rpc -to RPC_DEM="$1" -wo USE_OPENCL=TRUE   "$2" "$3"

I am attempting to orthorectify a large NITF which have I stored into an RPC model:

I have copied the different bands from GDAL info to see if I need to change the -wt or -ot to speed up the computation. I have recompiled GDAL with Nvidia's OPENCL and am running that recompiled version.

Why does GDALWARP not use the GPU in this case?

Band 1 Block=128x128 Type=UInt16, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=128x128 Type=UInt16, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=128x128 Type=UInt16, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=128x128 Type=UInt16, ColorInterp=Alpha


texts:
Does QGIS (or different software) have a method to extrapolate raster values to larger area (outer boundary)?

I have a boundary(.shp), which is 40 m far from the edge of raster. I want to fill the space between raster edge and polygon boundary with values extrapolated from edges of raster. Some cubic spline method would be ideal.


texts: I have a project in QGIS which involves interpolating several different variables on several different maps, then using the raster calculator to multiply some of the resulting interpolations together.  I'm now trying to automate it so I can create this map on a daily basis.  What specific libraries must be installed so I can run interpolation automatically (via a python script, or any other method)?  Can anyone share an example of how this can be done?


texts: British National Grid and WGS84 conversion doesn't match in reverse in GeoPandas/PyProj

texts:
When I try to convert Google Map's lat/long into EPSG:27700 and then add a 50-meter buffer, I'm getting unexpected results. I noticed that conversions in pyproj between WGS84 and BNG are not reversible:

import pyproj
#Google Map coordinates of Borders General Hospital: 55.595754, -2.741771

#Projections:
bng=pyproj.Proj(init='epsg:27700')
wgs84 = pyproj.Proj(init='epsg:4326')

#Google to BNG
pyproj.transform(wgs84, bng, 55.595754, -2.741771)
#Retuns
(8233971.054460906, -6088455.394871053)

#BNG to Google
pyproj.transform(wgs84, bng, 8233971.054460906, -6088455.394871053)
#Returns
(54.35845259494041, -3.0681078545288214)

What am I doing wrong?


texts: How to speed up queries for raster databases?

texts:
I have a raster database in postgresql/postgis with these columns:

(ID, rast, data_of_data).

'rast' is the column that has raster files in WKT format. An example query to find the DN value of a point in WGS84 system (30.424, -1.66) and for 2002-01-09 is the following:

SELECT
     st_value(rast,(st_GeomFromText('POINT(30.424 -1.66)', 4326))) as val
FROM
     my_table
WHERE
     date_of_data='2002-01-09'

Is there a method (eg. spatial index) to speed up those kind of queries?


texts:
enter image description here

I am following step by step the online tutorial on QGIS website to learn how to build a virtual raster using the images provided with the tutorial. I am on a Mac with QGIS 3.8, yet when running the algorithm 'Build Virtual Raster' I receive the following error "The following layers were not correctly generated" (see picture) and no layer is created...
Any clue?


texts:
I have a polygon layer that looks like this i.e. they overlap.

When I run an intersection between the layers, I get more polylines i.e. for the circle on the top, the three lines get split because of the circle on the bottom. So, I get the following output for the top 3 lines

So, the lines are split into 6 parts. I would like the intersection to be not be affected by overlapping polygons. Only what is contained in the polygon must be obtained as the output without the feature being split because of nearby polygons.

Is this even possible by intersection or is there another way to achieve this?

I don't really want to dissolve the polygons because I depend on FID for further processing.


texts: My external hard drive is corrupt - according to my IT Expert possibly because of the Fat32 format. Is FAT32 really not suitable for GIS data? ie large Tifs but also thousands of 1kb txt files (e.g World files). Are there any situations when FAT32 is prefered to NTFS?


texts: Cloud Removal through Sen2cor for sentinel 2 images?

texts:
I have this imagery for an area of germany. According to docs of Sen2Cor, it states it can only remove haze and light clouds, i wanted to ask would it be possible to remove something as dense as shown in the imagery with the processor?

How can i perform cloud correction on such imageries without losing data below?


texts: How to make Choropleth map using Openlayers 3?

texts:
I want to make Choropleth map in Openlayers 3 which is already implemented in Leaflet but not in Openlayers 3. how to integrate 'colorbrewer' in Openlayers 3?
Give the working example of it.


texts: Running processing workflow in PyQGIS for every selected polygon separately in for loop

texts:
I would like to run the processing chain written in PyQGIS (as stated below) for every polygon that is extracted by the extractbyexpression separately. But, the code below does not do this for every polygon separately, but rather n times for all polygons together which results in wrong output (with n the number of extracted polygons).
Any thoughts on how to solve this? So for each extracted polygon, the workflow is run separately? I am new to PyQGIS.
openruimte = iface.addVectorLayer(....)

expression_divide = ['Area >= 10 and Area < 20']

selectbyexpr = processing.run("native:extractbyexpression",
                              {'EXPRESSION' : expression_divide[0],
                               'INPUT': openruimte,
                               'OUTPUT': 'memory:'})['OUTPUT']

for feature in selectbyexpr.getFeatures():
    x = feature["Area"] / 5
    selectbyexpr.select(feature.id())
    random_points = processing.run(
        "qgis:randompointsinsidepolygons",
        {'EXPRESSION': '20',
        'INPUT': QgsProcessingFeatureSourceDefinition(selectbyexpr.id(), selectedFeaturesOnly=True),
        'MIN_DISTANCE': 0.01,
        'OUTPUT': 'memory:',
        'STRATEGY': 1})['OUTPUT']

    QgsProject.instance().addMapLayer(random_points)

    cluster= processing.run(
        "native:kmeansclustering",
        {'CLUSTERS': x, 'FIELD_NAME': 'CLUSTER_ID', 'INPUT': random_points,
        'OUTPUT': 'memory:'})['OUTPUT']

    aggregate= processing.run("qgis:aggregate",
                              {'AGGREGATES' : [{'aggregate': 'array_agg', 'delimiter': ',', 'input': '"CLUSTER_ID"', 'length': 0, 'name': 'CLUSTER_ID', 'precision': 0, 'type': 2}],
                               'GROUP_BY' : '\"CLUSTER_ID\"',
                               'INPUT' : cluster,
                               'OUTPUT' : 'memory:'})['OUTPUT']

    centroid_aggr= processing.run("native:centroids",
                            {'ALL_PARTS': False,
                            'INPUT': aggregate,
                            'OUTPUT': 'memory:'})['OUTPUT']

As can be seen from the picture, only 1 output of the QgsProject.instance().addMapLayer(random_points) is shown, but it includes the random points for all selected polygons and not only for 1 selected polygon.


texts: Combining two scenes with different extents without overlapping in QGIS

texts:
My AOI is split between two Sentinel-2 scenes. These rasters do not overlap. So function 'merge' in the Raster menu of QGIS cannot help - it says 'Layers do not have common extent'.

How to combine it to the one, bigger layer?

I tried in R as well, but no functions helped.


texts:
I am trying to generate opendrive (xodr) files from the OpenStreetMap data.
OpenStreetMap provides nodes and ways. In opendrive I need to compute the length of roads, heading angle, inertial x and y coordinates and the geometry records for the roads. Right now I am using osmnx module to visualize the road network.
Can I get some help regarding the generation of geometry records?


texts: Least cost path between various points QGIS

texts:
Is there an easy way to find the least cost path (LCP) between various points in an open source GIS?

This question has actually already been asked here: How to find a least cost path through several points?
https://gis.stackexchange.com/questions/172894/least-cost-path-through-3-points-start-destination-point-between-them
However this post aims to figure out if there´s really no direct way to get favoured results as there were no satisfying answers yet.

There are modules and plug-ins (e.g.: saga:leastcostpaths; grass:r.drain) for least cost path analysis (LCPA) that work well for the LCP between two destination points, but all of them I know and tried so far are not able to compute the LCP between various point.

The reason is that these LCPA-modules compute the path looking for the neighbor pixel with the lowest value of an accumulated cost raster. This works fine as long as there is only one direction (from one point to one other point) to go, however applying these LCPA-modules inputing a cost raster based on various points leads to unsatisfying results as the path got stuck immediatly or after a small amount of pixels.

A cost raster created by various points is like a hilly landscape with ups and down, though LCPA-modules are only able to find the way down they are not able to find the LCP between points, as these point are either located in valleys (so there´s a hills in between) or (at inverted raster) on top of hills (so there are valleys in between where LCP got stuck).

So, I´m looking for a solution how to connect various points by the least cost path based on a cost raster. Does anybody of you has experiences with this topic? Is there an easy way (module/plug-in that already is able to do that job)?

Using the LCPA-modules one cost layer (+LCP) for the distance between each combination of points would be needed. (resources and time consuming)

The created LCP shall build a network suitable for further network analysis.

Edit: In order to specify my question it should be added that LCPA with SAGA GIS (and GRASS GIS probably too) works fine to connect one point with several others but I´m looking for a way to connect several point within themselves.
The problem is that GIS tools for LCPA (at least th ones of SAGA and GRASS) are only able to way "downwards" at a cost raster from "peak" where there is the starting point the cost layer was created from, however creating a cost layer with various starting points results in a "hilly" layer where LCPA tools of SAGA and GRASS are not able to find favored paths as they are not able to "climb" but rather got stuck in local minimum pixel values.

Is there any GIS tool able to compute paths between point or polygons based on a cost raster where the points or polygons are located at several local pixel value peaks?


texts: Making intra-annual image composites for a series of years in Google Earth Engine?

texts:
Below is my code to build a time series collection of all ETM+ and OLI Landsat images intersecting my region of interest. Now I want to build an ImageCollection that represents the mean or median of all images within a year for all years. Can you provide an example of how to build intra-annual composites in Earth Engine.

var coefficients = {
  itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172]).multiply(10000),
  slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071]),
};

// Define function to get and rename bands of interest from OLI.
function renameOLI(img) {
  return img.select(
        ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
        ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']
    );
}

// Define function to get and rename bands of interest from ETM+.
function renameETM(img) {
  return img.select(
        ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
        ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']
  );
}

// Define function to apply harmonization transformation.
function etm2oli(img) {
  return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
    .multiply(coefficients.slopes)
    .add(coefficients.itcps)
    .round()
    .toShort()
    .addBands(img.select('pixel_qa')
    .copyProperties(img, ['system:time_start'])
  );
}

// Define function to mask out clouds and cloud shadows.
function fmask(img) {
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;
  var qa = img.select('pixel_qa');
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return img.updateMask(mask);
}

// Define function to calculate NDVI.
function calcNdvi(img) {
  return img.normalizedDifference(['NIR', 'Red']).rename('NDVI');

}

// Define function to prepare OLI images.
function prepOLI(img) {
  var orig = img;
  img = renameOLI(img);
  img = fmask(img);
  img = calcNdvi(img);
  return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}

// Define function to prepare ETM+ images.
function prepETM(img) {
  var orig = img;
  img = renameETM(img);
  img = fmask(img);
  img = etm2oli(img);
  img = calcNdvi(img);
  return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}

// Define a point on the study area
  var aoi = ee.Geometry.Point([34.0049766165165,-6.137332241818918]);

// Define AOI on the map.
Map.centerObject(aoi, 10);
Map.addLayer(aoi, {color: 'f8766d'}, 'AOI');
Map.setOptions('HYBRID');

// Get Landsat surface reflectance collections for OLI, ETM+ and TM sensors.
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol= ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');

// Define a collection filter.
var colFilter = ee.Filter.and(
  ee.Filter.bounds(aoi),
  ee.Filter.calendarRange(1, 365, 'day_of_year'),
  ee.Filter.lt('CLOUD_COVER', 50),
  ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
  ee.Filter.or(
    ee.Filter.eq('IMAGE_QUALITY', 9),
    ee.Filter.eq('IMAGE_QUALITY_OLI', 9)
  )
);

// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOLI);
etmCol= etmCol.filter(colFilter).map(prepETM);

// Merge the collections.
var col = oliCol
  .merge(etmCol);

Export.image.toAsset({image:col,
description:"collection",
assetId:"rungwa",
region:aoi.bounds(),
scale:10,
maxPixels:1e13});


texts:
mbpipe is a simple utility to batch process images in an MBTiles database (sqlite) by piping them to an image processing program.

But how does one install and use it?

It says here that mbpipe comes with the mbtiles node.js module, that is installed with the command npm install -g mbtiles.

However after doing that there is in fact something in /usr/local/lib/node_modules/mbtiles but the "mbpipe" command or script cannot be found anywhere on my machine.

What should one do...?

Edit: on the github page for node-mbtiles it apparently says that the mbpipe utility was removed from the mbtiles module, but no one seems to be following up on this issue.

Is there another utility to do the same? (It should be pretty straightforward: extract images from an sqlite db, send them to a program for processing, re-insert the processed image).


texts:
I am new to JSON, GitHub, and Leaflet.

I am trying to create a simple application which will render my points which are in GeoJSON up on GitHub.

From everything I have read this should work.  Here is what I have so far but I do not see my points displayed.

What am I missing?

<link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.2/leaflet.css" />


texts: How to get all the streets and their postal codes from a given state (New Hampshire) that has been loaded inside PostgreSQL

texts:
I downloaded .osm for state New Hampshire from here. I loaded it into PostgreSQL on Ubuntu 16.04 VB using osm2pgsql. When I use \dt on my db it shows this:

However SELECT COUNT(*) FROM county; returns 0. SELECT COUNT(*) FROM planet_osm_point; returns 55083.

I'd like to get a list of all towns and all streets in those towns. However, there are no counties, nevermind towns!

Also \d+ planet_osm_point shows this;

Then I tried SELECT amenity FROM planet_osm_point LIMIT 10 and it just returns empty space with 10 rows:

EDIT:

Managed to get towns, cities, counties, types of highway. But can't get streets for a given city.


texts: ArcGIS 10.0 Python script won't recognize null

texts:
Has anybody ran into in issue where Python field calculations and/or scripts won't properly query Null (None) values?

The calculation I'm working on is a little more complex, but even the simple calculation below doesn't produce the correct results:

codeblock:

def test(field):
    if field is None:
        return "value is null"
    else:
        return field

expression:

test(!field!)

The result from above is:

input row 1 = <null>, output = <null>
input row 2 = "", output = "value is null"
input row 3 = "test value", output = "test value"

I can set value(s) in a field to null using:

Expression

None

I broke down the calculation I was working on in the Python Command Line to try to isolate how values were being interpreted.

>>> arg1, arg2, arg3 = "Testing", None, "this string"
>>> mylist = [arg1, arg2, arg3]
>>> remnull = [str(arg) for arg in mylist if arg]
>>> sep = " "
>>> finalval = sep.join(remnull)
>>> print finalval
Testing this string
>>> anullval = None
>>> if anullval is None:
...     print "yes"
...
yes
>>> if not anullval:
...     print "yes"
...
yes

It would appear that the Python Command Line is interpreting None correctly. Python scripts however won't recognize None values (using calculatefield_management or cursors), either as standalone script or script tool.

I've been doing a lot of coding lately.  I'm wondering if a rogue line of code hasn't corrupted something.

I can do a uninstall/re-install or repair but I'd like to spare myself a repeat if somebody knows the cause.

I'm on ArcGIS 10.0 SP5, Basic.  All data inputs are File Geodatabase.


texts: Convert .grib2 in triangular grid to regular grid

texts:
I have ICON datasets (.grib2.bz2 file format) from Deutscher Wetterdienst (DWD) in triangular grid that I would like to transform to a regular grid. Currently, I run the following statement in my Python script:
subprocess.call(['./transform.sh', blobName])

In this bash script I have the following:
#!/bin/bash
bzip2 -d ${1}

export in_file=${1//.bz2/}
export out_file=${in_file//icosahedral_single-level/world_0125x0125}

time cdo -f grb2 remap,settings.txt,weights.nc ${in_file} ${out_file}

DWD provides details on how to achieve this here. I would like to improve my script in two ways:

Run everything in python itself
Do not use physical files

This means I would like to download the .grib2.bz2 file in memory and produce the output in memory so I can store it in Azure Blob Storage.
Is this feasible? Maybe with rasterio?
For 1) I could use python cdo but this does not solve 2) as far as I know. I was also not able to run cdo on Python 3.X version.


texts:
i have GPS data from 15 trucks for one month. The data include (location x-y, speed, altitude and vehicle_id). I have used tracking analyst to visualize the data but i need to group them by day to show each truck's trip in detail. Because not all the cars travel every day and same hours.. So i get no result if i put them altogether.

How can i check that my GPS data has no bias?

i observed that there are cases when a truck makes a trip from a to b it turns off the gps at b and the other day it starts again from a. So when i get ArcGIS to calculate distances it appears a wrong (and often too big) distance.

How can i fix that? Should i split my dataset by day?


texts: GeoServer v2.14 Responses to DescribeFeatureType Requests

texts:
A browser's (Chrome) DescribeFeatureType request to GeoServer version 2.12 returns an inline XML response in the browser's window.

Using that same browser to send a DescribeFeatureType request to a GeoServer version 2.14, the response is a document that is downloaded to the Windows' Downloads directory with the filename geoserver-describefeaturetype.application. This document contains the same XML content as the inline XML response when requesting from GeoServer version 2.12.

Coincidently, QGIS is not able to resolve the DescribeFeatureType request to the GeoServer version 2.14 but does resolve the DescribeFeatureType request from GeoServer 2.12.

We don't see any settings in GeoServer to modify how it responds to DescribeFeatureType requests.

Interestingly, a browser's GetCapabilities request to both versions of GeoServer yields an inline XML response from GeoServer.

The DescribeFeatureType request is in the form https://myserver.com/geoserver/mynamespace/wfs?request=DescribeFeatureType

and the GetCapabilities request is in the form

https://myserver.com/geoserver/mynamespace/wfs?request=GetCapabilities

Most importantly, we want to get the WFS connection from QGIS to GeoServer version 2.14 working.  But we also want the inline response to a DescribeFeatureType request.

Any suggestions?


texts:
I have raster function for computing the EVI index from a Sentinel-2 raster located in a mosaic dataset. The code is similar to https://github.com/Esri/raster-functions/blob/master/functions/NDVI.py
I added the block for the blue band:

            {
                'name': 'blue',
                'dataType': 'numeric',
                'value': 1,
                'required': True,
                'displayName': "Blue Band Index",
                'description': "Index for Blue band"
            },

The EVI formula is:

outBlock = 2.5 * (nir - red) / ((nir + 6*red - 7.5*blue) + 1)

Obviously this does not raise a DivisionError because it's impossible due the the '+1' in the formula.

When i run the function on the ArcGIS dataset, i get the error:

 outBlock = 2.5 * (nir - red) / ((nir + 6*red - 7.5*blue) + 1)

  FloatingPointError: divide by zero encountered in divide

How can I solve this?


texts: QGIS : copy/paste with "Avoid Intersection" option

texts:
I need to copy and paste polygons from one layer to another but with adjustments in the pasted geometries: I need them to be cut in the same way they would have been cut by a single edit with the "Avoid Intersection" option in the project snapping options.

In Arcmap I can nearly do this with the "construct feature" tool (topology toolbar) user the target layer existing objects. But in my QGIS project I need the use a third layer to cut the pasted polygons.

Tanks


texts: Openlayers WFS events: click or hover functionality not registering

texts:
I've successfully gotten this to work before so I am confused as to why it's not working this time. Everything seems to be in place: no mapfile errors, no js errors in the console.

Here's a gist of the js: https://gist.github.com/anonymous/7664194

Lines 177-217 being the ones of concern.

Note: the points are displaying correctly, but I just can't get any events to register on the points, whether it be hover or click.

var wfs = new OpenLayers.Layer.Vector("SW Airports", {
            wrapDateLine: true,
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.WFS({
                    srsName: "EPSG:3857",
                    version: "1.0.0",
                    url: "http://another_server/cgi-bin/mapserv.exe?map=/tmp/navhome/navcore/interactive/interactive_viewer.map&SERVICE=WFS&VERSION=1.1.0&srsName=EPSG:3857",
                    featureType: "Airports"
            })
    });
map.addLayer(wfs);

var selectControl = new OpenLayers.Control.SelectFeature(
            [wfs],
            {
                    clickout: true, toggle: true,
                    multiple: true, hover: true,
                    toggleKey: "ctrlKey",
                    multipleKey: "shiftKey"
            }
    );

    map.addControl(selectControl);
    selectControl.activate();

    wfs.events.on({
    featureselected: function(event) {
            var feature = event.feature;
            feature.popup = new OpenLayers.Popup.FramedCloud("box",
                    feature.geometry.getBounds().getCenterLonLat(),
                    null,
                    '
Airport: '+feature.attributes.name +'
 IATA: '+feature.attributes.iata+'
'+'
'+ 'Center on Location'+ '
'+'Center on Location and Zoom In',
                    null,
                    true
            );
            while( map.popups.length ) {
                    map.removePopup( map.popups[0] );
            }
            map.addPopup(feature.popup);
            }
    });

Any insight on the matter is appreciated.


texts:
How to count the number of lines that are inside a polygon?

The methods I have tried don't solve the problem because they either count the FID of the line instead of the segment or cut the line at its intersection.


texts: Batch create multiple buffers with variable size around point

texts:
For a statistical analysis, I need to calculate some landscape properties in several (rectangular) sample areas, which vary in size (extent ranging from 10 by 10 to 500 by 500). These sample areas are located (scattered) randomly across a larger study area. For this analysis, I want to create several vector layers, in which each layer contains all sample areas with equal extent (every sample area being a rectangular polygon). All starts with a point-vector layer, each point being the center of a (potential) sample area.
I think the easiest way to create those rectangular polygons is by creating a (rectangular) buffer around every point. Pretty straightforward and lean.
But I would like to do this "in batch" since I want to create 50 layers. And this is where I'm stuck in QGIS. There's a very handy "batch processing" function for every tool, but I can't figure out how to run this for this (rather simple?) process of creating buffers around points.
Any tips on this?


texts: Why does FeatureClassToFeatureClass work in ArcMap but not in Python script?

texts:
When trying to copy a feature class from an SDE connection file to a scratch GDB, the FeatureClassToFeatureClass_conversion works in the ArcMap Python window, but I can't get it to work in a Python script. In ArcMap, the code looks like this:

pathScratchGDB = arcpy.env.scratchGDB + os.path.sep
arcpy.FeatureClassToFeatureClass_conversion("C:/SmartFarm_Tools/smartfarm_connection.sde/smartfarm.master.zone_catalog", pathScratchGDB, "fcZones", 'fldid_id = ' + str(504) + " and type = 'BW2014'")

which results in:

<Result 'C:\\Users\\username\\Documents\\ArcGIS\\scratch.gdb\\fcZones'>

In my Python script I have the exact same code:

pathScratchGDB = arcpy.env.scratchGDB + os.path.sep
arcpy.FeatureClassToFeatureClass_conversion("C:/SmartFarm_Tools/smartfarm_connection.sde/smartfarm.master.zone_catalog", pathScratchGDB, "fcZones", 'fldid_id = ' + str(504) + " and type = 'BW2014'")

which results in:

Exception in Tkinter callback
Traceback (most recent call last):
  File "C:\Python27\ArcGISx6410.3\lib\lib-tk\Tkinter.py", line 1486, in __call__
    return self.func(*args)
  File "C:\SmartFarm_Tools\python_scripts\SmartFarm_Engine1.py", line 1228, in onLoad
    arcpy.FeatureClassToFeatureClass_conversion("C:/SmartFarm_Tools/smartfarm_connection.sde/smartfarm.master.zone_catalog", pathScratchGDB, "fcZones", 'fldid_id = ' + str(504) + " and type = 'BW2014'")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\conversion.py", line 1789, in FeatureClassToFeatureClass
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset C:/SmartFarm_Tools/smartfarm_connection.sde/smartfarm.master.zone_catalog does not exist or is not supported
Failed to execute (FeatureClassToFeatureClass).

In the same Python script I am able to successfully run this line of code:

arcpy.FeatureClassToFeatureClass_conversion("C:/SmartFarm_Tools/smartfarm_connection.sde/smartfarm.master.boundary_2015", pathScratchGDB, "fcBound", 'master_id = ' + str(inMasterID))

Has anyone seen this before?

EDIT

The data resides on a PostgreSQL/PostGIS database.


texts:
I have a hierachical List of attributes e.g. column (1) S (2) SD (3) SD0
or (1) A (2) AN (3) An1 etc. So there is one S with several SD s and several SD0 in the attribute table. When I display column (3) I wonder if there is an automatic way to set the colours in a stacked way i.e. all S different grey colours all A different blue etc. All I could find was the rule based way but then it is still a lot of work because there are so many types in column 3. Is there a way to use the other columns, in my case for example column 1 to classify?

Thank you for help!!


texts: Change marker within leaflet map when loading data from geojson file

texts:
This section of code adds markers to my map (it works), how can I change the marker to a different icon when dealing with geojson?

var cycleracks = $.getJSON('cycle-racks.geojson', function (geojson) {
  L.geoJson(geojson, {
    onEachFeature: function (feature, layer) {
      layer.bindPopup(feature.properties.type);
    }
  }).addTo(map);
});