texts:
I work at urban planning and I am trying to create an hexbin map from landscape use of a small city.

I have different polygons (vector layers) like buildings, roads, gardens, services, industry, and so on... that doesn't overlap.

I also have created a HEXAGONAL GRID (polygons), and each hexagon has 10.000 m2 (1 ha).

My idea is to create something like this:
https://pbs.twimg.com/media/CEiZ9eEVIAApcSO.png

A cool and very interesting alternative to point-based maps, in my opinion...

My first (1) option was to go to:
menu -> Vector -> Geometry tools -> polygon centroid and create a centroid point layer, and then "count points in polygons" but I came to the conclusion that much of the information I have, like what is the use that is predominant inside each hexagon, is lost...

So I am now searching for a method that counts area (% of the 10.000 m2) of the different uses... For example for green areas, what is its predominance inside each hexagons?

How can I "pass" that attribute to the hexagonal shape?

Should I clip building-green areas-industry layers, with the grid, and the perform a spatial join (join by location)?

Or is there another way to perform this?

I am using QGIS 2.14 ESSEN.

____________________EDITING____________________________

To specify my question:

I want to know if it is possible to calculate all red-pink-blue-green area inside each hexagon, knowing that every hexagon has 10.000 m2.

My objective is to know for hexagon ID = n:

Building area  = x %
Green area  = x %
Industry area  = x %
Roads area  = x %

And for all other 450 Hexagons.

Do I have to clip with the grid and then perform spatial join?


texts: Transferring layer files from ArcGIS Pro to QGIS

texts:
I have access to ArcGIS Pro at my workplace and will be transferring a project over to someone who will be using QGIS.
I have saved my layers as a layer file but the problem is that the other person is unable to import .lyrx files to QGIS.
Does anyone know how I can transfer the project to them and include all the symbology settings I have created?


texts: Is it possible to symbolize a layer with time from now?

texts:
Using Arc Desktop 10.3.1 and AGOnline, I have a polygon feature class showing different service zones and the dates of last service.  I would like to symbolize this using a time from now,  ie if it has been more than 1 month since last service date it will be yellow, if it has been more than 2 months it will be orange, etc...  Is there a way to put an expression into the symbology like

if (now) - [last date] > 30  & < 60 then yellow

if (now) - [last date] < 90 & > 60 then orange

etc...

This map will ultimately be used in an ArcGIS Online web map, so it has to be dynamic and constanly updating.


texts: Does leaflet support user defined Layer identifiers

texts:
It would be nice to be able to check for a layer presence or remove it without having to keep an explicit reference to it. To that extent I would like to use ad id for the layer that I can compute myself from my data. Is this feasible with leaflet? If no, do you have any clean workaround?

For example with a LayerGroup, the idea is to have a setId(<String> id) function that would allow leveraging the native removeLayer( <String> id ) and getLayer( <String> id ) methods to work on the LayerGroup content by user defined identifiers. Example below:

marker = L.circleMarker(...);
marker.setId('my_marker_id');
lg = L.layerGroup(...);
lg.addLayer(marker);
...
lg.removeLayer('my_marker_id');


texts: Overlaying spatial polygon with grid and checking in which grid element specific coordinates are located using R

texts:
How can one use R to

split a shapefile in 200 meter squares/sub-polygons,
plot this grid (incl. ID numbers for each square) over the original map below, and
evaluate in which square specific geographic coordinates are located.

I am a beginner in GIS and this is perhaps a basic question, but I haven't found a tutorial on how to do this in R.

What I have done so far is loading a shapefile of NYC and plotting some exemplary geographic coordinates.

I am looking for an example (R code) how to this with the data below.

# Load packages
library(maptools)

# Download shapefile for NYC
# OLD URL (no longer working)
# shpurl <- "http://www.nyc.gov/html/dcp/download/bytes/nybb_13a.zip"
shpurl <- "https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nybb_13a.zip"

tmp    <- tempfile(fileext=".zip")
download.file(shpurl, destfile=tmp)
files <- unzip(tmp, exdir=getwd())

# Load & plot shapefile
shp <- readShapePoly(files[grep(".shp$", files)])
plot(shp)

# Define coordinates
points_of_interest <- data.frame(y=c(919500, 959500, 1019500, 1049500, 1029500, 989500),
                 x =c(130600, 150600, 180600, 198000, 248000, 218000),
                 id  =c("A"), stringsAsFactors=F)

# Plot coordinates
points(points_of_interest$y, points_of_interest$x, pch=19, col="red")


texts: Does QGIS require a specific name for the geometry field in a PostGIS table

texts:
I seemed to have noticed that QGIS cannot detect the spatial type or SRID of a PostGIS table if the geometry column is set to anything other than "geom". You can manually tell QGIS what these fields should be but then upon opening a saved project the tables in particular seem slow to load.
Is it recommended that "geom" is used for the geometry column to speed up the table reads or could this be something un-related.


texts: Update field of attribute table of a shape in selected rows with arcpy

texts:
I read the post Accessing attribute table within shapefile and replace values?, but I don't know how to include a "Where clause".

I want to update the attribute table of a shapefile. In SQL the Code looks like this:

alter table point_shape
    add buffer number;

update point_shape
    set buffer = a/2
    where a >0;

alter table line_shape
    add buffer number;

update line_shape
    set buffer = a/2
    where a >0;

alter table point_shape
    add buffer number;

update point_shape
    set buffer = sqrt(b)/acos(-1)
    where b is not null and b >0;

alter table line_shape
    add buffer number;

update line_shape
    set buffer = b/shape_len
    where b is not null and a = 0 and b>0;

I translated the first update statement:

cur = arcpy.UpdateCursor(point_shape)

for row in cur:
    row.setValue('buffer', row.getValue('a') / 2)
    cur.updateRow(row)

But where can I include the where statement?


texts: Can SqlWorkspaceFactory connection files be created with ArcGIS Desktop?

texts:
With ArcObjects, I can use IWorkspaceFactory.Create I can create a connection file for a SqlWorkspaceFactory.

Is there an out-of-the-box method to create the same file in ArcGIS Desktop?

Update: I tried the Create SDE Connection GP tool approach.  I am unable to Open with the OpenFromFile method.  The files are the same length but different when compared.

Update 2:  I experience the same result in ArcGIS 10.2.


texts: gdalinfo --formats return different values in python3 than in OSGeo4W

texts:
I am working whit geopdf to obtain the layers. When I use driver = ogr.GetDriverByName('PDF') it return None, and obviously I can't open it.

I have checked with the command gdalinfo --formats and it is the result in OSGeo4W

and this another when I use python 3.6 whit anaconda prompt

someone know why is completely different?
I have downloaded gdal from anaconda.


texts: How to add a legend for mapfile (Raster Data)

texts:
What I would like to do is shown in the figure below. How can I add this color ramp legend in a mapfile which will also be linked to a raster. I don't know how to get it working while drawing a figure based on that mapfile.

Thanks !


texts:
I have an ArcObjects extension that creates a selection set and then zooms to the selected features.  This part works as expected.  The problem that we are having, however, is that the selection set does not display in the Attributes Window as it would if you were to use the normal selection tools in ArcGIS.

I looked at the IAttributeWindow interface but it seems that this is more concerned with showing that window rather than populating it.

Does anyone know which interface I need to look at for this, or if anyone can point me to an example?

We're running ArcGIS 10.2.1.


texts: I am looking for a way to loop through GRASS gis v.net.distance in QGIS. I want to iterate over all the 3 layers which has multiple features at the same time, and print the shortest distance for each of the layers in a CSV format. Can anybody please refer me to some resources or links with examples of how to do it in PYTHON? is it possible to do it with PYQGIS at all?


texts: Error when using Fill Sinks (wang & liu) in Saga in QGIS

texts:
I am using fill sinks (wang & liu) from Saga in QGIS.
When I run the tool I am having the following error:

  Loading resulting layers The following layers were not correctly
  generated.C:/AGMEE/Workspace
  Saga/DEM_FlowDirection.sdatC:/AGMEE/Workspace
  Saga/DEM_WatershedBasins.sdatYou can check the 'Log Messages
  Panel' in QGIS main window to find more information about the
  execution of the algorithm.

However, when I run the tool the Filled Dem parameter is executed but the Flow Direction and Watershed Basins are not.

I am using the last version of OSGeo4W Network.

Do you know something about it?


texts: ST_ClusterDBSCAN gives all features a cluster

texts:
I am using ST_ClusterDBSCAN to help cluster some building polygons. I am using the following SQL:

CREATE TABLE clusters AS SELECT ogc_fid, cid, wkb_geometry AS geom FROM (
    SELECT ogc_fid, ST_ClusterDBSCAN(wkb_geometry, eps := 5, minpoints := 1) over () AS cid, wkb_geometry
    FROM building) sq;

When I check the results every feature has a cid even though there are no other features with the same cid, i.e a building is not within 5 metres of another.

Reading the docs, ST_ClusterDBSCAN should give a feature not within a cluster the value of null.

So I am not sure why this is happening.


texts: Making 'persistent' ArcMap tool without creating Python add-in?

texts:
I want to make an addon that, while editing in ArcMap, will track where I pan and move the viewport. I have the pieces together, but I don't like debugging with the Python Add-In wizard, because I have to recompile and restart ArcMap each time I modify the code.

Is there another way to create this type of 'persistent' script without using an add-in? If I make a toolbox tool and run it, will it run in the background, and how would I stop it? All the toolbox tools I've ever created have a certain task they perform and, once completed, the tool stops running. Same with ModelBuilder.


texts:
Is there an analogue of ArcGIS Dissolve function in R? I need to dissolve the overlapping polygons and create a single part object.

I was trying to use gUnionCascaded, but it creates a multipart object. Is there a function in R to do Multipart to singlepart?


texts:
I have a point shapefile with "Lat" and "Lon" attributes with more than 1000 features.

I want to add street address data in the shapefile by Reverse Geocoding but the Geocoding plugin has no such option. Can someone suggest any other tool to create address data?

I have CSV-files as well.


texts: Problems configuring MapProxy to serve WMS layer

texts:
I'm having problems making MapProxy render a Mapserver WMS layer. My
configuration file looks like this:

services:
 wms:
   srs: ['EPSG:3785', 'EPSG:4326', 'EPSG:23028']
   image_formats: ['image/jpg', 'image/png']

layers:
 name: mylayer
 title: My layer
 sources: [my_wms_source]

sources:
 my_wms_source:
   type:wms
   req:
     url: http://mydomain/cgi-bin/mapserv?
     layers: mylayer_in_mapserver
     map: /path/to/my/mapfile.map

Direct queries to mapserver work fine. So, a query like this works:

http://mydomain/cgi-bin/mapserv?styles=&map=/path/to/my/mapfile.map&format=image/png&width=800&height=600&bbox=-18.4267,27.3264,-13.0047,29.5457&layers=mylayer_in_mapserver&service=WMS&request=GetMap&srs=EPSG:4326&version=1.1.1

The key, I think, is in the bounding box. If I deploy MapProxy with
debug activated and try to see my WMS source through the demo app,
queries like these are raised:

http://mydomain/cgi-bin/mapserv?styles=&map=map=/path/to/my/mapfile.map&format=image%2Fpng&width=800&height=600&bbox=-506.77734375,-266.8359375,370.01953125,190.1953125&layers=mylayer_in_mapserver&service=WMS&request=GetMap&srs=EPSG:4326&version=1.1.1

And that query doesn't work, because the bounding box is in a
different coordinate system (I guess epsg:3785).

How could I control that? The bbox is used for WMS queries. I've tried
the 'srs' parameter in the wms option of the 'services' section, but it
doesn't work. Any clues?

Thanks in advance, and best regards


texts: Saving view to different schema using QGIS DB manager?

texts:
When I ask to save a query to a view QGIS pops up a dialogue box asking the for a name and then saves the view in the public schema with that name.  I have tried "layers.name" but this just puts "layers." in front of the name and saves it in public.

The only way I have found to get the views created in the layers schema is to use postgres directly.

Am I missing something?


texts: sorting polygons extracted using Fishnet

texts:
I generated polygons using fishnet in ArcMap, since the polygons didn't fall exactly on top of my raster image, I had to edit and move some of them. Now the FID of my polygons are all messed up, and are not in the order that I want them to be (I want them to start from top right and end at the bottom left).

I tried calculating the centroid of the polygons and sort them based on that, but it doesn't work since not all the polygons are in line. I attach a picture.


texts: Computing convex hull of points using shapely

texts:
I have a points shapefile and I want to compute the convex hull of the points using shapely, and output the result as a shapefile using Fiona.

Does anyone know how I can do that?

This is what I have tried:

from shapely.geometry import MultiPoint
from fiona.crs import from_epsg
shape = fiona.open('.cache/town_points/ny_points.shp')
with fiona.open('.cache/mew.shp', 'w', 'ESRI Shapefile', shape.schema.copy(), shape.crs) as output:
    for feat in shape:
        ch_area = MultiPoint(shape).convex_hull
        output.write(feat)


texts:
While trying to round out my answer to Geometry Generator using polygon boundary expression only returning Exterior Ring QGIS I am trying to determine if the geometry passed into the function is a Polygon or a MultiPolygon and then to extract it's rings to pass to the zigzag drawing code.
My first attempt was to use the solution to Remove inner rings from QgsGeometry polygon in QGIS Python gives me 'list' object has no attribute 'interpolatePoint' as an error with my code looking like:
geoms = []
if geom.isMultipart() is False: # if only simple polygon, calculate only for this
    polyg = geom.asPolygon() # transform to list of points
    for ring in polyg:
        geoms.append(ring)
else: # is multipart
    multi = geom.asMultiPolygon()
    for polyg in multi:
        for ring in polyg:
            geoms.append(ring)
# interpolate points on linestring
for gx in geoms:
    points2d = [(lambda g: (g.x(), g.y()))(gx.interpolatePoint(d)) for d in distances]
    vertices = gx.points()

So one of the asX methods returned an empty list instead of a Polygon, I think?
So next I tried to go from basics with:
geoms = []

# assume it is a polygon
if geom.wkbType() == QgsWkbTypes.Polygon:
    poly = geom.asPolygon()
    geoms.append(poly.exteriorRing())
else:
    geoms.append(geom.boundary())

Which gives the same error message when applied to polygon features.
So then I tried:
geoms = []

geoms.append(geom.exteriorRing())
for i in range(0,geom.getNumInteriorRings()):
    geoms.append(geom.getInteriorRingN(i))

But it says, Eval Error: 'QgsGeometry' object has no attribute 'exteriorRing' so I thought it needed changing to a polygon.
Finally, I tried:
geoms = []
poly = geom.asPolygon()
geoms.append(poly.exteriorRing())
for i in range(0,poly.getNumInteriorRings()):
    geoms.append(poly.getInteriorRingN(i))

But that just says 'list' object has no attribute 'exteriorRing'
My test features are:
Polygon ((513366.96000000002095476 102956.35000000000582077, 513355.33000000001629815 102954.64999999999417923, 513356.32000000000698492 102947.89999999999417923, 513367.95000000001164153 102949.60000000000582077, 513366.96000000002095476 102956.35000000000582077))
Polygon ((513392.47999999998137355 102973.47999999999592546, 513401.75 102916.32000000000698492, 513424.30999999999767169 102919.85000000000582077, 513415.40000000002328306 102977.08000000000174623, 513392.47999999998137355 102973.47999999999592546),(513408.89000000001396984 102967.58999999999650754, 513415.45000000001164153 102926.5, 513407.88000000000465661 102925.28999999999359716, 513401.32000000000698492 102966.38000000000465661, 513408.89000000001396984 102967.58999999999650754))

So (currently) there is no problem with a multipolygon or anything complicated.


texts: Matching color of pie charts with another set of data?

texts:
I am working with QGIS and as part of the project I would like to make the pie charts give me to sets of information.

I work % of people commuting by bike and I would like the volume of the pie to represent the %age of people commuting in 2015 and the color of the pie to be the growth between 2005 and 2015. So for instance, for Paris, I would like the pie volume to correspond to 1%, while its color would be blue because it grew significantly since 2005. I have multiple cities so I cannot do it by hands.

I looked extensively and could not find anything. Any idea?


texts:
I am using health mapper to display health data. It uses ArcGIS to draw district boundaries as a shapefile. I am interested to to draw sub district boundaries, e.g County or neighborhood boundary shapefile.

As ArcGIS is not free, any idea how to do it using a free software applied on Google Maps to define the county?


texts: Radiometric Resolution of Landsat Sensors

texts: Landsat data is in 8-bit (0-255) format. Would the difference between the sets of visible-NIR bands for Landsat 4,5, and 7 vary significantly?


texts: GeoDatabase topology and rules with a *.rul file

texts:
I have been working for some time on a solution to a problem much like this question.
I have several layers that I want to resolve together.  I have them all in topology (geodatabase) and have been trying to decide which rules and in what order would be best practice or even would accomplish.
For instance the same data in the overlapping sections needs to not overlap or have gaps.
That data should not move where a county polygon is near and should be moved to it.
As should the school districts.
All data layers are polygonal.
Ending with all features being moved and made exact (within tolerance) to the plss layer which is first repaired for gaps and slivers.  The topology tool allows you to load a *.rul file, which I assume is an ordered XML.  I haven't found one anywhere to examine it.
To my question does anyone hold a *.rul file that would begin to give me insight into the rules and the order (the tolerances I will most likely have to tweak for my data)?

If user823 would ever come back by I might convince him to post his complex rule set.


texts: Getting point attribute from polygon layer using QGIS?

texts:
I have two layers, one (A) has point features and the other (B) has polygons.

How do I produce a new table with some (the ID) or all the fields from the A table and some (the area code) or all from the B table where the B polygon contains the A point?

I am using QGIS.


texts: Implementing UTFGRID with MapServer and OpenLayers?

texts:
I'm having troubles using utfgrid module with openlayers. The implementation on Mapserver's side is quite simple, even if i have doubts about the imagetype :

MAP
    IMAGETYPE   UTFGRID
    SIZE        2000 2000
    CONFIG "MS_ERRORFILE" "/ms4w/tmp/ms_error.txt"
    #DEBUG 4
    EXTENT -176.197 -31.0869 181.47 81.8381
    NAME "allexmaple"

    PROJECTION
        "init=epsg:4326"
    END
    QUERYMAP
        STATUS ON
    END
    OUTPUTFORMAT
        NAME "utfgrid"
        DRIVER UTFGRID
        MIMETYPE "application/json"
        EXTENSION "json"
        FORMATOPTION "UTFRESOLUTION=4"
        FORMATOPTION "DUPLICATES=false"
    END
    WEB
        IMAGEPATH "tmp/"
        IMAGEURL "tmp/"
    END
    SYMBOL
        NAME "mysymbol"
        TYPE pixmap
        ANCHORPOINT 0.5 0.5
        IMAGE  "icons/mysymbol.png"
    END

    LAYER
        TOLERANCE 10
        CONNECTIONTYPE postgis
        NAME "example"
        CONNECTION "host=127.0.0.1 dbname=ais_data user=postgres password=admin port=5432"
        DATA "geom from mytable using srid=4326"
        STATUS ON
        TYPE POINT
        UTFITEM   "id"
        UTFDATA   "{\"id\":\"[id]\",\"name\":\"[name]\"}"
        TEMPLATE 'query.html'
        HEADER 'header.html'
        FOOTER 'footer.html'
        CLASSITEM "ct"

        CLASS
            EXPRESSION "mysymbol"
            STYLE
                ANGLE [angle]
                SYMBOL 'mysymbol'
            END # style
        END # class
    END
END

Openlayers, on the other hand, doesn't have an efficient example. I don't know what to provide in my URL. Here's what I tried:

var gridSource = new ol.source.TileUTFGrid({
       url: 'http://localhost:5000/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/essai.map&layers=example&tiled=true&mode=map&type=utfgrid&format=application/json'
      });
 var mySuperLayer = new ol.layer.Tile({
           source: gridSource
      });

But nothing appears and my URL stays in pending until timeout


texts: If-statement in Google Earth Engine is not working properly

texts:
In the following link, is a code where I'm trying to build a simple Google Earth Engine app where a button control's a selector's status. The button when clicked, it should check whether a predefined list of dates is empty or not and reflect that on the selector's condition. In order to achieve this goal, an if-statement is used but it's not working properly (line 52-55). It always return as if the condition is false even if it's true.
https://code.earthengine.google.com/e498bf6f86dd1a8413e0e19e42d561e5
Here is a snippet of interest from the code:
//Create selector
var av_dates = ui.Select({placeholder: "",style: SELECT_STYLE});
av_dates.setDisabled(1) //Initally disabled

//var dates =[ee.Date(Date.now()).format('YYYY-MM-dd').getInfo()]
var dates =[]
print('Dates length' , ee.Number(dates.length))
print('Condition status' , ee.Algorithms.IsEqual(ee.Number(dates.length), 0))

//Create button which controls selector's status
var CloudButton = ui.Button({label: 'Check availabel dates',
style:BUTTON_STYLE,
onClick:function() {
  ee.Algorithms.If(ee.Algorithms.IsEqual(ee.Number(dates.length), 0),
    av_dates.setDisabled(1) && av_dates.setPlaceholder('Non!'), // if no dates availabel then leave selector disabled
    av_dates.setDisabled(0) && av_dates.items().reset(dates) //otherwise show available dates and show the first one in the list
    && av_dates.setPlaceholder(dates[0]))
}
})

Why if statement is not working properly? It always act like the case is false!
Update:
According to Gennadii Donchyts insightful answer, I have modified the code accordingly from server-side to client-side so it does the job as expected as shown below:
//Create selector
var av_dates = ui.Select({placeholder: "",style: SELECT_STYLE});
av_dates.setDisabled(1) //Initially disabled

//var dates =[ee.Date(Date.now()).format('YYYY-MM-dd').getInfo()] //For testing purposes
var dates =[]
print('Dates length' , dates.length)
print('Condition status' , dates.length == 0)

//Create button which controls selector's status
var CloudButton = ui.Button({label: 'Check availabel dates',
style:BUTTON_STYLE,
onClick:function() {
if (dates.length == 0) {
  av_dates.setDisabled(1) && av_dates.setPlaceholder('Non!') // if no dates availabel then leave selector disabled
} else {
  av_dates.setDisabled(0) && av_dates.items().reset(dates) //otherwise show available dates and show the first one in the list
  && av_dates.setPlaceholder(dates[0])
}}
})
infoPanel.insert(1,CloudButton)
infoPanel.insert(2, av_dates)


texts: Programatically add polylines from existing points in QGIS

texts: I am using QGIS 2.18 and I need help finding a python script or tool which will allow me to programmatically draw a line from existing GPS (sometimes as many as 500 or more) points that is 250 feet long at a bearing listed in an azimuth attribute field for each point.


texts: UK National Grid references and GeoDjango

texts:
Newbie question. I'm trying to insert some UK National Grid references into a Django PointField defined as follows:

oscode = models.PointField(srid=27700, null=True, blank=True)

My question is this: can I store them in their original form ('TR3241', for example) or do I have to convert them to lat/lng first?

I don't know how to format them as WKT. This is what I get if I try inserting them as POINT(TR3241), which is obviously wrong:

INSERT INTO places (placeid, structidx, subidx, county, name, oscode) VALUES ('10', '1', '1', 'Kent', 'Dover', 'TR3241');
LINE 1: ...'1', 'Kent', 'D1', 'Eastry', 'Bewsbury', 'Dover', 'POINT(TR3...
                                                             ^
HINT:  "POINT(" <-- parse error at position 6 within geometry

So: is there a way that I can store points like 'TR3241' in GeoDjango, or do I have to convert them to standard lat/lng first?

I do have code to convert National Grid > lat/lng, but I  was hoping GeoDjango would be able to store them in the original format, and do any conversion for me.


texts: Calculating distance between each pixel of rasterstack image?

texts:
I am new on R and I have an image of class "rasterstack" of only 3 colors: Red, Green and blue as illustrated:

I use this function:
image = stack(img.red, img.green, img.blue)

I would like to find the distance between each pixel of one color to each others pixel of another color. By example the distance between each red pixel, to each blue pixel:

Distance between red pixel 1 to blue pixel 1 ,pixel 2 ,pixel 3...
Distance between red pixel 2 to blue pixel 1 ,pixel 2 ,pixel 3...
etc

And create a table with all these distances. Then I could do some statistic as to know the distribution of the distance, the average min distance between red and blue, etc

I have no idea how to do it?


texts: How to merge single polygons to multi part polygons by maximum distances?

texts:
I have a polygon layer with a land-use class and ID for every polygon. How can I merge polygons (into multi part polygon) which do not adjoin, but are not further than x meters apart from each other.

As an example you can imagine two forest polygons which are separated by a road. I would like the two forest polygons to be one multi-part polygon.

Of course one can merge two polygons by hand, but this would take forever!

Clarification: In case two forest polygons are separated more than x meters, e.g. a 100m broad meadow, I would like to keep the two separate forest polygons.

I am using this map for the ecological modeling of animals. This is why habitat patches need to be defined. Basically this process decides if two polygons are still the same habitat or two separate habitats.


texts: I am attempting to merge two shapefiles together. Both are line features. I have merged these files before and it worked successfully however due to problems during editing I have to start again and merge the two files once again. However any time I try to merge this time it comes up as a fail and with the results as shown below.


texts:
I have several large shapefiles from various sources. All cover various sections of oceans and many overlaps, as such the coastline is repeated across the several shapefiles but to varying levels of details/accuracy.

The outcome I would like to achieve is one shapefile, with one coastline.

I can merge the shapefiles in QGIS but how would I align the coastline and drop its duplication?


texts: Why is VBS If Then Label Expression not working?

texts:
I am trying to label a field based on other field.

For
[TYPE] = "PRIMARY", the [LABEL] is having this structure: "TAN80901B"
For
[TYPE] = "SECONDARY", the [LABEL] is like this: "TAN80901C/D20b"

I want the [LABEL] to return "B" if "PRIMARY" and "D20" if "SECONDARY"

To delete from the right and left, I can use : "left(right( [LABEL], 5), 2)"

Any help with this expression?


texts: Measuring along elevated line using ArcGIS Desktop?

texts:
I am trying to measure various lengths along a line which has changing elevation, and splitting the line at the various end points.

I have the lengths that have been taken on along the line (a walking track) but as elevation is constantly changing I have no idea how I am to go about completing this as I haven't been able to find out if measuring along the line will automatically take into account the elevation, or if it will treat it like it is a flat surface which would obviously result in the line being split in the incorrect places.


texts:
I want to clarify something about feature references in WFS 2.0 standard (OGC 09-025r1 document)

Do I understand correctly, that "resolve feature reference" means "get (remote or local) feature by xlink and include it to response results subset"?

Accordingly, remote feature - feature linked by "http://(.*)#myfeature.1", local - feature located on the same (requested) server


texts: I downloaded Landsat image and used composite band toolbox. Now, I produced color map but I need to create boundaries across urban area. Is there any way to create polygons among that area and how?


texts: How to get the rotation of a square Polygon?

texts:
I've got a layer made of perfect square polygons. I would like to obtain the rotation of each polygon:

In this case, the rotation would be 45 degrees.

I suppose I would have to compare it with a 0 rotation polygon, which has all of its sides parallel to the X axis or the Y axis. Example:


texts: PostGIS SQL Editor plugin - path to pgsql2shp

texts:
I'm new to PostGIS/QGIS, and I was looking at using the PostGIS SQL Editor plugin (v.2.0.0) instead of the OpenJump tool recommended in PostGIS in Action 2nd ed. However, I am having problems following the plugin help advice stating that 'The Path to pgsql2shp must be added to the system PATH environment'.

I already have the relevant PostgreSQL paths (;C:\Program Files\PostgreSQL\9.6\bin;C:\PostgreSQL\9.6\lib) added to my system variables (psql works fine from the command line), but I am still getting the error ''pgsql2shp' is not recognized as an internal or external command,
operable program or batch file.' from the QGIS plugin.

I'm using QGIS 2.18.4, PostgreSQL 9.6, PostGIS 2.3.2 on Windows 7 Enterprise.


texts: Best GeoDjango MySQL workaround for distance queries

texts:
I'm using mysql 5.1.50, so the distance QuerySet method is not available to me.

Here is my model:

class Location(models.Model):
    name = models.CharField(max_length='100')
    point = models.PointField(spatial_index=False)

    objects = models.GeoManager()

I would like to find all the locations where the point is within 5 miles of another point (eg. test_point).  I've found some potential ways to do this but am unsure of all of them.

Location.objects.filter(point__within=test_point.buffer(1))

But it's not clear to me what the unit is for '1'.  It looks like degrees.  I can convert miles to degrees but it seems like there's a cleaner way to do this.


texts: How to return the first row of feature class as a list in ArcPy?

texts: I need to return the first record of a feature class as a list. Is this possible? I can't seem to figure it out using row.getValue() as it does not allow me to put a wildcard in there.


texts: Evaluating images in HTML widget using QGIS

texts:
I would like to have images displayed in my HTML widget.
Unfortunately, the initial script generated by QGIS doesn't work.

  Images

It comes with the text string only.
I found some solution at Display photo stored as blob in GPKG but it shows a Python function instead of a front-end solution.
I need my image related to the feature, but not in the HTML display section, but in HTML widget - see
Show images related to features in QGIS?
There are some solutions here but don't work in my case:

https://stackoverflow.com/questions/29910137/how-to-document-write-within-an-image-src-string-doesnt-get-parsed
https://stackoverflow.com/questions/22503908/html-is-img-src-scriptsome-javascript-script-allowed
https://stackoverflow.com/questions/116967/img-src-tags-and-javascript

Is there any way to show in the HTML widget the image related to my object?

Even if I place the link to static image, the problem is the same:

 Images


texts: Seeking Web based GIS for viewing only without need to install client side apps?

texts:
First of all I am new to GIS, I have used ArcGIS and QGIS for some basic engineering projects. Mainly I have looked for freely available GIS data online (solar radiation), downloaded the data, generated JPEG images and submitted engineering reports. I have also done some very basic tasks in ArcGIS, like aspect map, slope map, buffer zone, etc.

Here is my problems, I have a server, i want to log in to server via internet using google Chrome browser and accesses the GIS data that is on the server for viewing only!

I have in total 6 layers, and it would be nice to be able to turn the layers on and off.

Think of it as Google Maps, but only a person with a user/pass can log into it and the person can view a number of layers, turn them on and off, zoom in and out, that is it!

It would be nice if I can make a connection between the above layers and OpenStreetMap as well.

Is this possible and how?

Is there a GIS open source software (or commercial software) that allows for such task?

I have been looking at "web based QGIS" but i am very confused and don't think it is what i want.

Note that I don't want the client to install any apps on his computer, I want him to only use Google Chrome browser without any addons.


texts: Openlayers request header is not going to server

texts:
I have created a Proxy script on my server. This Proxy script is standing between client and geoserver. But Proxy requires Authorize header from user. Openlayers application will send requests to Proxy with Authorize header.

I am setting default headers

OpenLayers.Request.DEFAULT_CONFIG.headers= { Authorization: "BASIC" + "123456" };

But server request header does not contain Authorization header that I set.

        var lines = new OpenLayers.Layer.WMS(
            "My Lines",
            "http://localhost:43561/geoserver/wms",
            { layers: "GIS:my_lines", transparent: true},
            { isBaseLayer: false, layerId: "paths", isFilterEnabled: true }
        );

This is my wms layer.

How can I set header values of OpenLayers request?


texts:
I apologies if this is a very open ended question but I need a point in the right direction. I've got GeoServer downloaded and running nicely on my Tomcat server. I want to do something quite simple and make a track from a number of GPS coordinates to show the path of satellite across the earth and save this into a png/gif/jpeg so it can be used by another application.

I thought you'd sort of upload your background image and then plot the tracks onto it. But I don't think that's correct. I need to have the map as a shapefile too, right?

So what I'd like to do is:

- Have a map of the world
- Read GPS data from mySQL database.
- Plot cooridnates onto map to generate track
- Output track and map as image file.

Is there away to schedule this so it updates the track every few hours?


texts: understand .prj and what CRS is used (QGIS)

texts:
I am trying to devise what CRS is being used by a shapefile published by the Slovenian Ministry of the Environment and Spatial Planning here. Unfortunately, I cannot download the .prj file, but I can open it in the browser (?) here is what I can see:

  PROJCS["Gk-XY",GEOGCS["GCS_Bessel_1841",DATUM["D_Bessel_1841",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Gauss_Kruger"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",-5000000.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9999],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

I copied the above text in notepad and created a .prj file. When I open the shapefile in QGis and look at the info of the layer this is what I get:

When I try to overlay the shapefile to a georeferenced orthophoto with correct CRS (EPSG:3794 - Slovenia 1996 / Slovene National Grid - Projected) I can see that there is a translation:

Can you help me find out what CRS has been used to create the shapefile?


texts: How to filter TIFF images (mosaic/pyramid) in GeoServer?

texts:
I have total 7000 TIFF images that I want to display on application in OpenLayers. I didn't want to create 7000 stores so I created mosaic and hosted it. But how can I filter TIFF images based on extent? I don't have any time or elevation data attached to TIFF image.

Can you also suggest any other approach to tackle this problem? (Hosting 7000 images and displaying them based on filtering either by bbox or anything else.)


texts: Kriging error when writing ArcPy script (ERROR 000622 & ERROR 000628)

texts:
I am trying to write an ArcPy script that uses Kriging to interpolate various nutrients in the soil (i.e. nitrogen) using ArcGIS Pro 2.4.0. These nutrients were collected by sampling the soil, so it will be a very high-resolution output raster.

The Kriging tool will run when using ArcGIS Pro. Here are the parameters:

(Every field below this tool is blank).

However, when running an ArcPy code, get the following error messages:

Traceback (most recent call last):
  File "<pyshell#24>", line 1, in <module>
    outKriging = Kriging(inFeatures, field, kModelOrdinary, cellSize, kRadius, "outVarRaster")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.7\ArcPy\arcpy\sa\Functions.py", line 2804, in Kriging
    out_variance_prediction_raster)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.7\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.7\ArcPy\arcpy\sa\Functions.py", line 2796, in Wrapper
    out_variance_prediction_raster)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.7\ArcPy\arcpy\geoprocessing\_base.py", line 510, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 000622: Failed to execute (Kriging). Parameters are not valid.

ERROR 000628: Cannot set input into parameter semiVariogram_props.

I'm new to Kriging, so what can I do to get the tool to run?

Here's my code:

#Import modules
import arcpy
from arcpy import env
from arcpy.sa import *

#Set workspace
env.workspace = "C:\\Users\\asdf\\Documents"
env.overwriteOutput = True

#Set simple parameters
inFeatures = "ne10soilpnts"
field = "NO3"
outRaster = "krig"
cellSize = 0.000026
lagSize = 0.000026

#Set complex variables
kModelOrdinary = KrigingModelOrdinary("ORDINARY", lagSize)

kRadius = RadiusFixed(12)

#Error occurs after running this line
outKriging = Kriging(inFeatures, field, kModelOrdinary, cellSize, kRadius, outRaster)


texts: How to visualize a C-band SAR sensor measurement

texts:
Files: \measurement\s1a-iw2-slc-vv-20141205t171642-20141205t171705-003587-0043c5-002.tiff
Size is 24435, 13472
Coordinate System is `'
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2014:12:05 21:24:05
  TIFFTAG_IMAGEDESCRIPTION=Sentinel-1A IW SLC L1
  TIFFTAG_SOFTWARE=Sentinel-1 IPF 002.36
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,13472.0)
Upper Right (24435.0,    0.0)
Lower Right (24435.0,13472.0)
Center      (12217.5, 6736.0)
Band 1 Block=24435x1 Type=CInt16, ColorInterp=Gray

Its from teh sentinel 1 sat, and I am wondering about what the steps are to visualize this in somethign that makes sense to normal people. In the download, there is a preview image that makes sense, but I dont know how to create that myself.


texts: LAStools menu seems to be incomplete in QGIS?

texts:
I am using QGIS 2.18.2 Las Palmas on a MacBook Pro running OSX 10.13.1 High Sierra.
I followed Martin Isenburg's instructions how to install LAStools toolbox in QGIS
and the LAStools menu is available under the "Processing > Toolbox" menu as it should be.

However, the number of tools (algorithms?) available to me in the menu is, I think, smaller than it should be. My LAStools installed menu shows 13 tools or menu items; the example in Mr. Isneburg's instruction shows at least 32 items.

Here is a screen shot of what I have:

And here is what Mr. Isenburg's image, from his instructions, shows that I should have:

How do I get more of the menu items, as shown in the tutorial's screenshot?

Update 2017-11-28. I found that some files in the LAStools folder are zipped.

Should I go unzipping these files in order to have a proper installation of LAStools?


texts: I have a series of lat/long coordinates in NAD27 that i need to convert to NAD27 zone 17n eastings/northings. I have ARCGIS 9.3 but no direct and reliable means for the transformation. Any ideas?


texts:
I have ~125 analog aerial photographs from 1937 (examples below) and need to reconstruct their locations where the images were taken for a research project. Unfortunately, little is known about them: neither if they are all taken at the same altitude nor if the images are oriented to the north. They could have been taken anywhere within the Germany's borders at the time. And the images were printed in different sizes and thus differ in print/scan resolution. Image recognition won't work because the spatial structures have changed too much since then, but I would imagine that a large part of the land parcels and streets / squares have remained the same.
I found this paper by Dragos Costea and Marius Leordeanu describing "Aerial image geolocalization from recognition and matching of roads and intersections", and so far this looks to be the most promising approach for my problem, as it doesn't rely on image recognition. It would be a success for me if I can locate about 7-10% of the images.
However, since I am not a specialist and the implementation / testing means a lot of effort for me, I wonder if you think the project is realistic in principle or even know an alternative method that would be worth trying?

Update: With the help of a known subject (Monument to the Battle of the Nations) I can know determine that the images are not oriented to the north.


texts:
I made my first choropleth map using Folium which worked quite well! However I would like to change the opacity of my circles to be a bit more opaque. How does one do this? There does not seem to be an option.
m = folium.Map(location=[48, -102], zoom_start=3)
for geoid, (long, lat) in geoids.items():
    city_name = df_metro.loc[df_metro['geoid']==geoid, 'Geographic Area Name'].iloc[0]
    city_ratio = df_metro.loc[df_metro['geoid']==geoid, 'Ratio'].iloc[0]
    city_pop = df_metro.loc[df_metro['geoid']==geoid, 'Total'].iloc[0]
    folium.Circle(
        location=[lat, long],
        color=colormap(city_ratio),
        radius=50*math.sqrt(city_pop),
        fill=True,
        opacity=0.8,
        tooltip='{}: {}'.format(city_name, city_ratio)
    ).add_to(m)


texts:
I've got some data from "http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp" and it is in *.asc form. To process them I need to convert them to shapefile (need to merge them all into one shapefile and add attribute table, on and on ...), but I don't know how. Some web sites such as "http://converter.mygeodata.eu/ " do convert my data but I'm talking about 2.5 GB of data and it is impossible for me to upload this much data ... how can I convert them inside QGIS or ArcGIS ??? Is there another software that can convert my *.asc file into *.shp ???
A lot of thanks


texts: How to edit filegeodatabase with Python GDAL/OGR bindings?

texts:
I'm able to view layers, fields and iterate features but when I try to delete/create a feature or add a new field I run into errors or it does nothing. I'm using the FileGDB driver. If I throw this .gdb into QGIS I can easily delete, or add features (which confirms FileGDB driver works). Is there a way to do it with Python? Ideally I find a way to add/delete fields and layers and add/delete features.

from osgeo import ogr,gdal

test_gdb = r"C:\Users\52\Downloads\geodatabase\GIS\DdaS69f.gdb"
driver = ogr.GetDriverByName("FileGDB")
dataSource = driver.Open(test_gdb, 0)

fldDef = ogr.FieldDefn('DistFld', ogr.OFTReal)
fldDef2 = ogr.FieldDefn('Name', ogr.OFTString)
fldDef2.SetWidth(16) #16 char string width

for x in range(dataSource.GetLayerCount()):
    layer = dataSource.GetLayerByIndex(x)
    layer.CreateField(fldDef)
    layer.CreateField(fldDef2)
    print("Layer Name: {}".format(layer.GetName())) ## Get layer name
    layerDefinition = layer.GetLayerDefn()

    for field in range(layerDefinition.GetFieldCount()):

        print("\nField : {}\n".format(layerDefinition.GetFieldDefn(field).GetName()))
        for feature in layer:
            print(feature.GetField("{}".format(layerDefinition.GetFieldDefn(field).GetName())))
            layer.DeleteFeature(feature.GetFID())

        layer.ResetReading()

Is it possible to edit a .gdb using the FileGDB driver?


texts: Error "Unable to open data source" - import .shp into GRASS

texts:
I'm new to GRASS and GIS in general, so please forgive my ignorance...

I have problem with loading Vector Data (.shp) into GRASS

I used the following steps 1.1. and 1.2:
http://qgis.org/en/docs/training_manual/grass/grass_setup.html?highlight=grass

In the end I have choice loaded layer (.shp) and when I clicked the "Run" button to begin the import, appear error: Unable to open data source

I have all files like .shx, .shp in the same directory.
I don´t know what I´m doing wrong.


texts: Update database via osmosis and osm2pgsql too slow

texts:
I got a box with the current hardware config:

CPU: Intel Xeon E3-1220 V2 3.1 GHz
RAM: 10gb
Hard Disk: Western Digital WD10EZRX SATA Hard Drive

OS: Centos 6.5 64 bit
DB: PostgreSQL 9.3.2

and PostgreSQL got following config

shared_buffers = 1024MB
work_mem = 256MB
maintenance_work_mem = 4096MB
fsync = off
synchronous_commit = off
checkpoint_segments = 100

I imported the entire OSM Europe data into PGSQL DB with the following command:

osm2pgsql -c -d gps -U gps --cache 8000 --number-processes 4 --slim europe-latest.osm.bz2

and it took 6 days to complete... it's a lot but it's ok...

now i want to keep my DB updated with the combo of osmosis + osm2pgsql cmds:

osmosis --read-replication-interval workingDirectory=/osmosisworkingdir/ --simplify-change --write-xml-change - | \
osm2pgsql --append -d gps -U gps --cache 8000 --number-processes 4 -s -

Well... this takes MORE than one day to just align single day changes! That makes impossible to keep my local DB aligned!

Am i doing something wrong? Should i upgrade my RAM? Why updating is SOOOOOOOOOO damn heavier than load data for the first time?


texts: Creating 9.3 version of geodatabase in ArcMap 10?

texts:
Under the What's new for geodatabases in ArcGIS 10 Esri talks about an optional parameter that allows you to create a 9.3 geodatabase in version 10.

I cannot find any other info.

Any help on how to get to this "optional parameter" working?


texts: How to Convert Code from VB to Python in Map Algebra?

texts:
I have a problem with Map Algebra in ArcMap 10 to convert some codes made by VB, so i need to convert it to python to work in arcmap 10 without problem.
and can any one help me how do i deal with map algebra in arcmap10?

the code:

con([Direction to univer] > 0 & [Direction to univer] <= 180,[Direction to univer] + 180,[Direction to univer] - 180)

Thank you.


texts: ArcGIS Pro Python install causing issues with importing packages

texts:
I've looked at several posts and articles and I can't seem to find a solution to this. I need to use jenkspy for a specific task but my Python 3x was installed with a software called ArcGIS Pro that designates the default Python environment "Read Only" meaning that packages not installed with it must be installed in a virtual environment.
So, I opened command line and ran:
py -m venv env
py -m pip install jenkspy --user

If I continue to run python from command window I can import jenkspy no problem. However, when I attempt to import it in iPython or load it into pythonWin as a script I get:
import jenkspy
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ModuleNotFoundError: No module named 'jenkspy'

So I go back to command window and run python to locate path for jenkspy:
import sys
print (sys.path)

which returns:
['', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\Lib', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\Lib\\idlelib', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\python37.zip', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\DLLs', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\lib', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64', 'C:\\Users\\johnsmith\\AppData\\Roaming\\Python\\Python37\\site-packages', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\lib\\site-packages']

I navigate to the native idle.py file in this environment and import jenkspy. No problem. I print sys.path from the idle console. Same result as in command line.
I then go back to iPython and run:
import sys
paths = ['C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\Lib', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\Lib\\idlelib', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\python37.zip', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\DLLs', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\lib', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64', 'C:\\Users\\johnsmith\\AppData\\Roaming\\Python\\Python37\\site-packages', 'C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Python37_64\\lib\\site-packages']

for p in paths:
    sys.path.append(p)

and then:
import jenkspy

And I get this error:
ImportError                               Traceback (most recent call last)
<ipython-input-19-88ab181dc114> in <module>
----> 1 import jenkspy

~\AppData\Roaming\Python\Python37\site-packages\jenkspy\__init__.py in <module>
      3 __version__ = "0.1.5"
      4
----> 5 from .core import jenks_breaks
      6
      7

~\AppData\Roaming\Python\Python37\site-packages\jenkspy\core.py in <module>
      7     from math import isinf, isnan
      8     isfinite = lambda x: not isinf(x) and not isnan(x)
----> 9 from . import jenks
     10 try:
     11     import numpy as np

ImportError: cannot import name 'jenks'

Tried to do conda install. Works fine in iPython if I open the venv jupyter notebook. Still, not seeing how this will work in a stand alone script which I need it to do because this will be a batch job run daily.
from anaconda prompt these are available envs:
base                  *  C:\Users\reguser\AppData\Local\Continuum\anaconda2
agolupdate               C:\Users\reguser\AppData\Local\Continuum\anaconda2\envs\agolupdate

This is from print (sys.path) in ArcGIS Pro python window:
C:\\Users\\adm.user\\AppData\\Local\\ESRI\\conda\\envs\\arcgispro-py3-clone2

What do I need to do in order to get jenkspy or other post Python installation packages/modules to be imported using stand alone scripts and or/iPython?


texts: Mapillary Vector Tiles are not shown at every zoom levels in a Leaflet map

texts:
I'm trying to load the Mapillary vector tiles layer in a Leaflet map.

I'm using Leaflet.MapboxVectorTile (ref.  https://github.com/SpatialServer/Leaflet.MapboxVectorTile).

Here you are my code ...

  var map = new L.map('map').setView([44.907852, 7.673789],16);

  L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
    attribution: '? OpenStreetMap contributors'
  }).addTo(map);

  var config = {
    url: "https://d2munx5tg0hw47.cloudfront.net/tiles/{z}/{x}/{y}.mapbox"
  };

  var mapillarySource = new L.TileLayer.MVTSource(config);
  map.addLayer(mapillarySource);

The code is working and I can see the vector tile layer on the map but when I zoom in on my map I can't see the layer after 14 zoom level.

I've tried to use also

minzoom = 0
maxzoom = 16

in the layer configuration but nothing change.

Trying to use the Mapillary viewer, for example

http://www.mapillary.com/map/search/44.90213653940583/7.661219012127532/13

I see that the vector layers at every zoom level.

Using Firebug I see that their code is quite similar to my code  .....

      config = {
        url: 'https://d2munx5tg0hw47.cloudfront.net/tiles/{z}/{x}/{y}.mapbox?access_token=pk.eyJ1IjoibWFwaWxsYXJ5IiwiYSI6ImZCRUVhYm8ifQ.MwcTuq2Q0KHcJxC_Whra7w',
        zIndex: 1000,
        style: function(feature) {
          var style;
          style = {};
          style.color = 'rgba(53, 175, 109, 0.5)';
          style.size = 3;
          return style;
        }
      };
      mvtSource = new L.TileLayer.MVTSource(config);
      map = mapCtrl.getMap();
      return map.addLayer(mvtSource);

For test I've tried to change my code using the Mapillary original layer definition above (with  and without the access_token ....), but nothing change.

How may I see the Mapillary vector tiles at each zoom level?

Any suggestion, example, jsfiddle?


texts: There used to be a plugin that allowed GIS data (.shp file) to be exported and then opened in SketchUp.  This plugin no longer exists.  I have wireframe data of buildings that is 3D in ArcGIS that I need to be able to import into SketchUp so that I can then build the model from the wireframe.  Does anyone have any idea how to go about doing this?  I have ArcGIS 10.3 and SketchUp 2015.


texts: Using QGIS modules in standalone application with PySide?

texts:
I want to make small standalone application using PySide and QGIS Python modules to display satellite images and do standard GIS operations, like creating new a shapefile and placing polygons or points inside it.

I installed QGIS 2.2 and added paths to the Python folder and the QGIS .dll's.

PYTHONPATH: I added ;C:\Program Files\QGIS_2.2\apps\qgis\python; - path to python folder inside QGIS.

path: ;C:\Program Files\QGIS_2.2\apps\qgis\bin; - path to .dll's like qgis_core.dll and qgis_gui.dll.

When I try to import qgis.core, I get following:

Traceback (most recent call last):
  File "F:/Downlods/PyCharm_Projects/Project_1/Testing_Area.py", line 1, in <module>
    import qgis.core
  File "C:\Program Files\QGIS_2.2\apps\qgis\python\qgis\__init__.py", line 26, in <module>
    import sip
ImportError: No module named 'sip'

So, its telling me there is no 'sip'. I don't know, is there is some problem with my 'PYTHONPATH' or 'path'? Or its just will not work with PySide?


texts: I want to import  AURA-OMI ozone level3  data which is available from http://giovanni.gsfc.nasa.gov/giovanni/ .The data is in hdf5 format. The extension is .he5. Do anyone know how to import in grass. Is there any tool in gdal for importing to geo-tiff or bil.


texts: Find the minimum elevation pixels from an elevation raster using arcpy

texts:
I want to create a python script  to find the minimum elevation pixels from an elevation raster (aster,srtm) using arcpy and arcmap 10.1.  The purpose is to reclassify these pixels.
I have tried but I do not have enough knowledge in scripting.

Can somebody please help me ?


texts: Using Mapnik's OGR plugin for Oracle Spatial Database

texts:
Mapnik wiki page of OGR plugin says that I can use the plugin for MSSQL server in the following manner:

<Layer name="mssql" srs="+proj=latlong +datum=WGS84">
    <StyleName>mssql_style</StyleName>
    <Datasource>
        <!-- mssql database must contain geometry_columns and spatial_ref_sys metadata tables-->
        <Parameter name="type">ogr</Parameter>
        <Parameter name="string">MSSQL:server=localhost;database=gis;trusted_connection=yes</Parameter>
        <Parameter name="layer_by_sql">SELECT * FROM dbo.planet_osm_line</Parameter>
        <Parameter name="extent">-180,-90,180,89.99</Parameter>
    </Datasource>
</Layer>

I want to use it for Oracle Spatial Database. How should I modify the parameter with name = "string" in the XML so that I can?

Also, do I need to install any additional software and make any configurations on my machine so that I can connect to the database?

By the way, the database is not on my localhost but in a remote location.


texts:
I would like to make a buffer for a polyline and get as Output a Polygon with the same angles as the polyline, having not arcs, like in ordinary buffers.

i tried the FLAT option in ArcPy, but it's applied to just the both ends.


texts: Densify tool equivalent for Basic ArcGIS License?

texts: I'm looking for a tool/script/quicker process that does the same thing as the standard license's densify tool. I'm only using a basic license though. Just wondering if there was anything out there that did something similar.


texts: Clipping shapefile to interpolated surface (contours) in R?

texts:
R novice here, trying to work through the ever-daunting task of clipping my data. I've posted several questions regarding how to clip various types of data, so here we go with another, complete with all of my data to reproduce results. :)

Data:
CSV data,
Shapefile

My map is of rainfall return periods in Tanzania. I used the IDW interpolation method to achieve this plot. Instead of having the entire surface over/around TZ contoured, I only need the contours to be clipped to the TZ shapefile. Here's my code, which was provided by this excellent tutorial:

setwd("D:/Mapping-R/Returns-Practice")

library(ggplot2)
library(gstat)
library(sp)
library(maptools)

tz <- read.csv("TZ_ReturnPeriods.csv", header = TRUE)

temp <- tz
temp$x <- temp$Lon
temp$y <- temp$Lat

coordinates(temp) <- ~x + y

plot(temp)

x.range <- as.numeric(c(28.61, 41.19))
y.range <- as.numeric(c(-12.61, -0.19))

tz.grid <- expand.grid(x = seq(from = x.range[1],
                               to = x.range[2],
                               by = 0.1),
                       y = seq(from = y.range[1],
                               to = y.range[2],
                               by = 0.1))
coordinates(tz.grid) <- ~x + y

plot(tz.grid, cex = 1.5, col = "grey")
points(temp, pch = 16, col = "blue", cex = 0.5)

idw <- idw(formula = Return.10 ~ 1,
           locations = temp,
           newdata = tz.grid)
idw.output <- as.data.frame(idw)
names(idw.output)[1:3] <- c("Longitude", "Latitude", "Return")

ggplot() + geom_tile(data = idw.output,
                     aes(x = Longitude,
                         y = Latitude,
                         fill = Return)) +
  geom_point(data = tz,
             aes(x = Lon, y = Lat),
             shape = 21,
             color = "red")

tz.contour <- readShapePoly("TZA_adm0.shp")
tz.contour <- fortify(tz.contour, region = "NAME_0")

ggplot() + geom_tile(data = idw.output,
                     alpha = 0.8,
                     aes(x = Longitude,
                         y = Latitude,
                         fill = round(Return, 0))) +
  scale_fill_gradient(low = "cyan", high = "orange") +
  geom_path(data = tz.contour, aes(long, lat, group = group), color = "grey") +
  geom_point(data = tz, aes(x = Lon, y = Lat), shape = 16, cex = 0.5, color = "red") +
  labs(fill = "Rainfall", title = "10 Year Return Period in Tanzania")

This code produces this image:

BUT that was only half of the battle. I need to clip the Tanzania shapefile to the contoured surface so that ONLY the contours within Tanzania's borders are shown. I've found many almost-solutions, but nothing that quite works.


texts: Can Google Maps Android API display WMS layer?

texts:
I am currently making an app for my final career project, and I am trying to combine the layer of Google Maps and WMS layer that displays a map from a WMS server, helping me with ArcGis SDK for Android,

The MapView object of the ArcGIS API for Android allows you to add WMS layers (with the object WMSLayer), while Google Maps Android API v2 does not allow it, and I can't find a way of combining these layers.

Someone could tell me if this is possible?


texts:
I have created a folium map with markers and popups in jupyter notebook and saved it as map.html
I am looking for a way to share this map with others or host on the web, such that the markers and popups show. Any help?


texts: Is there a way to find the size (in bytes) of the header information in a GeoTIFF file using Python? I am not trying to read the header (I know I can use gdalinfo for this) but rather to figure out what position I need to skip to in order to read the file as binary. This is because I am interested in testing numpy memory-mapped arrays for reading portions of an image from the disk. I would like to test this versus the gdal builtin ReadAsArray() method because I enjoy the flexibility of numpy indexing.


texts: OpenLayers 3, PostgreSQL, WFS: Extent jumps back to initial extent if zoom out too far

texts:
I'm working on an app using OpenLayers 3 with PostgreSQL/PostGIS tables and views, and an Esri basemap.

The web app loads a set of features using an array of attributes from PostGres table (a view) and a cql filter. then zooms to the extent of the collected features.
All this works as expected.  The problem happens when I zoom out after the features are loaded...if I zoom out more than a few levels, the extent suddenly jumps back to the initial extent of the features when loaded.

Here is (what I think is) the relevant code below... I am not getting any errors in console; the fetching and loading of the features looks ok from what I can see.

var map = new ol.Map({
  layers: [
  new ol.layer.Tile({
    source: new ol.source.TileArcGISRest({
      url: 'http://server:6080/arcgis/rest/services/basemap/MapServer',
  })
})
],
  overlays: [overlay],
  target: document.getElementById('map'),
  view: new ol.View({
    center: ol.proj.fromLonLat([lon, lat]),
    minZoom: 11,
    zoom: 15,
    maxZoom: 20
  })
});

...

var locExtent = ol.extent.createEmpty();
layerWFS = new ol.layer.Vector({
      source: new ol.source.Vector({
          loader: function (extent) {
              $.ajax('http://geoserver:8080/geoserver/wfs', {
                  type: 'GET',
                  data: {
                      service: 'WFS',
                      version: '1.1.0',
                      request: 'GetFeature',
                      typename: 'all_features_search',
                      srsname: 'EPSG:3857',
                      outputFormat: 'application/json',
                      CQL_FILTER: "attribute IN " + tostring
                  }
              }).done(function (response) {
                  layerWFS
                  .getSource()
                  .addFeatures(new ol.format.GeoJSON()
                  .readFeatures(response));
                  layerWFS.set('name', 'search_layer');
                  //
                  var radius = 100;
                  var  locExtent = layerWFS.getSource().getExtent();
                  locExtent = ol.extent.buffer(locExtent, radius);
                  ol.proj.get('EPSG:4326'), ol.proj.get('EPSG:3857'));
                  map.getView().fit(locExtent, map.getSize());
                  var theseFeatures = layerWFS.getSource().getFeatures();
                  for (var i = 0; i < theseFeatures.length; i++) {
                    var featureWKT = theseFeatures[i].get('wkt');
                    var featureType =
                    featureWKT.substr(0,featureWKT.indexOf('('));
                    // there are both point and line features in view
                    if (featureType == "POINT"){
                      theseFeatures[i].setStyle(pointStyle);
                    } else { // MULTILINESTRING
                      theseFeatures[i].setStyle(lineStyle);
                    }

                  }
              });
          },
          strategy: ol.loadingstrategy.bbox,
          projection: 'EPSG:3857',
          crossOrigin: 'anonymous',

      })
    });

  map.addLayer(layerWFS);


texts: Visualize feature class with many relationship classes?

texts:
I have a feature class in a geodatabase that I want to review.

It has many relationship classes - possibly a hundred or more.  If I load the feature class into ArcMap and identify one of the features, it takes several minutes to expand the relations.

What is a good way to visualize this feature class so that I can understand it?  Specifically, I want to know if having a relation to a relation will obstruct editing the feature class.


texts:
Is it possible to create a style file to complement a shapefile using ArcEngine?

I'm using C# VS2010 ArcEngine 10.


texts: Export CSV in Google Earth Engine from multiple point

texts:
I want to export the csv file which will contain the data associated with multiple points from the following code. Here I am trying to export one csv file which will contain the datewise  data of the different points. But the exported csv shows the data for a single point not for multiple points. Please suggest me what to do.
// Make a list of Features.
var features = [
  ee.Feature(ee.Geometry.Point(93.881640625, 37.391718430950796), {name: 'city1'}),
  ee.Feature(ee.Geometry.Point(-106.27421875, 33.71035637052483), {name: 'city2'}),
  ee.Feature(ee.Geometry.Point(-117.436328125, 39.7618477958936), {name: 'city3'}),
  ee.Feature(ee.Geometry.Point(-95.5515625, 29.74978297462196), {name: 'city4'}),
  ee.Feature(ee.Geometry.Point(-115.151171875, 36.15978520358036), {name: 'city5'}),
];

// Create a FeatureCollection from the list and print it.
var fc = ee.FeatureCollection(features);
print(fc);

// get the data
var y2019 = ee.ImageCollection("COPERNICUS/S5P/NRTI/L3_NO2")
              .filterDate("2019-01-01","2019-01-30");
// Select the NO2 column number density.
var times2019 = y2019
    .select("NO2_column_number_density")
// Create a time series chart.
var TimeSeries = ui.Chart.image.seriesByRegion(
    times2019, fc, ee.Reducer.mean(), "NO2_column_number_density", 200, 'system:time_start', 'label')
        .setChartType('ScatterChart')
        .setOptions({
          interpolateNulls: true,
          title: 'NO2 Timeseries',
          vAxis: {title: 'NO2'},
          lineWidth: 1,
          pointSize: 4,
          series: {
            1: {color: '0000FF'}
}});

// Display.
print(TimeSeries);

// filter to region
var filtCol = times2019.filterBounds(fc);

// export to CSV
var data = ee.FeatureCollection(filtCol.map(function(image){
  return ee.Feature(null, image.reduceRegion(ee.Reducer.mean(), fc, 200))
             .set('system:time_start', image.get('system:time_start'));
}));
print(data)

// export directly as CSV
Export.table.toDrive({collection: data,
                      description: 'YOURDESCRIPTION',
                      fileFormat: 'CSV',
                      selectors: ['NO2_column_number_density', 'system:time_start']})


texts: Using a length attribute and a start coordinate to get the end coordinate on a vector in FME

texts:
I have a road network (with vectors all in correct direction of travel).
I  have a number of coordinate pairs (call these the start lat longs) that are found on the network.
Each start lat long has an assigned Lenght attribute (ranging from 2m - 300m),
What I'm having trouble with is locating an end lat long coordinate pair (on the network) using the Length attribute.
I know how to get the length to point and extract coordinates, but is there someway in fme to combine the Length to point transformer and add the Length attribute to then be able to extract the end lat longs?


texts: Showing annotations in exported map in QGIS 3.4

texts:
I have a problem to export a map with annotations in QGIS 3.4.
I create annotations in the visor, but when exporting the map, this shows with other size.
What is the problem? How can I solve it?


texts:
For example I have some large TIFF file, then I build tiles from it using gdal2tiles, then I need to share this tiles with my colleagues.

One of the options is to comress tiles to .zip archive and send it to colleagues, for example here is solution how to view tiles in Google Earth, but I need solution without any additional software on client side, so just browser without plugins allowed.

So my question is there any free service where I can publish my maps? If not, I assume I need to rollout some kind of server on cloud patform by myself, any suggestions?


texts: How do I load a CBG GeoJSON into Python for Census Block Group level maps in Plotly?

texts:
I made county level maps a while ago and I'm trying to modify it to CBG level.
Am I limited to state and county?
with open('cbg.geojson') as response:
    cbgfile = json.load(response)

data = pd.read_csv("roundone.csv")

fig = px.choropleth(data, geojson = cbgfile, locations='FIPS',
                        scope = 'usa',
                        color = "bids_for_item",
                        #color_continuous_scale = ["white", "lime", "green"]
                        )

fig.show()

I'm having issues loading the JSON (invalid argument). According to some other posts it might be because the GeoJSON is 3 GB. How could I get around that?

File
"/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/geojson/codec.py",
line 43, in load
**kwargs)   File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/json/init.py",
line 296, in load
return loads(fp.read(), OSError: [Errno 22] Invalid argument


texts:
I wonder what the unit is of the float being returned from ST_Distance.

In the documentation it says:

  ...cartesian minimum distance (based on spatial ref) between two geometries in projected units.

What are these projected units?

The geometry is stored in a field: geometry(Point,4326).


texts: Calculating area of geometry using Geopandas?

texts:
Is there an easy way in geopandas to calcuate polygon areas after an intersection analysis?
I am switching from ArcGIS to open source and try to translate my code.
In ArcGIS I used:

arcpy.AddGeometryAttributes_management("df1", "AREA", Area_Unit="HECTARES")


texts: Unable to calculate values in a joined table

texts:
Using Arc 10.4.1, I have joined a geodatabase feature class table to a standalone geodatabase attribute table.  This is a many:1 join (standalone attribute table:feature class).  After making the join, I need to do a series of queries on the standalone attribute table, then calculate a field in the feature class table.  Because of this many:1 situation, I am forced to join the feature class table to the standalone attribute table, rather than joining the standalone attribute table to the feature class table.

However, when I try to change values in the joined (feature class) table, the field calculator is greyed out. Sure enough, the Arc Help documentation states:

  "You can only calculate values for fields in the origin table."

But I could swear that I've successfully used the field calculator on a joined table in the past (Arc 9.x, methinks).

Is there a way to create calculations on joined fields?


texts: QGIS / pyqgis error: wrapped C/C++ object of type QProgressBar has been deleted See log for more details

texts:
I have this error when I run some of my python scripts.  It does not always occur.  It always happens at the END of the script, (all script actions finish!!)

I mainly need to fix it because it is annoying and it prevents me from stringing separate scripts together.

The traceback I get is:

Traceback (most recent call last):
          File "C:/QGIS/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py", line 203, in execute
            self.processAlgorithm(progress)
          File "C:/QGIS/apps/qgis/./python/plugins\processing\script\ScriptAlgorithm.py", line 378, in processAlgorithm
            exec((script), ns)
          File "<string>", line 146, in <module>
        RuntimeError: wrapped C/C++ object of type QProgressBar has been deleted

However the line it refers to does not seem to make any sense:

error2 = QgsVectorFileWriter.writeAsVectorFormat(bufflayer, buffdef, "CP1250", None, "ESRI Shapefile")

Note that this and all of the python after this command always run just fine.
Any thoughts?


texts: Does MOD021km images from MODIS needs atmospheric corrections

texts:
I have downloaded some mod021km modis images and I want to work with them.

I want to know if these data need atmospheric corrections?

Perhaps something like FLASSH in ENVI or other methods.

I can convert hdf files using modis conversation toolkit and convert them to Radiance,Reflectance or brightness temperature but I don't know if it needs any atmospheric corrections before converting.


texts: How to fix a broken plugin (it cannot import name NavigationToolbar2QTAgg)

texts:
I want to use the the pluging ProfileFromPoints but I am getting an error saying cannot import name NavigationToolbar2QTAgg. The solution seems to link with an old library name as mentioned here: https://github.com/semiautomaticgit/SemiAutomaticClassificationPlugin/issues/2

Can I amend it without wayting for the bug fix?


texts: ArcGIS - if clause in the field editor with python

texts:
I have 2 columns in my table.
"dB1" and "dB2"
Now I want to calculate a third columnn "dB3" depending of the entries in column "dB1" and "dB2"

if "dB1" or "dB2" have an entry of "0" I want to calculate: "dB1" + "dB2"

if neither "dB1" nor "dB2" have an entry of "0" I want to calculate: 10*log(10^(0,1*"dB1")+10*(0,1*dB2))

I struggle with the syntax and the if clause here...
Can someone give me the code as I theoretically know how it should be done but can't work around the syntax in python...


texts: Calculating habitat openness from digital elevation models (DEM)?

texts:
I have digital elevation models (DEM) of an area and want to calculate habitat openness.

Can anyone help calculating habitat openness from a digital elevation models (DEM)?

I tried the Package 'horizon' - cran.r (R: The R Project for Statistical Computing), but got some strange results. I'm not familiar with this type of analysis, hence I do not know the scale. I'm about to analyse habitat selection in a owl inhabiting coastal and maritime heath habitats. I have GPS-positions and want to check if the owl are more likely to stay at elevations in the terrain (place with a good overview of the landscape). The owl has a home-range of about 30 to 50 km2, the accuracy of each owl GPS position is about +/- 25 m

The two DEM files could be downloaded here: https://www.dropbox.com/sh/ct02cns54c9i1wo/AAAaAT-h-F8L-FA8npm6flFoa?dl=0


texts: Using two different reference scales in one MXD?

texts:
I'm putting two MXD's together. They don't have the same reference scale, so I was wondering if it is possible to have to different reference scale in one MXD?

I'm working with ArcMap 10.5.


texts: Leaflet.js and CartoDB.js giving event conflict

texts:
I have a layer from a CartoDB viz that I’m listening to the featureclick event. When a feature is clicked I highlight it by applying cartoCss. However, I want to unhighlight it when I click on an area of map with no feature. I did this using the Leaflet map.on(‘preclick’) event. I chose preclick so it fires before the featureclick and applies the default cartoCss first.

Problem is the two events don’t play nicely. The map click event always fires, but the layer featureclick event only fires some of the time. So it sometimes takes two or three clicks before the feature is highlighted – even when the cursor has changed to a pointer so the map knows there is a feature there.

What is the best way to unhighlight a feature after a mouse click on the basemap?

Can I ignore my map click event when a layer feature is clicked?

Can I change which event fires first i.e. layer featureclick before map click?

Code:

   map.on('preclick', function(e) {
          //remove project highlighting and close popup
          J.Dashboard.highlightProject(null);
   });

layer.getSubLayer(0).on('featureClick', function(e, latLong, pos, data, layerIndex) {
          J.Dashboard.highlightProject(data.name);
   });

   highlightProject : function(projectName) {
          //Update the style
          J.Map.projectsLayer.getSubLayer(0).setCartoCSS(J.Style.getSelectStyle(projectName));

}


texts:
Im just wanting to understand my options for producing an offline map using Mapbox gl js. I have investigated several option with little signs of progress.

First, I have tried using Studio by uploading geojson and styling and then downloading said style as json meeting Mapbox GL Style Specification.

I am unsure of what to do with this json?

Secondly, without success I have tried supplying the original geojson to the data parameter and as expected receive cors errors.

map.on('load', function () {
    map.addSource("markers", {
        "type": "geojson",
        "data": 'wpi.geojson'
    });

    map.addLayer({
        "id": "markers",
        "type": "symbol",
        "source": "markers",
        "layout": {
            "icon-image": "{marker-symbol}-15",
            "text-field": "{title}",
            "text-font": ["Open Sans Semibold", "Arial Unicode MS Bold"],
            "text-offset": [0, 0.6],
            "text-anchor": "top"
        }
    });
});

Error

XMLHttpRequest cannot load file:///C:/Project/wpi.geojson. Cross origin requests are only supported for protocol schemes: http, data, chrome, chrome-extension, https, chrome-extension-resource.

Is this possible with the Mapbox.gl API?

It seems as tho the documentation is limited in directions for this.


texts: is there a function which closes a linestring to a linering? I couldn't found such a function in the documentation.


texts: random forest image classfication on python

texts:
I am new to python, I would like to do a rf classification on an multispectral image which I applied the PCA.
After applying acp on different bands including NDVI I got negative values, after that my training file contains negative spectral values can this be correct?
import pandas as pd
df = pd.read_csv("training.csv")
print("Shape: ", df.shape)
print(df.columns)
print(df.head())
df_cat = pd.DataFrame ()
for i in list(df.columns):
    df_cat['{}_cat'.format(i)] = df[i].astype('category').copy()
    df_cat['{}_cat'.format(i)] = df_cat['{}_cat'.format(i)].cat.codes
    print(df_cat.head())

    X = df_cat.drop('class', axis = 1)
    y = df_cat['class']

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import KFold
kf = KFold(n_splits=2, random_state = 24)
results = []

for train_index, test_index in kf.split(X):
     X_train, X_test = X[train_index], X[test_index]
     y_train, y_test = y[train_index], y[test_index]

     for train_index, test_index in kf.split(X):
         model = RandomForestClassifier(n_estimators = 100, random_state = 24)
     model.fit(X_train, y_train)
     y_pred = model.predict(X_test)
     results.append (f1_score(y_test, y_pred))
     print("Accuracy: ", np.mean(results))

error:

key=key, axis=self.obj._get_axis_name(axis)))

KeyError:

u"None of

here is my training sample:
| Class | ID | X | Y |  Map Y  |  Map Y   |   Lat   |   Lon   |   B1   |   B2   |  B3   |
|-------|----|---|---|---------|----------|---------|---------|--------|--------|-------|
|      1|   1|973|967|713458.25|3953393.25|35.701378|-0.640678| 99.8778|-27.7799|18.7059|
|      1|   2|974|967|713468.25|3953393.25|35.701376|-0.640568| -8.6351|-79.4471|46.5504|
|      1|   3|975|967|713478.25|3953393.25|35.701374|-0.640457| 37.6848|-59.7302|36.0935|
|      1|   4|972|967|713448.25|3953393.25|35.701380|-0.640788|153.5402| -8.8041| 8.3142|
|      1|   5|972|968|713448.25|3953383.25|35.701290|-0.640791|-28.4907|-71.9933|43.7319|
|      1|   6|973|968|713458.25|3953383.25|35.701288|-0.640681|-20.4926|-72.5833|43.9850|
|      1|   7|974|968|713468.25|3953383.25|35.701286|-0.640570| 51.7772|-52.6125|32.2387|


texts:
I want to make a tool in java that can configure import data from OpenStreetMap to PostGIS using Imposm3.

Can I do this in windows?


texts: Python/GDAL--Handling Big Rasters and avoid MemoryErrors?

texts:
I am using python 2.7.10.
I am trying to read a big Raster with gdal and numpy array.
But it show In MemoreErrors when i use this:
band. ReadAsArray()

Note: this error occurs when input raster is very big not for all of the rasters

Here is my Code:

    import sys, os, gdal, struct
from osgeo import osr
from gdalconst import *
import numpy
import time
import subprocess

startTime=time.time()
gdal.UseExceptions()
##OutCellSize = 1000
##OutFormat   = "GTiff"

# Command line arguments
# change to known values if you prefer
InRas= r"F:\Geodaten\Eigene\Luftbild\Orthobilder\20120502\Aschhorn_mosaik.jpg"

src_ds=gdal.Open(InRas)
if src_ds is None:
    print 'Unable to open %s' % InRas
    sys.exit(1)
print "gdal can open the file"

srcband = src_ds.GetRasterBand(1)

try:
    bandArray = srcband.ReadAsArray()
    print bandArray
except:
    print "cannot convert raster to array.\n" "Need to compress raster"

# obtaining reference info from input raster
geotransform = src_ds.GetGeoTransform()
originX = geotransform[0]
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = src_ds.RasterXSize
rows = src_ds.RasterYSize       ###we are here
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(src_ds.GetProjectionRef())

# creating raster in mepory  as outraster from input
try:
    dst_ds = gdal.GetDriverByName("MEM").Create('', src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Byte)
    dst_ds.SetGeoTransform(geotransform)

    # Create for target raster the same projection as for the value raster
    dst_ds.SetProjection(raster_srs.ExportToWkt())
except:
    print "Unable to create in memory raster"

How could I solve this issue. Other wise I am not able to do further analysis.


texts: Why do my font choices not appear when setting symbology in ArcGIS

texts: I have noticed this, it seems to be a bug.  I will set the typeface for the symbology for one legend/feature class, but as soon as I want to apply it to another feature class the font grays out or is not available to use.  However, if I restart ArcGIS, it works fine.  The problem was initially only occurring for a downloaded font, now it is happening with system fonts also.  I am using ArcGIS 10.2.


texts:
I wrote a little tool in ArcObjects/C# for planarizing a line feature class. I am using IFeatureConstruction::PlanarizeLinesFromCursor because I need to apply a query filter.
Now, I figured out a bug: when executing the operation on ArcSDE featureclass the method adds the planarized version of the geometries, but does not delete the origin (unsplitted) features. The effect only appears on ArcSDE feature classes when the filter contains more than 1000 features.

Does anyone know about this bug? Is there any workaround?

I am using ArcGIS 9.3.1 SP 2 and .Net 3.5 on Win7.


texts:
I am trying to build an HTML legend for my web application but I am having an issue with the tag [leg_icon] not being recognized. Seems like MapServer is not substituting any URL for this tag. Other tags like [leg_class_name] and [leg_layer_name] work fine. Has anyone encountered similar issue or know a possible fix?

The legend.html template:

[leg_class_html opt_flag=15]
    [leg_class_name] [leg_layer_name] [leg_class_index]

[/leg_class_html]

And the result of MODE=LEGEND request from MapServer. As you can see the icons for each of the classes are not loading.


texts: Converting double to date field using QGIS?

texts:
In ArcMap I have a table with a column with the date type yyyy-mm-dd, but after opening SQlite in QGIS, this column is written as a double type.

How do I convert data in the form of a number to a date in QGIS?


texts: pgrouting - multiple travelling salesmen?

texts:
I've trying to solve a routing problem to optomise the use of multiple pieces of machinery.

Part A - I have two machines, and thirty distributed sites.  Ten sites require machine A, the other twenty sites can use A or B. A and B can work at the same time. The desired result is shortest path to touch all points once with the correct machine.

Part B - To build on the initial problem, B sites will take twice as long to complete as A sites.  Can I solve the problem to minimum total duration?

I have basicalyl no experience in these types of problems, other than some playing around with pgrouting.  All data is currently stored in PostGIS if it has any relevence, but I'm open to all options


texts: Getting an empty geometry when retrieving the envelope of a GroupElement in ArcObjects: what might be wrong?

texts:
I am working in ArcGIS 10.0 (ArcMap) and Python 2.6.5 and in the Data View. I have an IGroupElement3 to which I add a TextElement and up to several RectangleElements. I want then to retrieve the envelope of the GroupElement, recenter that envelope, and set the envelope as the new geometry of the GroupElement. The intent is to reposition (translate) the GroupElement. Here is a snippet:

pGrpElem = CType(pGroupElementFull, esriCarto.IElement)
pEnv = pGrpElem.Geometry.Envelope
pEnv.CenterAt(pCtrPoint)
pGrpElem.Geometry = pEnv

When the code runs, I get an error on the CenterAt method, saying I am performing an operation on an empty geometry. I checked, and indeed the pEnv is empty. So I have a few questions:

(1) Is this a valid method for moving a GroupElement?

(2) If not, how should I be doing it?

(3) If so, then I have a problem elsewhere. Do I have to add the component elements - and the GroupElement for that matter - to the display's GraphicsContainer to get them to "gel" and have valid envelopes? I was trying to avoid the GraphicsContainer step, because I really don't want the graphics to draw to the screen. The GroupElement is really for loading into an annotation feature. When I use the GraphicsContainer the screen drawing can get messed up unnecessarily. The graphics draw anyway when they are added to annotation, so I don't need two copies onscreen.

(4) If, as described above, I need to do something in addition to simply adding elements to the GroupElement to make their geometries visible to the GroupElement, how might I do that and avoid drawing the graphics to the screen at the same time?

Thanks for the help, all.


texts: vector layer, no more editing is allowed

texts: I have a polygon vector layer with already 529 objects but without any attributes so far and now I get the message that shapefiles without attributes are treated as write protected. What can I do to go on working with my shapefile?


texts:
Is it possible to add a new SRID to sql server 2012?

I've tried

insert into sys.spatial_reference_systems values (4938, 'GEOCCS["GDA94",DATUM["Geocentric Datum of Australia 1994",SPHEROID["GRS 1980",6378137.0,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0.0,0.0,0.0,0.0,0.0,0.0,0.0],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["m",1.0],AXIS["Geocentric X",OTHER],AXIS["Geocentric Y",EAST],AXIS["Geocentric Z",NORTH],AUTHORITY["EPSG","4938"]]', 'metre', 1)

But I get

  Ad hoc updates to system catalogs are not allowed.

Haven't had much luck on google with any answers as well


texts: Changing the ArcGIS Server 10.3 context root

texts:
I have a arcgis 10.3 server running on a machine (dev)
and i want to something like following

http://server1/devarcgis/rest....... (dev)

I have changed the the context root of arcgis apps (dev) in

C:\Program Files\ArcGIS\Server\framework\runtime\tomcat\conf\server.xml

file from arcgis to devarcgis but its not working and its context root is still arcgis instead of devarcgis. Any idea what might be wrong here?

EDIT:

I edited C:\Program Files\ArcGIS\Server\framework\runtime\tomcat\contexts\rootapp\WEB-INF\urlrewriter.xml file and added a rule

<rule enabled="true">
    <from>^/devarcgis/rest(/*)(.*)$</from>
    <to context="arcgis/rest" last="true">/$2</to>
</rule>
<rule>
    <from>^/devarcgis/rest/services(/*)(.*)$</from>
    <to context="arcgis/rest/services" last="true">/$2</to>
</rule>

now  a call to

http://localhost/devarcgis/rest/services/myservice/MapServer/export?f=json......

returns

ags_jsonp.ags_jsonp_7_48333({"href":"http://localhost/arcgis/rest/directories/arcgisoutput/myservice/_ags_map784f1c72c2a24befa3e0519617d2ef80.png................

in response.


texts: Error using con in Map Algebra in Python script

texts:
I am experiencing a painful problem in the middle of a Python script I am writing. When I get to a Map Algebra function, I have tried two different options and both return a different error.

At first I wrote:

facteur_Pente = Con((surf_Drainee > 12)  &  (surf_Drainee <= 80)  &  (pente_Rclass <= 24),pente_Rclass + 1, Con((surf_Drainee > 80)  &  (pente_Rclass < 24),pente_Rclass + 2,Con((surf_Drainee > 80)  &  (pente_Rclass == 24),pente_Rclass + 1,pente_Rclass)))

The error I get from this is: "TypeError: cannot concatenate 'str' and 'int' objects"

This seems normal, as I have strings (surf_Drainee, pente_Rclass) so I convert the values to str using the following:

facteur_Pente = Con(("surf_Drainee" > str(12))  &  ("surf_Drainee" <= str(80))  &  ("pente_Rclass" <= str(24)),"pente_Rclass" + str(1), Con(("surf_Drainee" > str(80))  &  ("pente_Rclass" < str(24)),"pente_Rclass" + str(2),Con(("surf_Drainee" > str(80))  &  ("pente_Rclass" == str(24)),"pente_Rclass" + str(1),"pente_Rclass")))

The problem is that it doesn't work, and produces the following error: "RuntimeError: ERROR 000732: Input Raster: Dataset pente_Rclass1 does not exist or is not supported."

I don't even have a pente_Rclass1 in my script!

Does anybody know how I could fix this problem in order to continue my script?


texts: Biweekly time series for raster stack in r


texts:
I have got a set of layers, for example, villages, districts and towns.
Their borders are the same on edges of town layer and need to be smooth. I used to perform it manually using topology edit tool. But is there something like advanced topolohy edit tool that can smooth selected lines in several layers.

I tried feature polygon to line and then smooth lines. But it is complicated to convert it back to polygons with original attributes.


texts: How to link to external files with relative path in QGIS actions?

texts:
Objective: I have point layers with attributes contain paths to external files (visio, autocad dwg...). I want to put them all in one folder together with the project file, compress the folder into a zip file and send it to my colleagues for reviews. When my colleague de-compress the zip file, they would be able to open those files with simple actions like the posts linked below.

My observation: what I want to accomplish is very similar to those discussed in this post and this post with one difference: when the zip file is de-compressed, absolute path to referenced files will be different. But I think that RELATIVE PATHS to the project file will be the same, so why not use this? The problem is: I can not find a way to open external files with relative paths inside QGIS.

Any idea how to solve this problem?

Thanks a lot!


texts: Changing Raster Resolution to match other Raster in QGIS

texts:
I have a raster layer of 0.125 degree resolution and another raster layer of 0.0833 degree resolution. Both the raster layers are to be overlapped. So, I want to change the resolution of 0.0833 resolution raster to 0.125 resolution such that the pixels exactly overlap over each other i.e. corresponding pixels signify the same point on the map.
P.S.- The 0.125 resolution map is exactly as it should be and is not to be changed.

I tried Gdalwarp simply to change the resolution to 0.125 but the pixels are not overlapping. Any help!!


texts: Arcpy Intersect analysis function returns an error when using certain types of strings for input

texts:
In an arcpy for the intersect_analysis function, certain syntax works and others don't. What is the reason for this?

Works:

arcpy.Intersect_analysis("'D:\\Program Files (x86)\\ArcGIS\\data\\uscities.shp' #;'D:\\Program Files (x86)\\ArcGIS\\data\\riverbuffers2.shp' #", Cities_Near_Rivers, "ALL", "", "INPUT")
arcpy.Intersect_analysis("'D:\\Program Files (x86)\\ArcGIS\\data\\uscities.shp' #;C:\\Users\\t.thomas\\Documents\\ArcGIS\\Default.gdb\\usrivers_Buffer #", Cities_Near_Rivers, "ALL", "", "INPUT")

Does not work:

arcpy.Intersect_analysis("'D:\\Program Files (x86)\\ArcGIS\\data\\uscities.shp' #;D:\\Program Files (x86)\\ArcGIS\\data\\riverbuffers2.shp #", Cities_Near_Rivers, "ALL", "", "INPUT")

It seems that whenever a 'Program Files (x86)' is there in the path, the tool requires a single quote around the input. What is the reason for this?

This is the error I get when running this as a tool script.

Executing: testscript "D:\Program Files (x86)\ArcGIS\data\uscities.shp" "D:\Program Files (x86)\ArcGIS\data\usrivers.shp" "D:\Program Files (x86)\ArcGIS\data\testoutput23.shp"
Start Time: Mon Oct 05 13:43:07 2015
Running script testscript...
Failed script testscript...

Traceback (most recent call last):
  File "D:\Program Files (x86)\ArcGIS\data\oriscript.py", line 43, in <module>
    arcpy.Intersect_analysis("'D:\\Program Files (x86)\\ArcGIS\\data\\uscities.shp' #;D:\\Program Files (x86)\\ArcGIS\\data\\riverbuffers2.shp #", Cities_Near_Rivers, "ALL", "", "INPUT")
  File "d:\\program files (x86)\\arcgis\\desktop10.3\\arcpy\\arcpy\\analysis.py", line 289, in Intersect
    raise e
RuntimeError: Object: Error in executing tool

Failed to execute (testscript).
Failed at Mon Oct 05 13:43:12 2015 (Elapsed Time: 5.40 seconds)


texts:
I had a dwg file without spatial reference. I converted it to shapefile and gave it a projection: WGS84, but when a try to open the basemap doesn't fit with the shapefile.

I think there's a problem with the projection, but I don't know how to fix it.


texts: Can't send the string output (from a Python tool) to work as a input for a tool in a model (ArcGIS Pro)

texts:
I have a model that uses the Feature Class to Feature Class tool. It needs 3 inputs (Input layer, location and output name). My goal is that a user would choose a layer, the model would create a new name based on the layer, and then it would export the layer to somewhere else. For testing purposes, the input layer is static, because I have no issues with that part.  Here is my model so far:

The change name is a custom tool that take a layer name, and return a string.  Again, for testing purposes, I deleted almost all the code, and left only the important part:

import datetime
import arcpy
import os
name = arcpy.GetParameterAsText(0)
text = "lol"
arcpy.SetParameterAsText(1, text)

Of course, on the tool properties, I made sure that the data type is a string for the output. I also tried Any type, and also tried changing the type to from Derived to Optional and Required.

But I always get the same error. Output Feature Class is not valid on the Feature to Feature tool:

and

I'm not sure what I'm doing wrong.


texts:
When producing a hill-shade from Raster_>Analysis->Hill-shade in QGIS, I cannot stop my hill-shade from becoming too dark, leaving massive parts cities to be black.

I have tried the extremes and in between of light altitude and azimuth, as well as change the z factor to 11120, 111120, and 1111420. These are all suggestions I found on the internet.

I am not sure if it is the fact that I am doing to large of an area? I would like to get the shade as it the light was coming from directly above so that I made see all features.

EDIT: I have added a few screen shots below. The hillshades are followed by their settings. I included the stock settings, and then the best one I have come up with. The problem with the latter is still the amount of black when zoomed in.


texts: Are there any options/plugins to allow labels in busy screens to be visible in QGIS?

texts:
I have quite a busy screen, showing a large number of polygons in a small area. I would like each of these to be labelled neatly, without the 'crowded' feel of overlapping labels, as suggested in the article "When labelling polygons, not all labels are displayed - QGIS 1.7.0".

Instead, I would like to see something as shown on the right (quickly sketched in Paint).

Can anyone recommend a solution?


texts: Projecting OpenStreetMap raster base in other projections

texts:
I know that OSM WMS service works with EPSG:4326 and EPSG:900913
I need to have an OSM base layer projected in another Projection, for example EPSG:32632
Do you know some way to do this?
Even Reproject WMS service, or download new tiles for Geoserver.


texts: Calculating all distances from GPS route using date time and coordinates for php

texts:
I have big GPS data which have 4 columns for id, date time, lat, long. I want to implement Haversine's formula to get distance between 2 nodes.

If I want to find all distance from start to last using GPS data, then how do I set the parameters for $lat1, $lat2, $lon1, $lon2 in my code?

The final result is I can show all routes with the details of maximum and minimum's average speed on map. This is my code :

<?php
$con=mysqli_connect('localhost','root','');
$db=mysqli_select_db($con,'taxi');
$rs=@mysqli_query($con,"SELECT * FROM taxi WHERE id=1 ORDER BY date ASC");
//the id is from ID each Taxi which have 100 ID and want to show each day routing
$lat="";
$i=1;
while ($row = @mysqli_fetch_array($rs, MYSQLI_ASSOC)) {
if($i%2!=0){
//save the start node for latitude and longitude
$awal_lt=$row['lt'];
$awal_lg=$row['lg'];
}else{
 echo getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2);
}
$lat .= "{lat: ".$row['lt'].", lng: ".$row['lg']."},";
$i++;
}
$lat=substr($lat,0,(strlen($lat)-1));

 //haversine formula
function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2) {
var R = 6371; // Radius of the earth in km
var dLat = deg2rad(lat2-lat1);  // deg2rad below
var dLon = deg2rad(lon2-lon1);
var a =
Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2)) *
Math.sin(dLon/2) * Math.sin(dLon/2)
;
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c; // Distance in km
return d;
 }

function deg2rad(deg) {
return deg * (Math.PI/180)
}
?>

//this is my google map route
<!DOCTYPE html>

    <meta name="viewport" content="initial-scale=1.0, user-scalable=no">
    <meta charset="utf-8">
    <title>Simple Polylines</title>

This is the result without average speed :