texts:
I have added the GetFeatureInfo Tool to my Viewer but the PopUp does not have a grid in it to display the Feature Info neatly.

How do I get the grid layout in the Feature Info PopUp, like the one that we are able to see in the popup box in the feature edit tool?


texts: PostGIS ST_Difference function returns wrong geometry in result

texts:
I have a strange st_difference behaviour. I’ve got two overlapping lines (each one has two points), but st_difference result is complete geometry of the first line.

The code to reproduce is here:

WITH
mylines AS
(
SELECT
 st_geomfromtext('LINESTRING(4377625.44 7399187.54,4377633.47 7398278.09)') as ro_geom,
 st_geomfromtext('LINESTRING(4377625.44 7399187.54,4377630.67957507 7398594.12386704)') as lk_geom
)

select
  1 as gid,
  st_difference(ro_geom, lk_geom) as diff_geom
from mylines

Visually, it looks like this:

So the st_difference between red and blue lines returns a green line, which isn't expected. I'd expect the green line to be shorter and below the blue one (without overlapping it).

This is just a test example I've made for demonstration of the problem. In practice, I'm using st_difference for two huge datasets and every difference looks unexpected.

Any clues?


texts:
I would like to create a map, that shows a vehicle traveled path. The map should shows arrowhead for each turning.

I need to show the map like this.

I can create the map by following code

 var flightPlanCoordinates = [
        new google.maps.LatLng(10.51178,76.197739),
        new google.maps.LatLng(10.512117,76.261253),
        new google.maps.LatLng(10.554309,76.322365),
        new google.maps.LatLng(10.56646,76.372833)
    ];
    var flightPath = new google.maps.Polyline({
      path: flightPlanCoordinates,
      strokeColor: '#04770C',
      strokeOpacity: 1.0,
      strokeWeight: 2
    });

    flightPath.setMap(map);

I am using Google Maps API V3.


texts: Arcpy options to copy feature class to new feature class (fastest runtime)

texts:
I have a simple script that copies four feature classes from a sde geodatabase and places them in a file geodatabase.  Three of the feature classes run quickly but the forth feature class has a few hundred-thousand records.  The larger feature class takes about five minutes to process.  I'm currently utilizing the feature class to feature class tool.  I've also tested with the copy tool and experienced similar runtimes.  I noticed in ArcMap the fastest way to copy the data is to right-click on the layer and select Data-->Export Data.

Is there a specific tool in the "system tools" that works off the same algorithm or is there a tool that is quicker than the "feature class to feature class" and "copy" tools and also functionally equivalent?


texts: recalculate latitude and longitude from a pixel in a cropped map

texts:
i want to recalculate a geographic coordinate from a pixel in a map.
For this i have the size and the coordinates from the map

to convert a coordinate to a pixel is use the algorithm, you find at https://gist.github.com/gagern/6636176

But how can i convert this back to the orginal geographic coordinate?

for the longitude i simply use

 Math.toDegrees((x+from.getLongitudeAsRad()*xFactor)/xFactor)

but how can i calculate the latitude?

best regards


texts: What's the keyboard shortcut for select/deselect (on/off) multiple layers in QGIS?

texts: As in ArcGIS, CTRL+CLICK, ALT+CLICK, SHIFT+CLICK works for multiple selection of layers, I didn't find a similar one in QIS. Is there any?


texts: GDAL (Python) open/load raster with offsets

texts:
I am using the Python GDAL API to perform operations on geotiffs.  I have a good understanding of what needs to be done, but the how and particularly interpreting GDAL's documentation is a barrier.  Basically I have two geotiff arrays with different extents but they overlap ~ 85%.  I just want to load each image as a numpy array in Python clipped to their common extent.
I have used this great code from another post (How to add different sized rasters in GDAL so the result is only in the intersected region) but I am hung up on the arguments for band.ReadRaster() at the end of the accepted answer.
I can find the offset number in pixels (hypothetical example) - xoff1 = 20, yoff1 = 100, xoff2 = 1800, yoff2 = 2000 if hypothetical tiff1 size = 2000 x 2250. So for an equivalent numpy array I would clip it like this: np_array[20:1800,100:200]
With two datasets loaded as ds<1 and 2> = gdal.Open(filepath/to/tiff) and all the metadata, methods, objects, etc. HOW do I finally use my extents to open clipped bands?  I want to use band.ReadRaster() I just don't understand the arguments.
Here is example from API tutorial of using the method: scanline = band.ReadRaster(xoff=0, yoff=0, xsize=band.XSize, ysize=1,buf_xsize=band.XSize,buf_ysize=1, buf_type=gdal.GDT_Float32.  And I have referenced the gdal documentation which looks the same as this example.  I have TWO offsets per northing/easting why only one value for offset?  I am confused.


texts: What is the correct way to get road segments from OSM based on a route?

texts:
I'm using a Postgresql+postgis database with OSM data (just the Netherlands). I would like to get a list of all road segments based on a set of LAT/LONG coordinates.

For my test I use (part of) a random GPX track which you can easily visualize on for example https://clydedacruz.github.io/openstreetmap-wkt-playground/:

LINESTRING(4.62925 52.66297,4.62873 52.661519999999996,4.62967 52.66099,4.62994 52.66038,4.63432 52.66063,4.63506 52.66077000000004,4.63619 52.661460000000005)

Below is my query:

select st_astext(st_transform(way, 4326)), highway from planet_osm_line straten where st_dwithin(straten.way, st_transform(ST_GeomFromText('LINESTRING(4.62925 52.66297,4.62873 52.661519999999996,4.62967 52.66099,4.62994 52.66038,4.63432 52.66063,4.63506 52.66077000000004,4.63619 52.661460000000005)', 4326), 3857), 0.1) and not straten.highway is null

This query returns 12 road segments of which one is the Torensduin street which isn't part of the actual route but is part of the results just because the linestring between two point crosses this street.

Road segment from the query result for the Torensduin street:

LINESTRING(4.62970101935771 52.660935187484,4.62952629703495 52.6608663721515,4.62933181177594 52.6606963765665)

Also this (nameless) street  is added to the road segments but isn't part of the route:

LINESTRING(4.62970101935771 52.660935187484,4.62952629703495 52.6608663721515,4.62933181177594 52.660696376566506)

What I have tried so far is:

Use MULTIPOINT instead of LINESTRING, but that query only returns a
very limited list of road segments because not every road segment on
the route has a point in the GPX track
Use ST_INTERSECTS instead of ST_DWITHIN, but (obviously) with the same result.

I think I'm using the wrong approach with my query but haven't got a clue on what approach I should use to get the actual road segments.

So the question is, how can I get a 99,99% correct list of road segments based on a list of LAT/LON coordinates?


texts: Choosing geometry to store tin format data to support spatial queries?

texts: What is the best way to store 3D TIN data in postgis, to support queries like st_intersects, st_within, st_overlaps, etc. I can't use TIN since it only supports st_intersects.


texts: How to replace every element in a matrix or raster with the maximum of the values within a variable moving widow of size = n?

texts:
So essentially, what I want to do is replace every element in a matrix with the maximum of neighboring cells within a window that is determined by the value in that cell.
The window size would be determined by this function (fitlwr), where Tree_Height calls a linear model that was fit to a dataset of Tree Height and Crown Diameter data:
RoundOdd <- function(x) {2*floor(x/2)+1} #makes sure window size is an odd number

fitlwr <- function(x){for(i in x){
  if(i > 13){
    m <- RoundOdd(Tree_Heights[Tree_Heights$Tree_Height == i, "fit.lwr"])
  return(matrix(1, nrow = m, ncol = m))
    }
  else {
    return(matrix(1, 3, 3))
    }
}}

I then want to replace every value in that matrix with the maximum of the values within that window.The matrix was derived from a raster layer and the values represent the height above ground for a given cell. The dimensions are 6,571 x 5,764. A section of the data might look like this:
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    9   47  103   58   80   55   72   56   14    52
 [2,]   68   49   49   43   62   80   62   23   55    82
 [3,]   58   10   79   70   75   49   68   60   74    79
 [4,]   78   19   51   26   61   77   57   70   51    43
 [5,]   47   88   57   80   25   33   24   30   56    63
 [6,]   73   36   53   25   63   30   19   59   17    63
 [7,]   95    9   49   95    6   13   21   75   60    34
 [8,]   36   65   47   64   22   66   52    9   71    20
 [9,]   45   53   31   47  114   55   44   42   44    44
[10,]   47   23  102   34   67   60    5   23   61    32

The raster focal functions were my go-to, but they don't let you use a variable window size (see below).
RoundOdd <- function(x) {2*floor(x/2)+1}

fitlwr <- function(x){
  RoundOdd(Tree_Heights[Tree_Heights$Tree_Height == x, "fit.lwr"]/2)
}

m <- raster::focalWeight(x = CMM, d = fitlwr(), type = "circle")

CMM <- raster::focal(x = CMM, w = m, fun = max)

This returns the following error:
Error in `[.data.frame`(Tree_Heights, Tree_Heights$Tree_Height == x, "fit.lwr") : argument "x" is missing, with no default
6.`[.data.frame`(Tree_Heights, Tree_Heights$Tree_Height == x, "fit.lwr")
5.Tree_Heights[Tree_Heights$Tree_Height == x, "fit.lwr"]
4.RoundOdd(Tree_Heights[Tree_Heights$Tree_Height == x, "fit.lwr"]/2)
3.fitlwr()
2..circular.weight(x, d[1])
1.raster::focalWeight(x = CMM, d = fitlwr(), type = "circle")

If I try instead to use the function in the argument for window size, I get this error:
Error in .local(x, ...) : is.matrix(w) is not TRUE
5. stop(simpleError(msg, call = if (p <- sys.parent(1L)) sys.call(p)))
4. stopifnot(is.matrix(w))
3. .local(x, ...)
2. raster::focal(x = CMM, w = fitlwr, fun = max)
1. raster::focal(x = CMM, w = fitlwr, fun = max)

I am open to using another language or software tools to accomplish this task, including GRASS, Python, QGIS, or ArcGIS if necessary.


texts: Interactive bearing measurement in Esri Production Mapping?

texts: I have ArcGIS Desktop Advanced 10.4 and I want to use ARC Production Feature Builder to measure the bearing between two points just entered on the map, whilst there is a description of this function in the ESRI info, I believe it is only in the ARC Maritime extension. Am I mistaken with this assumption? if not how do I access it , it is not in the Customise>toolbars list. Otherwise is there a similar function available as an add in ?


texts: How to add rows from a Dataset to another?

texts:
The DB I have is divided in many datasets (shape format, each .zip I send to CartodDB consist of 4 files, .dbf, .prj, .shp, .shx). The problem is that the free account only allows 4 datasets, each one associated to a layer.

To solve this problem I'd like to join all my datasets in just one. All datasets share the same structure (same columns) due to the fact they form part of the same DB. I hope there is a way to add rows from one dataset to another.


texts: How to perform spatial interpolation using PostGIS?

texts:
I need to interpolate some point data and would like to use the ST_InvDistWeight4ma function.

As I am new to raster PostGIS, I would like to get some help about this.


texts: Creating a polygon outside another but using existing as boundary

texts:
I need two polygons. I already have one (the river edge boundary as seen below)

What i want to do is create another polygon which covers the part of the raster which does not fall inside the river polygon. It will therefore need to snap to the outside of the river polygon which is a complex shape. I have tried the auto-complete polygon tool but not having much luck.
Furthermore, i have part of the river polygon which is separate from the rest. This would need to be classed as the river polygon and not outwith.


texts:
I have some recent data that is at the Lower Super Output Area (LSOA) level for Greater Manchester. I wrongly assumed that the data was for the latest LSOA boundaries, i.e. for the 2011 boundaries.

Acting on this assumption, I downloaded the 2011 boundaries from ONS and used mapshaper to simplify it to 1% of discardable vertices. I imported the data in to CartoDB and got this nice looking map:

https://gmpagdata.cartodb.com/tables/gm_2011_lsoa_boundaries_1pc/public/map

Once I realised that my data did not fit this LSOA version, I tried to do the same with the 2001 ONS data. This was less successful. Even retaining 2.5% of discardable vertices still leaves a lot of holes in the map, as you can see here:

https://gmpagdata.cartodb.com/tables/gm_2001_lsoa_boundaries_2_5pc/public/map

I used the following command to do the simplification and conversion before using a simple ruby script to select just the Greater Manchester features:

mapshaper -i LSOA_2001_EW_BFC_V2.shp encoding=utf8 -simplify 1% keep-shapes -o uk-2001-lsoa-boundaries-1pc.geojson format=geojson

Other notes:

The encoding parameter in the command above is a guess. Both sets of data required an explicit encoding and utf8 worked fine on 2011 data.
I think I used the same command for the 2011 data (there are not a lot of parameters to play with) but I may be wrong.
The results that I get with the full, clipped 2001 data look like the results I got when I unknowingly used the 20m generalised 2011 data. In fact, the first thing I did when I saw the results was to check that I had not used 20m generalised data.


texts: Extracting spectral values from distinct points in ENVI Classic

texts:
I have 354 distinct points (as a shapefile) that are very close together i.e. 2-3 points have the same lat but different long or vice versa as the sampling was done in a very small area and the field conditions were difficult.

I have imported this points shapefile as a new vector and have checked the attributes to see if the points have distinct lat & long (together). I have then exported the vector layer into roi and then exported the roi into an ASCII sampling the RGB bands.

Ultimately, I get 354 points having duplicate lat long and there are only 37 distinct points, instead of 354. Are the points too close to be treating each point as a distinct roi?

If yes, then how do I extract the spectral values at these points. (Exporting to Arcmap and sampling is not an option)


texts: How can I extract the measurelength, highmeasure and lowmeasure of routes in ArcGIS

texts:
I'm trying to create road Routes in ArcGIS, which traditionally have been created using workstation. I've managed to write a Python script to do this, however, unlike ArcInfo, the attributes MEASURELENGTH, LOWMEASURE and HIGHMEASURE are not automatically created as you could do with the WorkStation 'routestats' command.

Anyone know how to do this in ArcGIS or Python?


texts: Can Hazus tables be downloaded into table or Excel format?

texts: Can a Hazus user download results in table format to be sent to a user who does not have Hazus on their ArcMap?


texts: Connect Qt DoubleSpinBox to PyQGIS script

texts:
I have done a python script and I want to change the '1.0' value on pas15 with the value that is set in the double spin box, but I can't find the way to do it.
Any idea?

self.dlg.doubleSpinBox.valueChanged[str].connect(self.value)

def value(self):
    self.value=QDoubleSpinBox.getValue()
pas15=processing.runalg('saga:contourlinesfromgrid',self.input,0.0,5000.0,1.0,None)


texts: Feeding direction of the road segments in QNEAT3 plugin

texts: I am going to use the Iso-Area algorithm of the QNEAT3 plugin to generate Isochrones. I need to determine which parts of the road network are one-way and which parts are two-way. In the way that for example assign value 1 to the records in the attribute table of road layer that are one-way and assign the value 2 to the records that are two-way then feed the direction filed, value for forward direction, and value for backward direction with the related field and values. Is this a correct way to generate Iso_Areas based on the one-way and two-way road segments?


texts: How to create cost maps based on OSM (landuse/...)?

texts:
I'd like to create cost maps (georef raster) for cross-country hiking based on OSM from pgsql/postgis plus height data.
Maybe by creating multiple raster-files and the later merging them (lake=unpassable, brige/ferry will later overwrite with a path).

Can someone please guide me to helpful threads or a "howto" for creating georef raster from psql/postgis OSM data without use of a tileserver.

Would gdal_rasterize be the right approach, or is there a better/faster way ?


texts: Convert meta tiles names to an x,y,z notation in OpenStreetMap

texts:
Using mod_tile/renderd to serve OpenStreetMap tiles, I would like to list all tiles already generated in the cache but with an z, x, y notation. Currently, I only get paths for my list of *.meta files

cd /var/lib/mod_tile/ && find `pwd` -name '*.meta' |sed -e 's/\/var\/lib\/mod_tile\/default\///g' | sed -e 's/\.meta//g'

The remaining problem is to convert meta tiles paths to the z, x, y notation.

I already tried to:

find an existing Python function using a search engine
to reverse the logic of function xyz_to_meta that do the opposite  (in the Github mod_tile repository) but I'm really bad at bit manipulation (e.g the definition of the meta tiles pattern )
convert this available C code to Python

Python or Bash answers are preferred because I don't want to add more packages on the server.


texts: Keep Column Order when renaming Attribute

texts:
Im renaming an attribute by

layer.startEditing()
    layer.renameAttribute(layer.fieldNameIndex("Old Name"), "New name")
layer.commitChanges()

Which works fine. But in the attribute table, the renamed column is then displayed to the very right, after all other columns.

How do I keep it at its position?

Its just the display order and not re-ordered in the file behind, but still unhandy.

I've seen this on other changes to Attributes, too. However, the only one I can reproduce is this one.


texts: Standard for publication-ready R-based geographic images?

texts:
I've been using the ggmap package to create figures that will go in journal articles for many years. These maps are a little harder to create than Leaflet maps, and as OpenStreetMap and Google have both shut down their tile APIs, getting background tiles that work for me has become more difficult. (Stamen is okay, but the selections don't always work).

I've seen ways to export leaflet maps as static images, but these often leave artifacts like zoom buttons that I definitely don't want in my paper. What is the best tool now (2019) for creating clean, publication-quality geographic images in R?

As a further complication, I usually write my articles in rmarkdown so that my figures can update with changes to the data, research, etc. Taking a screenshot of a leaflet map where I just crop out the things I don't want isn't ideal.


texts:
I did a Network Analysis with ArcMap, executed by a python-script (I am building an ArcGIS-Tool).

In the very end, I want my directions-output to be shown as a displayed map, so the user can navigate with the map.

But right now I only get an HTML/XML/Text-File, which is nice, but not perfect.

# Process: Directions
arcpy.Directions_na(Route1__6_, "HTML", Directions2102_html, "Kilometers", "REPORT_TIME", "Minutes", "en", "NA Desktop", "")

Is there a way to display a map like ArcMap does, when I run the Network Analyst tool Find Routes?


texts: Data smoothing of Sentinel-2 NDVI time series chart in Google Earth Engine?

texts:
I am exploring Sentinel-2 time-series NDVI with Google Earth Engine. In another post I learned how to plot a NDVI time series chart.

I know that there are many smoothing methods for reducing NDVI noise, as per this article. I would like to apply one of them directly in GEE (not R), for example the Whittaker filter.

Can I use any of these methods in GEE? How?


texts: Map of Russia's administrative divisions

texts:
I'm looking for a map of the administrative divisions of Russia, on a rather local level. In Seeking polygon files corresponding to Russian Census? is presented a map on different levels, the most local being just one level under the oblast level I guess :

I would like to know if it possible to find a more detailed administrative map, on the municipality level if possible.


texts:
Just trying to figure out how to delete pages from a layout in QGis 3.0

I have found where to add pages. Maybe there is no function yet for deleting pages?


texts: It is possible in MongoDB do a query order by the near distance to a given GeoJSON point?

texts:
I have in MongoDB a collection with GeoJSON Points that represent places of interest like this:

{
    "name" : "Malba",
    "type" : "Museum",
    "location" : {
        "type" : "Point",
        "coordinates" : [
            -58.474184,
            -34.550208
        ]
    }

}

For example, I could have many points of type "Museum" with different location. It is possible to make a query in MongoDB to recover all this point order by the near distance of a given location?


texts: Compiling on CentOS fails due to ignored LIB flag

texts:
I need help with installing PostGIS.
In general: I try to install it on a cluster (CentOS) and do not have admin rights. There is a sqlite3 installation (connected to the error), but this seems corrupt. My own one usually works and has (also by install-process) been used previously.
I cannot use yum.
In detail:
For having the MakeFiles build I set flags as I believe to possibly need them:
./configure --prefix="$HOME/.local" CFLAGS="-I$HOME/.local/include" LDFLAGS="-L$HOME/.local/lib" LT_SYS_LIBRARY_PATH="$HOME/.local/lib" LIBS="-L$HOME/.local/lib" SQLITE3_CFLAGS="-I$HOME/.local/include" SQLITE3_LIBS="-L$HOME/.local/lib"
But the following error is thrown when running make:
/usr/bin/sed: can't read /usr/local/lib/libsqlite3.la: No such file or directory
libtool:   error: '/usr/local/lib/libsqlite3.la' is not a valid libtool archive
Right after (but I don't see the connection):
libtool: link: ranlib .libs/liblwgeom.a
If I run make again, it jumps that directory. But gets stuck with a similar error (which is not jumped by running make another time):
/bin/sh ../../libtool --mode=link gcc -Wall -Wmissing-prototypes -std=gnu99 -I<exp. $HOME>/.local/include -fno-math-errno -fno-signed-zeros  -fPIC -DPIC -I../rt_core -I../../liblwgeom -I<exp. $HOME>/.local/include   -I/usr/include/gdal -I<exp. $HOME>/.local/include   ../rt_core/librtcore.a raster2pgsql.o -L<exp. $HOME>/.local/lib -static ../../liblwgeom/liblwgeom.la -L/usr/lib64 -lgdal  -L<exp. $HOME>/.local/lib -lgeos_c  -lc -lm -o raster2pgsql
libtool:   error: cannot find the library '/usr/local/lib/libsqlite3.la' or unhandled argument '/usr/local/lib/libsqlite3.la'
I don't get why libtool searches in '/usr/local/lib/'. I didn't find that hardcoded in libtool or the MakeFiles. I'm clueless - and this is above my skill (apperently).
Can you help me?


texts: Iterator and List of Values in ModelBuilder?

texts:
My model iterate through many shapefiles:

First it checks if, the field I want to add (model parameter), already exists in the table (check if the field name is the same).

If it doesn't exist:

Great! Just add the field I want.

If exists:

Creates a temp field
Store the the value from the existing field in the temp field
Drop the existing field
Create the field I want, the way I want (string, with exactly 137 length)
Copy the values from the temp field to the new field
Drop the temp field

It works perfectly, always. As many shapefiles I want.

The problem is I want to run it several times, checking for a lot of fields, in a single run in the tool.

If I try to bath it, doesn't work. And if I try to put 'List of Values' in the first parameter 'Multiple Value' it gets error, because model builder doesn't allow both 'model iterator' and 'list of values' working together.

Can you guys help me?

download model - wetransfer


texts: Looking for pedestrian-specific routing engines

texts:
Most of the applications providing pedestrian routing I see are sufficient for many needs but not optimal: you cannot cross open spaces like town squares, routing is in the middle of the street etc.

Do you have examples of websites or applications with routing really specific to pedestrians?


texts:
I have a mosaic dataset (Raster) in an SDE database. As expected, when I click start editing, ArcMap says that I have to do "Register as Versioned" before, but I can't have a mosaic dataset in a Feature Dataset. Can I?? Can I register this mosaic dataset as registered in another way?

Normally, when a mosaic dataset is in an File geodatabase, you can start an edit session and edits the Boundary and the Footprint. I just want to do the same operation but with my mosaic dataset located in an ArcSDE server.

Is this possible? If yes, how?

Thanks!

Please note: If it can be useful, I use direct connections to my SDE database.


texts: A patch of water that is part of a shapefile is painted as land

texts:
So I've set up OpenMap in a Java application, but the default world map shapefile (dcwpo-browse.shp) that comes with OpenMap wasn't sufficient enough, since the area I specifically needed wasn't fully drawn in this default file.

Which is why I ended up using a shapefile that I found somewhere on the internet. But OpenMap seems to want a .ssx file accompanying the .shp file, but this other shapefile doesn't. And many other shapefiles I still have to use don't have a .ssx file either, and some of them don't even display on the OpenMap MapBean at all.

And even though this new map shapefile works, it seems like it just fills every shape, even the one that is actually supposed to stay transparent because it's water. This is probably because it's an enclosed area due to a bridge that is drawn, which makes that particular area appear like another piece of land. And I can't help but wonder whether it might be because I lack the .ssx file, or because of something else.

What's the best way to fix this problem?


texts: Get distance between multipatch features in ArcGIS

texts: I've created 2 multipatch feature classes in ArcGIS, which represent the outer wall of two pipe systems. Now I want to check the smallest 3D distance between all features in the 2 feature classes. I was trying to use Near 3D, but that results in a NEAR_DIST of zero.


texts: How to move labels in QGIS using the old labelling method

texts:
I would like to use the new labelling method in QGIS, however it does not display labels for all of my polygons. If I use the old labelling method (in the properties dialog) the labels display fine.

Is there an easy way to move labels using the old labelling method other than to specify a label coordinate in the table? I would like to be able to click on the label to move it like in the new label tool?

This question asked about the problem I am having with the new labelling tool, but in this question I am interested in how to move labels using the old method.


texts:
I have a challenge with advanced labels of polygons in ArcMap 10.1. I wish to make the label take the longest direction(either curved, horizontal, vertical or diagonal) inside the polygon.
The python code I used is below. I need a modified python code to specify and orient the labels along the longest direction of an ArcMap polygon.
def FindLabel([BLOCK], [CLONE], [PLANTED_DATE], [Area]):

    labels1=' '
    if [CLONE] is not None:
        labels1 =("<CLR red='150'><FNT size = '6.7'>" + [CLONE] + "</FNT></CLR>")

    labels2=' '
    if [PLANTED_DATE] is not None:
        labels2 = ("<CLR blue='255'><FNT size = '6'>" + [PLANTED_DATE] + "</FNT></CLR>")

    labels3=' '
    if [Area] is not None and [Area] !=0:
        labels3 = ("<CLR red='255'><FNT size = '5.5'>" + [Area] + "</FNT></CLR>")

    lab =  labels1 + " " + labels2 + " " + labels3
    return lab


texts: Using multiple fields as parameters in QGIS Processing Algorithm Script

texts:
I have a script looks like this:

##Input_Layer=vector
##Date_Field=field Input_Vector

# Do Something...

So normally it would like one vector layer and one field, and the choice of field would be depend on the Input_Layer. However I want to make it to take multiple fields.

I should be looking for something like "multiple raster" or "multiple vector" as shown below. It allows user to input as many as they want. But there is no "multiple field". Is there a way to do that?


texts: Can I use SELECTION mode on FeatureLayer from ArcGIS.com Webmap in ESRI JS API

texts:
I am using an ArcGIS.com Webmap in an ESRI Javascript API application I am developing.  One of the layers in my Webmap is rather large and cannot be rendered all at once when zoomed out, however I would like to avoid showing holes in the data on screen when data does actually exist there.

I have seen samples and demos before where it is possible to use a tiled map service to cover the entire large area without rendering all the individual features, but still be able to show an InfoWindow when the map is clicked.  I am assuming this would be done using a FeatureLayer with its mode set to SELECTION and the onClick handler for the map would query the FeatureLayer based on where the map was clicked, then the selected features would appear on the map and the InfoWindow could be visible.

I have not seen where it is possible to set the mode of a FeatureLayer that is provided via an ArcGIS.com Webmap.  If I simply create a new FeatureLayer in my application that points to the Hosted Feature Service on ArcGIS.com, then I lose all of the configuration (symbology, popup, etc.) that I have done when publishing the Webmap.

Is what I have described above possible?  If not, is there another way to render a large number of features on the map and still have access to attribute information?


texts:
I love the new Display Filters feature in ArcGIS Pro. I've used them a few times and they're great! However, recently when I've tried to utilize them they don't seem to work anymore. Below is a screenshot of the settings I have and the data. The data is mileposts on a pipeline. I want to show only every tenth milepost on the first filter at a scale range of 1:50,000-1:5,000 (SQL query is (MRL_MEAS%1<>0)), every whole milepost on the second filter at a scale range of 1:100,000-1:50,000 (SQL query is (MRL_MEAS%1=0)), and only every milepost divisible by 10 on the third filter at a scale range of 1:20,000,000-1:100,000 (SQL query is (MRL_MEAS%10=0)). In the screenshot below I'm at a scale of 1:500,000, so only mileposts divisible by 10 should be shown, but ALL mileposts are being shown. I'm not sure what's going on, because I haven't had this issue before, but maybe I'm doing something wrong.

Any help?


texts: Geoserver + Boundless - web app icons being clipped

texts:
I have created a map app, that displays a few points with their attribute data.
The problem I am facing is that the point icons get clipped  when serving the map.
I have read this answer as well as this one, but none of them seem to be working for me.

This is what I get when debugging and deploying my app.

As you can see, the icons get clipped and I cannot seem to "fix" this.


texts:
I'm in the process of attempting the set up of a GeoNetwork instance and have been given an xls with metadata for 1200+ maps.

Has anyone seen docs for a good workflow for batch import? Maybe xls templates which match the GeoNetwork ISO19139 templates? Or a postgres DB template I could copy them into which I could then connect to GeoNetwork?

I tried to export the GeoNetwork ISO19139 templates as xml and bring them into excel but the formatting gets lost and their not usable.

Ideas?


texts: Replacing geometry of a polygon using the geometry of a polygon from another layer

texts:
I am trying to update a polygon shapefile by redoing the geometry on certain polygons to match the geometry of polygons in a separate, more accurate shapefile. Some of these polygons are extremely complex, and the replace geometry tool is taking much too long.

Is there a tool or method that would allow me to alter a polygon to match exactly the geometry of a polygon in a separate layer while retaining the attributes of the original polygon?

I vaguely recall doing this in ArcMap, but ArcGIS Pro is throwing me for a loop.


texts: Counting number of layers selected in ArcMap TOC using ArcObjects?

texts:
Has anyone ever tried to get the number of layers that are selected in ArcMap?  If I have a map document that has six layers and three of them are selected in the TOC, how does one count them using ArcObjects?  I have looked for a while at IContentsView.SelectedItem.  Documentation says it can be an enumerator, but I am having a hard time finding out actually what it's type is.  In C#, if I do

string myType = IContentsView.SelectedItem.GetType().ToString();

myType is "System._ComObject".  Not very descriptive...  Trying

if(IContentsView.SelectedItem is IEnumLayer)

and similar checks has been unintuitive... I wasn't sure if multiple layers were selected that they would be stored in there but that was my first guess.


texts: ogr.BuildPolygonFromEdges fails on my multiline

texts:
Using python ogr, polygon = ogr.BuildPolygonFromEdges(multiline) fails for the following multilinestring:

Does anyone know why?
It throws the error: RuntimeError: OGR Error: General Error with exceptions enabled. Without exceptions:
ERROR 6: OGR Error: General Error
Edit: I think I should be using polygon = multiline.Polygonize() but this returns an empty geometry collection...
Edit2: The red line has 2 points, the green has 3. Do these need to be segmentized?
Edit3: Code that generated multiline:
import numpy as np
import netCDF4 as nc4
import matplotlib.pyplot as plt

#ogr.UseExceptions()

ds_shp = None
drv_shp = None

# GoM
wkt_boundary = (
    'LINESTRING (-67.25 44.75 , -66. 44. , -66. 38. )',
    'LINESTRING (-71. 48. , -71. 38.)',
)
shapefile_boundary = (
    'topography/etopo1_GoMSS_0005m_edit.shp',
    'topography/etopo1_GoM_0200m.shp',
)

drv_mem = ogr.GetDriverByName('Memory')
ds_mem = drv_mem.CreateDataSource('source')
drv_shp = ogr.GetDriverByName('ESRI Shapefile')

srs_in = None
multiline = ogr.Geometry(ogr.wkbMultiLineString)
for shapefile_in in shapefile_boundary:
    ds_shp_in = drv_shp.Open(shapefile_in)
    layer_in = ds_shp_in.GetLayer()
    if srs_in is None:
        srs_in = layer_in.GetSpatialRef()
        layer_out = ds_mem.CreateLayer('GoM_lines', srs=srs_in, geom_type=ogr.wkbPolygon)
        layer_defn = layer_out.GetLayerDefn()
    for linestring_feature in layer_in:
        geometry_ref = linestring_feature.GetGeometryRef()
        multiline.AddGeometry(geometry_ref)

for wkt in wkt_boundary:
    linestring = ogr.CreateGeometryFromWkt(wkt, reference=srs_in)
    multiline.AddGeometry(linestring)

With:
print(multiline.GetGeometryName())
MULTILINESTRING


texts:
I am creating a Replica from an ArcSDE database using python/arcpy
(CreateReplica_management (in_data, in_type, out_geodatabase, out_name, {access_type}, {initial_data_sender}, {expand_feature_classes_and_tables}, {reuse_schema}, {get_related_data}, {geometry_features}, archiving))

and the tables which are being replicated are empty, only the schema is replicated. Is it possible to replicate also all the data in the tables using arcpy?


texts: I have downloaded landsat satellite data, then I merged the 4 bands I am interested in into one file (raster-merge). The problem is that I merged bands 2,3,4 and 8, and in the merged file the bands are named 1 to 4. The question is which band in the merged file corresponds to which individual raster band file.


texts: Field data inspection and collection using GeoODK?

texts:
In my workplace we are looking at collecting and inspecting GIS data using free mobile apps. The one that is fit-for-purpose from a data collection point of view is GeoODK. I have been able to prepare custom forms and get them to work in the field really well.

We already have GIS data in a proprietary GIS file format (shp, tab etc.) and we would like to inspect them on a mobile app.

Ideally, I would like to have this loaded onto an app and then have the inspections done by changing a few fields such as 'Inspection date', 'Inspection Officer' etc. and then capture new data along the way if these did not exist in the original GIS file.

Can this be achieved using GeoODK?


texts: Load dojo from local and ESRI packages separately

texts: Is there by any chance a way to load dojo packages from my local and then load ESRI javascript api packages which i need?


texts:
I have four files: .prj, .qpj, .shp and .shx.

They are properly read by Luciad (it seems that only .shp and .prj are necessary for that).

Now I would like to import shapes into QGIS taking into account projection (that doesn't happen). I'm using "Add Vector Layer" option for this, maybe that's the problem?

I searched for it a lot and from what I understand .qpj is some QGIS equivalent of .prj. Mine look quite similar. But I'm still not sure which CRS should I use and how, I guess it should be somehow possible with those files.

My .prj:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

My .qpj:

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

How do I know that projection isn't applied? I can only suppose that - this is example comparision between QGIS layer and Luciad map result (expected):

Maybe I shouldn't mention Luciad here, because it's not that related, but that's how I know that my data are good.

EDIT: I see that .qpj/.prj file is used (if none of them is present, QGIS asks me for a CRS, but it doesn't seem to make any difference when I choose just default WGS 84). However, for some reason something (CRS?) from those files isn't applied in the way I want - like in Luciad.


texts: Using Intersect to find where line crosses itself in ArcGIS Desktop?

texts:
Can I not generate point features where a line feature intersects itself by using the Intersect tool?

I try using it with just one line input but when I run it I get a warning saying the output generated is empty. There are visible intersections in the GPS data so it should produce some sort of output.

I'm only using a basic license so unfortunately I can't use topology rules to find these intersections.


texts: Is it possible to calculate Euclidean distance but by using a road network?

texts:
Using ArcGIS I want to output a raster for the distance from a set of points, however I want to do this using a road network, as opposed to just using the Euclidean distance function i.e. using Euclidean distance but along a road network as opposed to as the crow flies. Can anyone tell me away of doing this?
Ultimately I'm trying to calculate a cost surface of how much it will cost to transport goods from  supply points to an undefined market i.e not a specific point but anywhere within a region.


texts:
Volume of water in basin with reference to incremental height of DAM?

I know there's already a link that answers this question. I am not quite sure about the results from Surface Volume Tool. It seems that the results are way above what I expect.

Supposed the results show:

275, 256, 949

but the Tool provides me I am using ArcGIS 10 and I my raster surface is ASTER DEM.

I would like to do, basically, the same tasks as indicated in the link: To calculate volume from ASTER-DEM with increasing height.


texts: osm2po: import multiple .osm files for disparate regions?

texts:
Does osm2po support importing multiple OpenStreetMap files into the same Postgres/PostGIS database (i.e., to be used for routing)?

If yes, can the OSM files be for geographic regions that are not next to each other? For example, could I import something like NewYork.osm and SanFrancisco.osm into the same routing database?

Thanks!


texts: subset of shapefile changes value when defined "as.numeric" in R

texts:
In R, I have a shapefile with polygons. I subsetted by a secondary ID. Now, I want to sum up all area values of the subset.
These are the are values of the subset:

shapefile.subset$area
## 7767 1597 947  626  552
## 1424 Levels: 0 1 10060 10070 10124 1013 1017 10175 1019 10191 1022 102361
## 10237 10298 1030 10324 10327 1035 1039 10403 10405 1042 1044 104453 104533
## 10508 1051 10562 ... 9982

Summing returns the following:

sum(shapefile.subset$area)
## Error in Summary.factor(c(1247L, 275L, 1388L, 1101L, 1025L), na.rm = FALSE) :
  ‘sum’ not meaningful for factors

So I tried to convert into numeric

as.numeric(shapefile.subset$area)

which returns

1247  275 1388 1101 1025

Thus five values that do not have anything in common with the original values. Note, that these also appear in the error message of sum(shapefile.subset$area). I am really stuck here - how can I avoid this strange behaviour?


texts: When using case when in QGIS 3, even without else, all other previously filled cells are replaced by NULL

texts:
I would like to update a column with different strings based on other strings in another column.
For instance : if A, B, C in column 1, write Cat1 in column 2. So I use "case when". But even without writing 'else' when it always replace all other values by NULL. Which means that if I then want to update again the column with "if X, Y, Z in column 1, write Cat2 in column 2", "Cat1" gets deleted in all the cells where it had been added. If I don't want to update all cells at the same time, how can I conditionally update certain cells while ignoring others (i.e. not deleting their content)?


texts: Capture coordinates of a rectangle drawn on the canvas

texts:
I am writing a process that will take parameters as input, and then process part of a raster layer. One (well, four) of the parameters should specify a rectangle on the canvas. Ideally the user should be able to draw a rectangle - left click and drag - and the ulx, uly and lrx, lry coordinates be passed to the process.
Now, I have found some code which allows the drawing of a rectangle on the canvas but I do not know how to capture the ulx, uly, lrx and lry for input to the process. I am using QGIS 3.16.
Can anyone point me at some existing PyQGIS that does this, or give some pointers as to how to do this?


texts: How to link pgrouting edges with street names to route user based in street names?

texts:
I have installed postgreSql + postGis + pgrouting and I have loaded OSM data into my database using osm2po and shp2pgsql. How can I search for points in the map using the street names and then route the user from point 1 to point 2. pgrouting divides the streets in edges and any given street will have multiple edges and nodes, how could I link those edges and nodes to the actual street names? Any ideas?
Thanks


texts: Can I supply my own REST API token to the ArcGIS Javascrpit API?

texts:
I am using the ArcGIS Javascript API to connect to the arcgis.com portal. To do so I am calling portal.signIn(). The issue is that I don't want the user to be responsible for performing the login, bo I want to generate the REST API token on my server and then use that token from the Javascript API.

Mock code

portal = new esri.arcgis.Portal('http://www.arcgis.com');
dojo.connect(portal, 'onLoad', function() {
    portal.id.token="<server-side writes token here>";
    portal.signIn().then(function(loggedInUser){
        console.log(loggedInUser);
    });
});

EDIT:

From what I can tell in this sample http://help.arcgis.com/en/webapi/javascript/arcgis/help/jssamples_start.htm#jssamples/widget_identitymanager_client_side.html

I might be able to build my own credentials object to run esri.id.initialize() (which from what I understand is an implementation of esri.IdentityManager).  I will try this and report back if it works.


texts: QGIS: Plugin, that can import a textfile, to open a QGIS-project

texts:
Is there a plugin for QGIS that is able to do the following steps??
QGIS imports an external text-file, opens the project and zooms to the (only) polygon that contains „xyz“ in „field1“.
the external file contains the following information:
name of the qgis-project
name of the layer (e.g. polygon-layer) – for example layer1
name of a field in this layer – for example „field1“
content of this field - for example „xyz“


texts: How to make black areas tranparent in GeoWebCache when using Geoserver to cache multi band imagery?

texts:
I am using Geoserver to cache 4 band (RGB,Alpha) satellite imagery.
The native projection of the imagery is NAD83 and I want them in WGS84, when I do the reprojection,using Arcmap, depending on the extent of the imagery, black borders or areas appear in the results. I can make the black areas transparent in Arcmap, but no matter what I do (Set Null, change no data value to zero, clip the black borders,etc) in Arcmap, when I cache the imagery using Geoserver, the black areas re appear in the results both in Open layer viewer of Geoserver and when I serve them as WMTS in Arcmap or QGIS.


texts: How to remove a specified distance from the end of a polyline?

texts:
I am relatively new to writing scripts in QGIS. I was wondering if pyqgis has any inbuilt functions that could take a distance parameter and remove that distance from the end of a target polyline.

I wasn't able to find what I was after in the pyqgis cookbook. I have used Shapely before and could manufacture a function, but I was hoping someone might know of pre-existing one.


texts:
I am using ArcGIS 10.2.

Does anybody know where I can obtain a shapefile or similar for the Irish National Grid? I need a 1 km vector grid with the associated grid references given within an attribute field.


texts: Transforming two different Geodetic CRS in pyproj using python?

texts:
I failed to transform two different Geodetic CRS

from pyproj import Proj, transform

#Transform from Arc1950 to WGS84
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:4209')
x1,y1 = 27.04892, -13.30552
x2,y2 = transform(inProj,outProj,x1,y1)
print x2,y2

but the return:

27.04892 -13.30552

this is not a correct result, it should be  27.04906, -13.30400.(I got this correct result from Mapsource.)


texts: How to find the cell location index of a raster using Lat long information?

texts:
I have a raster for which I want to find the Cell Location (Row, Col) if I have the Lat-long information.

I know how to get the value of a particular cell using lat long.

extract(ras,SpatialPoints(cbind(-82.8,35.2)))

But I need the row-col of this cell, rather than it's value. I am using the raster package in R.

Raster: https://www.dropbox.com/s/8nhfirxr2hm3l4v/fdr_fb.tif?dl=0


texts: Edit exporter.py file in plugin qgis2web?

texts: I want to change a small piece (a string) of "exporter.py" code in qgis2web Plugin, but when I change and save the .py file on "C:\Users...\python\plugins\qgis2web", changes doesn't apply. I also tried compiling the .py file with OSGeo4W Shell, but changes still doesn't apply.


texts:
How can run this script with Sentinel 2 images?

// Use a DateSlider to create annual composites of this collection.
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1');
// Use the start of the collection and now to bound the slider.
var start = ee.Image(collection.first()).date().get('year').format();
var now = Date.now();
var end = ee.Date(now).format();

// Run this function on a change of the dateSlider.
var showMosaic = function(range) {
  var mosaic = ee.Algorithms.Landsat.simpleComposite({
    collection: collection.filterDate(range.start(), range.end())
  });
  // Asynchronously compute the name of the composite.  Display it.
  range.start().get('year').evaluate(function(name) {
    var visParams = {bands: ['B4', 'B3', 'B2'], max: 100};
    var layer = ui.Map.Layer(mosaic, visParams, name + ' composite');
    Map.layers().set(0, layer);
  });
};

// Asynchronously compute the date range and show the slider.
var dateRange = ee.DateRange(start, end).evaluate(function(range) {
  var dateSlider = ui.DateSlider({
    start: range['dates'][0],
    end: range['dates'][1],
    value: null,
    period: 365,
    onChange: showMosaic
  });
  Map.add(dateSlider.setValue(now));
});


texts:
I am new to geodata processing.

For my sample project I am trying to solve next problem.

I have array of geo points. If we place it on plot we can visually see that there are areas where density is higher. I need to find this areas.

I heard about R-tree's and bounding box, but can someone advise me on what way should I go and what material should I read?


texts:
I'm new to GIS. Is machine learning used much in this field, and where?

I'm investigating a job opportunity with a mapping company working on traffic awareness. My background is in supervised machine learning (in text classification) and some use of Hadoop, and they assure me that my skills would be appropriate. I've been interested in GPS and GIS for some time, but I want to make sure that my machine learning skills (which I've invested in heavily -- including an MS degree) would be put to good use.


texts:
I want to use the importer extension for the geoserver. I 'installed' this library and it works so far. My challenge is that i want to upload several shapefiles from a directory and I found this example.

But every time I run this command:

C:\WINDOWS\system32>curl -u admin:geoserver -XPOST -H "Content-type: application/json" -d @import.json "http://127.0.0.1:8080/geoserver2110/rest/imports?async=true"

I get this error message:

Warning: Couldn't read data from file "import.json", this makes an empty POST.
A JSONObject text must begin with '{' at character 0 of

For me is unclear, where I have to place the import.json file?


texts:
I know that I can use arcGIS to display the USGS national map viewer (aerial, raster, topographical, etc), but I will be using API for commercial use and it would cost $4,000.

1) Can I use openlayers to display the USGS national map?  If so, please provide a link.

2) Does USGS have an API?  If so, please provide a link.

3) If there are no free methods of displaying the USGS national map, what is the cheapest method to accomplish this?


texts: Spatial join - allocate values from joined polygon based on input polygon

texts:
I'm trying to determine population inside Special Flood Hazard Areas. To do this I'm doing a spatial join of census shapefiles with the FEMA National Flood Hazard Layer, using the Boolean field SFHA_TF to identify areas inside the Special Flood Hazard Area. The map shows that the SFHA cuts through tracts, taking pieces of some, all of others. But when I join the two tables, the resulting attribute table contains multiple copies of the entire tract; it does not split the tract into pieces (for example a quarter of the tract if the SFHA occupies 1/4 of the tract). Is there a way in QGIS to do a spatial join between two polygons that will get proportionate shares of the joined polygon?
.


texts: Python script to set another layer inactive & select point by outside of boundary line?

texts:
I have two layers, "IDN_ALL" as Polygon and "Point_Link" as point.
"Point_Link" is the active layer. I want to select point data outside polygon data.
Code:
Layer1 = iface.activeLayer()
Layer2 = ??
layer1.selectByexpression("\"?..")

How can I select points outside polygons?


texts: How to resize the Legend Item Wndow in Map Composer?

texts:
I am just wondering if there is any simple way (no messing with codes) to adjust this Legend Item Window in Map Composer. I am using QGIS 2.14 in windows 8.

A much larger (taller) window would really help, rather than this 2-cm window that is so hard to scroll when you have 20+ legend items to edit.

Here is the screenshot of the Map Composer


texts: QGIS Measure area below contour to develop elevation-area table

texts:
Is there a plug-in or methodology to measure areas under contours, which can be used to develop elevation-area tables?
I know you can create contours as lines from DEMs, but you need to convert them to polygons to measure areas. You also need to trim/crop them with a mask to confine the query to the area of interest.
What I would like to do is draw a polygon around an area of interest and ask a plug-in or method to generate a table of elevations v areas that can be used in hydraulic or hydrological modelling to determine the storage volumes of reservoirs, lakes or dams. (Tables of areas are more useful than tables of volumes because you can apply evaporation to the surface area of the water, although you can differentiate from volume to area, or integrate from area to volume.)


texts: get feature on mouse click event on active layer

texts:
I am trying to write a python plugin and want to provide functionality similar to the "Identify feature" tool that already exists in QGIS. Plus the possibility to open the document relied to a feature by getting the path assigned to it in a column. How does one get, after launching the plugin, the clicked feature on mouse click event?

I've used some code from here but I couldn't get it to work.

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *

class IdentifyGeometry(QgsMapToolIdentify):
 def __init__(self, canvas):
  self.canvas = canvas
  QgsMapToolIdentify.__init__(self, canvas)

 def canvasReleaseEvent(self, mouseEvent):
  results = self.identify(mouseEvent.x(),mouseEvent.y(), self.TopDownStopAtFirst, self.VectorLayer)
  if len(results) > 0:
   self.emit( SIGNAL( "geomIdentified" ), results[0].mLayer, results[0].mFeature)


texts: Working with polar projections in Mapbox

texts:
My problem relates to how polar projections can be 'recreated' within Mapbox, although officially Mapbox only supports the Web Mercator projection. A working example of how it can work can be found here:

https://blog.mapbox.com/mapping-arctic-sea-ice-in-a-polar-projection-a8fd68a03c53

However, I am unclear from their blog description how GeoJSON files were successfully incorporated onto the map. For background layers such as a coastline I am able to reproject and assign a Web Mercator projection (ogr2ogr -a_srs EPSG:3857) and upload these shapefiles to Mapbox as tilesets fine.

Applying the same process to convert ice extent shapefiles into GeoJSON files basically outputs a json file with coordinates written in metres (e.g., ogr2ogr -f GeoJSON -a_srs EPSG:3857 B3C__2273_WM.json B3C__2273_PS.shp). When uploading to Mapbox an error message indicates that the units are not acceptable and should be degrees (Input failed. Datasets don't support features outside longitude +/-180, latitude +/-90.).

Example data:
{
   "type": "FeatureCollection",
    "features": [
    { "type": "Feature", "properties": { "Id": 0 }, "geometry": { "type":
    "MultiPolygon", "coordinates": [ [ [ [ -313218.45268004813,
    -2070401.2551701278, 259.16299999999649 ], ...

My question then is how is it possible to import GeoJSON files into Mapbox that are projected in something like EPSG:3995 (polar stereographic), such as in the link above? I would like to manipulate the json data with a time slider, so I don't think uploading as shapefiles is a viable option in this case.


texts: Troubleshooting creation of equal-area fishnet inside rectangular polygon

texts:
I have a set of polygons that are rectangles (10 km by 1 km dimension), and I'm using Fishnet to create smaller polygons (which I'll call "boxes"). The boxes need to be equal area, so I am working in a projected datum and defining the cell sizes in Fishnet. I'm doing test run with 100 meters by 100 meters, but I'll ultimately being looping through a bunch of sizes. The test case here is also only for one polygon, but finale shapefile for analysis will have multiple polygons.

For some reason, Fishnet is not fitting the boxes in my polygon as I expected. It's overshooting on the north side by 5 boxes (500 m) and undershooting the east side by 100 meters (screenshot below; north is up.) I don't understand why this is happening? There must be something about how Fishnet works that I am missing.

In the code below, I first extracted the polygon vertices ABCD to define the Y axis for the Fishnet. This was done because I'll ultimately have many polygons with different Y axis orientations.

NOTE: this question is edited from an earlier question. Earlier question has been deleted and this question has been cleaned up, because commenters complained about defining variables. So, code below is cleaned up with fewer objects.

import arcpy
arcpy.env.overwriteOutput = True

mainBox2 = r'C:\Users\X\_temp working\box.shp'
centerLine2 = r'C:\Users\X\_temp working\center_line.shp'
working_epsg = 26911
outpath = r'C:\Users\X\_temp working'

############################################################

arcpy.FeatureVerticesToPoints_management(mainBox2, outpath + '\\' + 'temp_mainBox_vertices.shp', "ALL")

# Define polygon vertices to match vertices as shown in screenshot
with arcpy.da.SearchCursor(outpath + '\\' + 'temp_mainBox_vertices.shp', ['OID@', 'SHAPE@XY']) as rows:
    temp = dict(r for r in rows)

D = temp[0]
A = temp[1]
B = temp[2]
C = temp[3]

############################################################

### Set local variables
outFeatureClass = outpath + '\\' + 'TempTest.shp'
cellSizeWidth = '100'
cellSizeHeight = '100'
originCoordinate = '%s %s' %(A[0], A[1])
yAxisCoordinate = '%s %s' %(B[0], B[1])
oppositeCorner = '%s %s' %(D[0], D[1])

### Execute CreateFishnet
from arcpy import env
env.outputCoordinateSystem = arcpy.SpatialReference(working_epsg)
arcpy.CreateFishnet_management(outFeatureClass, originCoordinate, yAxisCoordinate, cellSizeWidth, cellSizeHeight, "", "", oppositeCorner, 'NO_LABELS', '#', 'POLYGON')

#


texts: Seeing changes made in connected folder using ArcGIS Desktop?

texts:
I am working on a project in ArcGIS Desktop and I have a connected folder.  In it are multiple files. If I change a file in the folder (add or delete) it seems that I need to close my project and then re-open it to see the changes in the folder.

Is there a way to by-pass this?


texts: Using kernel density estimates as attributes for point data?

texts:
I have spatial point data (say addresses of all inhabitants of a metropolitan area or larger region) with context data for each point like age, gender, employment status, and so on.

For each spatial point, I would like to generate aggregate information about the "neighborhood" environment (with regard to age, gender, employment status), using kernel density estimation (KDE).

E.g., I would like to generate the share of unemployed persons in the vicinity of the spatial point.

So far I can do the kernel density estimation, generate a polygon map and attach the corresponding estimate to each coordinate pair of the point data. My research aim is use the characteristics of these small-scale neighborhood Environment as explanatory variables in regression models at the Level of individual persons.

Two issues are troubling me, though:

I have point data from denser city centers but also more rural areas. How can I generate comparable KDE-estimates for the whole study area?

Using a fixed bandwidth for the whole study area will be suboptimal for all contexts.

Will using adaptive bandwidths generate comparable estimates, though?

So far I only receive the estimated number of unemployed from the KDE.

Since this is a function of population density, I guess I have to divide the estimated number of unemployed by the estimated number of the total population, right?

If I do this, do both KDE have to use the exact same bandwidth to be comparable?

Some events like unemployment might be rare (= less dense point patterns) and thus would need to have a larger bandwidth, while the spatial patterns for the whole population should be denser and thus need a lower bandwidth. A common bandwidth could thus be suboptimal for one or both.

Should I compromise with a mean bandwidth or pick the largest common bandwidth or are KDE from varying bandwidth still comparable?

I hope I could convey what I want to do and why I want to do it.
MDo you think my approach is feasible or do you have an alternative proposal?

I am new to KDE.


texts: Issue with a point insertion into a Geoserver points layer. WFS-T

texts:
What is the best way to insert a feature from the client side into an existing GeoServer layer? Maybe using an AJAX call and JSONP? Or OpenLayers (i can't figure out how to have the proxy working)?

If someone have an example, I'll be very grateful

UPDATE

Actually what I'm trying is:

1- Build WFS-T formatted string in this way (reduced version):

<wfs:Transaction
 <wfs:Insert>
  <sde:uruguayos>
   <sde:the_geom>
     <gml:Point srsDimension="1" srsName="http://www.opengis.net/gml/srs/epsg.xml#3857">
       <gml:coordinates decimal="." cs="," ts=" ">-33.7291,-56.3353820</gml:coordinates>
     </gml:Point>
   </sde:the_geom>
  <sde:Name>Uruguayo Nuevo</sde:Name>
 </sde:uruguayos>
</wfs:Insert>

After that, i'm sending a POST with the data above via ajax to GeoServer and now i fall out into a cross-domain-issue (and make sense).

jQuery.ajax({
  async: false,
  dataType: 'text/xml',
  type: "POST",
  url: wfs,
  data: postData,
  //TODO: Error handling
  success: function(xml) {
    //TODO: User feedback
    alert("Ready...");
  },
  failure: function (xml) {
          alert("Something went wrong in the request");
  }
});

UPDATE 2

I could solve the Cross domain issue. And now I'm with a strange problem making the insertion over my PostGIS Layer (through GeoServer).

As I said, I have this xml string:

var wfs = "http://"+ip+":"+port+"/geoserver/ows?service=wfs&request=Transaction";
var wfsEncoded = encodeURIComponent(wfs);

// INSERT point

var postData =
    '<wfs:Transaction service="WFS" version="1.0.0" xmlns:wfs="http://www.opengis.net/wfs"     '+
  'xmlns:sde="http://geoserver.sf.net" '+
  'xmlns:gml="http://www.opengis.net/gml" '+
  'xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" '+
  'xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.0.0/WFS-transaction.xsd '+'http://www.openplans.org/topp '+ wfsEncoded +'/DescribeFeatureType?typename=sde:uruguayos">'+
  '<wfs:Insert><sde:uruguayos><sde:geom> '+
        '<gml:Point srsDimension="1" srsName="http://www.opengis.net/gml/srs/epsg.xml#3857"> '+
         '<gml:coordinates decimal="." cs="," ts=" ">-33.7291,-56.3353820</gml:coordinates> '+
         '</gml:Point> </sde:geom> <sde:Name>Uruguayo Nuevo</sde:Name> </sde:uruguayos>          </wfs:Insert> '+
    '</wfs:Transaction>';

I make two tries with a different wfs var (url in the code)

var wfs = "http://"+ip+":"+port+"/geoserver/ows?service=wfs&request=Transaction";

And

var wfs = "http://"+ip+":"+port+"/geoserver/wfs";

The first one got a response that says:

 <ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows/1.1"  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="2.0.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://192.168.38.229:8080/geoserver/schemas/ows/1.1.0/owsAll.xsd">
  <ows:Exception exceptionCode="MissingParameterValue" locator="request">
    <ows:ExceptionText>Could not determine geoserver request from http request org.geoserver.platform.AdvancedDispatchFilter$AdvancedDispatchHttpRequest@1d08ec4</ows:ExceptionText>
  </ows:Exception>
 </ows:ExceptionReport>

And the second response from Geoserver is like that:

 <ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows/1.1"         xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="2.0.0"      xsi:schemaLocation="http://www.opengis.net/ows/1.1     http://192.168.38.229:8080/geoserver/schemas/ows/1.1.0/owsAll.xsd">
   <ows:Exception exceptionCode="NoApplicableCode">
     <ows:ExceptionText>Could not find request reader (either kvp or xml) for: net.opengis.wfs20.TransactionType</ows:ExceptionText>
   </ows:Exception>
 </ows:ExceptionReport>

The URL encoding tip was taken from Data request to GeoServer: "Could not determine geoserver request from http request GET /geoserver/ows?service=WFS HTTP/1.1"


texts: Mapbox markerLayer from GeoJSON with prompt as input not working

texts:
I am trying to display a pin on a map. It works fine if I code the GeoJSON data into the script of my file. If I try to insert the GeoJSON data via a user prompt, it doesn't work. The variable must be somehow different when it comes from a prompt then when it is coded in the file. What is the difference and how can I fix it?

GeoJSON pin

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Point",
        "coordinates": [
          -104.0625,
          39.36827914916011
        ]
      }
    }
  ]
}

HTML file

  <link href='http://api.tiles.mapbox.com/mapbox.js/v1.6.0/mapbox.css' rel='stylesheet' />

    <button onclick="myFunction()">Enter GeoJSON</button>


texts: It it possible to uninstall IDLE from a ArcGIS 10.0 install?

texts:
I don't use the IDLE environment that was installed with ArcGIS's python and in particular the shell extension for .py files (Edit with IDLE) is an eyesore.

I was on my way to removing this shell extension from showing up on the context menu, but I was curious if I could just get rid of the entire IDLE installation. From what I understand its part of Python, but its rarely included with Linux package of python.

Is this possible to do away with IDLE in Windows/ArcGIS or are there unintended consequences?


texts:
I am trying to join attributes from 2 shapefiles. 1 shapefile is a polygon layer indicating police precincts and the second is a point layer with crime statistics for each police station. The police station and the police precinct are the same.
The police station layer that contains the crime statistics has roughly 40 points per police station. See attached police station image

.

I want to join this layers table to the police precinct layer (See attached precinct layer image)

The trouble I am having with this, is that when i use the Join by location tool, it only takes the first 17 records from the police stations layer and appends the data to the existing 17 police precincts.
I need it to append all 40 crimes points to the precincts layer without creating 40 polygons that lie ontop of oneanother. Just 1 polygon with multiple crime records.
Does anyone know how to do this in QGIS?


texts:
In my practice I use a lot WMS and WFS in Geoserver. I have not found a document describing the concept of a bounding box(bbox). ArcGIS has map extent properties to identify map boundary. Google Map/Bing map has getbound() method to identify the map view. However, when defines a bounding box in WMS or WFS, it seems to involve projections parameters and map view extent. The parameters in an HTTP GetMap request looks confusing to me since if I use different projection systems the units in the bounding box varies significantly, and the incorrect setup of bbox would cause rendering fail So, could anyone point me to the document explaining the bbox?
Thank you!


texts: Avoid multiple labels for the same layer in QGIS2WEB export

texts: When I export my QGIS generated shapefile with layers to a web map using qgis2web, the output map is showing repeated number of labels for the same layers. In my QGIS desktop, the label appears only once for each layer but once exported into a map, it shows multiple labels for the same layer. How can I solve this issue?


texts: How to find out whether a satellite image is atmospherically corrected

texts:
I have a satellite image and I am not sure whether an atmospheric correction has been done with the image. Is there a way to find that out by looking at some values?

A atmospheric correction was supposed to be done with the semi automatic classification plugin. The output sould be a layerstack of each Band. So I am not sure wether it is atmospherically corrected. I haven't done the atmospheric correction, I just got the outcome Image.

I added a spectrum view of the satellite image of which I am not sure whether it has been atmospheric corrected or not:


texts:
I have remote GeoServer with published GeoTIFF's which have data security via login/password. Initial size GeoTIFF about 5-10Gb. Is there any way download raw GeoTIFF file from GeoServer (using Python script, web URL link ...)?

Updated 26.02.2016

Example WCS URL for download raster files:

var workspace = layer.params.LAYERS.split(':');
var downloadurl = "http://localhost:8080/geoserver/" + workspace[0] + "/wcs?service=WCS&version=2.0.1&request=GetCoverage&CoverageId=" + layer.params.LAYERS + "&format=image/tiff"
window.open(downloadurl);


texts: Changing cell size output for Proximity (Raster Distance) tool in QGIS

texts:
I am trying to create a raster layer with distance to shore that has a resolution of 1km.

I've created a raster distance layer using the Raster > Analysis > Proximity (Raster Distance) tool in QGIS which worked however the difference between each section is 15km. I need the cell size to be more accurate as I will be using it in a species distribution model.

Is there any way to do this on QGIS?


texts:
I am trying to solve a network generation problem, and I'd be happy about inputs. First of all, my problem description:

I have a list of connections between locations with respective distances. For example:

LocationA <-> LocationB: 4500 m
LocationB <-> LocationC: 3000 m
LocationC <-> LocationA: 2000 m
...and so on...

What I do not know, but try to find out: The actual geographic position of my locations. I do know the position of a few nodes in the network, and given a large list of connections with distances, I should be able to approximatively find the positions of my unknown locations.

Note: I am not trying to actually to optimize any paths within this network (e.g. Dijkstra etc.). I simply would like to know what geographic coordinates the nodes in my system have.

My question is: What sort of "node placement algorithm" am I looking for? I'd be happy for any keywords in this regard. I am 100% sure there are some algorithms for this problem, but I do not even know what to search for.


texts: Geoprocessing Service Cannot Download Result?

texts:
I have implement the ArcGIS JavaScript example just like the ESRI official example (clip and ship) .Everything was going right except the final step ,when the gp job have finish and my browser lead me to the result URL it have show the 404 Error.
I checked out the gp job directory which do contain the job results(output.zip).At the same time i had add the server's hostip and hostname in my hosts file.
My  Environment: Window 2003 Sp2 ,ArcGIS Server 10, Chrome 18;
I am using the java version which the web server is tomcat 6.


texts: Using con statement in raster calculator to include raster as condition in existing raster?

texts:
I'm doing an assignment where we have to use a Weighted overlay to generate a raster of suitable places to plant crops. I've put all the factors in the weighted overlay and generated my output but now I need to consider the constraints. These are built up areas and conservation areas that need to be classified as 'not suitable'. The lecturer says that we need to use a CON statement in the raster calculator to get the '0' values from the constraint's raster put onto the weighted overlay raster.

But, I have no idea how to do that. All the resources I have found use the format CON (in raster = 0, x, y) but I don't know how to fill the Y value with an entire raster (the Weighted raster)... if that makes sense?


texts: Resizing point on zoom out using OpenLayers

texts:
I need to resize a point on zoom.  When you deep zoom the point get larger and on zoom out the point is smaller.

Resizing point on zoom out seems to have exactly the answer.

@CaptDragon posted an image of exactly what I need and then a demo, HOWEVER the demo is now dead.  I've looked into the API that @Mapperz mentions but for whatever reason cannot get my point to resize larger when zooming in.

Can anyone provide a working example similar to that of @CaptDragon?


texts: How to get the height according to EGM96 for a point?

texts: For a program I have to calculate the height according to EGM96. I found http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html - this site delivers files with parameters and a fortran-program. But I feel unable to understand what is going on in this Fortran-program and have no idea how to use the parameters to generate the correct height. Can someone explain the algorithm to me?


texts:
I have Shapefiles without a CRS (e.g. without .prj file).
I would like to set the CRS (EPSG:31468) before I pass the Shapefiles.
I tried the following:

Setting a .prj file manually with the string

PROJCS["DHDN_3_degree_Gauss_Kruger_zone_4",
    GEOGCS["GCS_DHDN",
        DATUM["D_Deutsches_Hauptdreiecksnetz",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]
        ],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]
    ],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",12],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",4500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]
]

Using ogr2ogr

ogr2ogr -s_srs EPSG:31468 -t_srs EPSG:31468 out.shp in.shp

If I import those shapes into QGIS, QGIS always recognizes CRS EPSG:2167.
I have set in QGIS the correct CRS and exported it in a new file. The .prj of the new file is identical to that one of the input file.
So what's wrong with my files?!


texts:
I want to implement a sort of content negotiation so as to export my map in formats like pdf jpg etc.

Is there any way to export my current map view (with all the enabled options so far) in various formats? With the term current map view i mean "what the end user sees at a specific moment".

My map is generated by combining map Services from different web servers in the Web.

As a result an option such as combining my services in arcMap, publish my map as a MapService and use the export capability provided cannot help:
Ex. http://sampleserver6.arcgisonline.com/arcgis/rest/services/USA/MapServer/export?bbox=-241.20512466608835,-27.80850864596666,51.688012394341115,67.48206676752923

I want to get the data live directly from the servers.

In fact would like to use a function like myMap.export("format") where myMap is the current map object shown by my browser.


texts: Selecting enclosed features within same layer using QGIS?

texts:
My problem: I have a single shapefile from a habitat-survey and need to select all hedge-features that are enclosed in meadows.The features that are to be selected were created by filling rings.
Select by attribute does not work since not all hedges are within meadows.
I've tried the Spatial Query Plugin's 'within' and 'contains' option after duplicating the layer, but since the features do not overlap, the plugin finds no result. The 'touches' option gets results, but unfortunately also selects neighbouring features which I don't need.

How can I select only the enclosed objects without having to do it manually?


texts: Updatecursor returns with 'variable not defined'

texts:
I am fighting a bit with the cursor function since this morning.
Thanks to the help of the people here and some tutorials I found I got this far

import arcpy

arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Ecodis_clip_split" #set workspace to folder with values
value = [] #create list for the returning value of searchCursor
eco = arcpy.SearchCursor ("Anagance.shp", "", "", "area_ha", "")
#search for the value in attributetable
for row in eco:
value.append(row) #write the value into the list

arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Wetland"
#set workspace to folder with values
diss = arcpy.UpdateCursor ("diss_1.shp", "", "", "area_eco", "")
for row in diss:
    row.setValue(area_eco, row.getValue(value)) #identify all the rows and Update die field area_eco with the extracted value from the search cursor out of the list
    cursor.updateRow(row)

del eco
del diss

It bounces back with

area_eco is not defined

Which means - if I understood the concept of python correct - that this is not a variable. Ture, its my field which needs to updated.
What am I doing wrong?


texts: Multiplying rasters by integers in processing modeler

texts:
I need to multiply several rasters by integers as a part of a multiple steps process. For this reason I'm using QGIS' Processing Modeler.

I'm looking for something similar to the Raster calculator, but it doesn't seem to be available from the Processing Modeler.

I've also tried to use the GRASS r.mapcalculator algorithm (with A my input raster layer, and applying the formula A*100 for example), but it always returns the following error message:

Oooops! The following output layers could not be open
Output raster layer: /tmp/processing/ad0f2935615a472eb6b6b6feb44f3e72/outfile.tif


texts: Is there limit on size of GeoJSON files that can be handled by QGIS?

texts:
I'm trying to load a GeoJSON file with about 5K features in QGIS, and it fails, saying Invalid Data Source without further details. If I remove features down to 1K, it works. Is there a limit on the size of GeoJSON files that can be handled by QGIS? How can I understand what's going on?

I'm using QGIS version 2.8.2-Wien on Mac.


texts: Set or change the geometry type of an empty point shapefile using Python with GDAL/OGR

texts:
I'm automatically creating a point shapefile from a postgreSQL point table (having a geometry field called "geo" defined as geo geometry(Point,4258)) for a spesific flatenr (%%a) using

ogr2ogr --config GDAL_DATA "C:\OSGeo4W64\share\gdal" -f "ESRI Shapefile" C:\kmfelt\shape\utm32\2017\%%a_p.shp PG:"host=%Server_host% user=%Schema_user% dbname=%Database% password=sl_nab" -sql "SELECT * FROM kmfelt_o2_p where flatenr=%%a AND komid < 1800" -t_srs EPSG:%SRID% -geomfield geo

followed by some editing of the shapefile using python

call python %prog%\05_deletefields_p.py C:\kmfelt\shape\utm32\2017\%%a_p.shp

In general, this works fine and results in an output point shapefile.

However: Some times it turns out that the shapefile is empty, e.g the shapefile does not contain any points since the postgreSQL point table does not contain any data for the specific flatenr.

Example: flatenr = 2060 does not exist in the table

An empty shapefile is created, but the geometry is set to linestring?!?

From 05_deletefields.py

...
# Set POINT-shapefile
fn = r""+sys.argv[1]+""   # Read filename from input
print('filnavn er ' +fn)

#get the driver
driver = osgeo.ogr.GetDriverByName('ESRI Shapefile') # will select the driver for our shp-file

# open the data source
dataSource = driver.Open(fn, 1) # 0 means read-only. 1 means writeable (read/write)

# get layer
layer = dataSource.GetLayer()
geotype = layer.GetGeomType()
print(geotype)

This results in (output on the screen)

  filnavn er C:\kmfelt\shape\utm32\2017\2060_p.shp
  2

The 2 indicates that the geometry type is wbkLineString. Why?
(This is confirmed in QGIS: the shapefile is shown with a line symbol in the Layers panel and when loading a point style file an error message appear telling QGIS cannot apply style to layer with a different geometry type.)

How may I change/set/force the geometry type to be point (wkbPoint)?

So far I have be using this code the py-file in case of no points

...
if layer.GetFeatureCount() == 0:
    km = ogr.FieldDefn('KMNR', ogr.OFTInteger)
    km.SetWidth(4) #4 integer width
    print "Legger til egenskapen KMNR som siste kolonne"
    layer.CreateField(km)
    # After adding fields, get the layer definition
    featureDefn = layer.GetLayerDefn()
    # Create a feature from scratch
    feature = ogr.Feature(featureDefn)
    # Populate the attribute table with a default (dummy) value
    print "Setter KMNR lik 9999"
    feature.SetField("KMNR", 9999)
    # Finally, feature must be written to the layer
    layer.CreateFeature(feature)
    # Clean up memory
    feature.Destroy()


texts: QGIS 2.4 layout view PDF/raster file export moves points

texts:
I created a map, using a Google basemap (OpenLayers plugin) and imported a table as a CSV file (containing a name of a city, latitude, longitude for wgs84 and a row named "amount"). It displays correctly also in the layout view but when I export the Layout the points are moved to the east while the base map stays where it should. I tried to convert the wgs84 coordinates of my CSV file by saving it as a shapefile (defining the coordinate system as wgs84 Pseudo Mercator) but that did not work either.

The points in my map are still moved after exporting from the layout view. I am working on this for 6 hours now and I am getting crazy.

Could anybody please please help me with this problem?


texts: Is there an equivalent "why-upgrade PostgreSQL" tool for PostGIS?

texts:
I rarely keep up with the PostgreSQL and PostGIS releases and it can a challenge deciding on when to do so for the major/minor releases.  The why-upgrade site (http://why-upgrade.depesz.com/show?from=9.5.1&to=9.6.7) helps show what I am missing or risking by not upgrading.

Is there an equivalent site for the PostGIS extension?


texts:
I need to draw an telephone  or energy lines as is done on the polish-core map. Thelephone line was drown by the polyline. I tried the SVG symbols, when I have a lot of single-lines, for each I 'put' symbol svg on the  first and the last vertex. For polylines, unfortunately, line markers located on the central vertexs do not rotate to the next segment, but are suspended 'in the air'.

If someone has already tested this, I'd like to know the rules.
I enclose a snapshot of my unsuccessful endeavors.

And below is the effect that I need, (now its generated by Inkscape)

Thanks for the reply.


texts: How to export vector tiles from a PostGIS raster database?

texts:
I have a PostGIS raster database with 5 tables loaded from 5 geoTIFF images, tiled 100x100. I am trying to setup a map tile server to display this data on an openlayers map. I need to transform the raster data into geometries and then export as vector tiles. I have the following query to select the data expecting a geometry:
WITH
            bounds AS (
                SELECT ST_Segmentize(ST_MakeEnvelope(xmin, ymin, xmax, ymax, 3857),seg_size) AS geom,
                       ST_Segmentize(ST_MakeEnvelope(xmin, ymin, xmax, ymax, 3857),seg_size)::box2d AS b2d
            ),
            mvtgeom AS (
                SELECT ST_AsMVTGeom(ST_Transform({some_geom}, 3857), bounds.b2d) AS geom,
                       gid, name, type
                FROM table t, bounds
                WHERE ST_Intersects({some_geom}, ST_Transform(bounds.geom, 3857))
            )
            SELECT ST_AsMVT(mvtgeom.*) FROM mvtgeom

How can I get the correct slice of the raster to input as a geometry into this query? I know I can use ST_PixelsAsPolygons to get polygons from the individual pixels, but how can I bound the data and give that query what it needs?
I tried adding:
fullrast AS (SELECT ST_Union(rast) As rast FROM table),
rast_clip AS (
    SELECT ST_Clip(fullrast, 1,
        bounds.geom)),

But it looks like the unioned raster exceeds the size limit:
ERROR:  invalid memory alloc request size 18446744072322304320


texts: How to exclude/eliminate invalid areas (such as NoData, cloud, water, etc) from raster so that it does not effect the accuracy assessment?

texts:
I have two classified raster (snow cover % 10-100) one is prepared from coarse resolution and another is from finer resolution. I am doing accuracy assessment and matrix calculation of snow cover %. However, both classified raster has invalid classes like NoData (-9999), cloud (value=205), water (210), permanent ice(215).

Now both the raster has same pixel size after resampling and covering same area by using gdalwarp

These classes/values are effecting the accuracy assessment. I only want to comapre both raster in terms of snow cover %. How can I eliminate those value from both raster??

My code

#reading in array Snow cover % raster from true raster (finer resolution)
img= rasterio.open(filepath+'20020310_ndsi_reclass.tif')
img_r=img.read(1)

#Reading in array Snow cover % from prediction raster (coarser resolution)
img_pre= rasterio.open(filepath+ '20020310_ndsi_reclass_pred.tif')
img_pre_r=img_pre.read(1)

#making it 1dimesional
y_true= np.ravel(img_r)
y_pred= np.ravel(img_pre_r)
y_true = y_true.astype('int')
y_pred = y_pred.astype('int')

from sklearn import metrics
print(metrics.classification_report(y_true, y_pred))

Matrix calculation

from sklearn.metrics import confusion_matrix

cm=confusion_matrix(y_true, y_pred)

import seaborn as sn
ax= plt.subplot()
sns.heatmap(cm, annot=True, ax = ax); #annot=True to annotate cells

# labels, title and ticks
ax.set_xlabel('Predicted labels');ax.set_ylabel('True labels');
ax.set_title('Confusion Matrix');
ax.xaxis.set_ticklabels(['nodata','0','10', '20', '30', '40', '50', '60', '70', '80', '90', '100', 'cloud', 'water', 'ice']); ax.yaxis.set_ticklabels(['nodata','0','10', '20', '30', '40', '50', '60', '70', '80', '90', '100', 'cloud', 'water', 'ice']);

How to exclude  Nodata, cloud, water, and ice from both the report?


texts: How to run Python script in QGIS via plugin

texts:
I am an absolute beginner in QGIS no idea in this field.
I just created a Python plugin with pluging builder in QGIS successfully

PYTHON_PLUGINS_PATH/ MyPlugin/ ----------------- inside it contain these files

       ____init__.py
       mainPlugin.py
       mainPlugin_dialog.py
       mainPlugin_dialog_base.ui
       plugin_upload.py
       README.txt
       metadata.txt
       resources.qrc
       resources.py
       resources_rc.py
       .
       .
       etc

I just want to write some code like print("hello world").
I want to know should I write it inside the mainPlugin.py file If yes then ok, otherwise tell me where should write my code.

After writing my code how can I run  It through my plugin which I created
in QGIS. I am seeing only my plugin name in QGIS my plugin has no option to import file to take this file  mainPlugin.py and run it.

What is the procedure to run my mainPlugin.py file through my plugin inside QGIS?


texts: How to calculate new point coordinates from distance and angle using ArcGIS SDK for iOS?

texts:
I want to calculate the coordinates of a new point.

What I have: The source point coordinates, the distance (meters) and the angle

For example: I want to create a point within 100 meters and 45 degree from my source point

I use wgs84 as coordinate system. My question now is if it's possible to calculate a new point with the given attributes above and when its possible how I can do that.

I use the ArcGIS iOS SDK, but I haven't found any method in there to solve my problem.


texts: Improve performance when iterating features in QGIS and using "for" loops

texts:
I am trying to modify attributes of features in QGIS comparing attributes with list'elements.
First I define new attributes that need to be completed:

layer_intersected.dataProvider().addAttributes([QgsField("ID_acc",QVariant.Int,
    QgsField("Nb_acc", QVariant.Int)])
layer_intersected.updateFields()

Then, I iterate over my features:

layer_intersected.startEditing()
for feature in layer_intersected.getFeatures():

I have two lists containing elements.
Now, I want to test if the first attribute ("NOMBRE") is in one of the extracts of "list1" (the list is divided into a certain number of elements, each part of elements is  an 'extract' of list) and if the second attribute ("nombre_2") is also in the same extract of "list2":

for i in xrange(0, len(list1), len(list_streets_network)):
    l_streets_t1 = list1[i:i+len(list_streets_network)]
    set_l_streets_t1 = set(list1[i:i+len(list_streets_network)])
    set_l_streets_t2 = set(list2[i:i+len(list_streets_network)])
    l_id = list_id_acc[i:i+len(list_streets_network)]
    l_acc = list_nb_acc[i:i+len(list_streets_network)]
    if feature["NOMBRE"] in set_l_streets_t1 and feature["nombre_2"] in set_l_streets_t2:

Then, if it is the case, I complete the attributes created

         index1 = l_streets_t1.index(feature["NOMBRE"])
         feature["ID_acc"] = l_id[index1]
         feature["Nb_acc"] = l_acc[index1]
         layer_intersected.updateFeature(feature)

Finally:

layer_intersected.commitChanges()
return layer_intersected

My problem is the speed of this code. I have to do this for thousands of features, and it takes days to calculate... I am sure there is a faster solution, but I do not know how to improve it...

I had a look on Numpy doc, but did not find the best solution...
If someone has an idea, please help me ! :).

Here is the full code of the function:

> def create_id_for_right_matches(layer_intersected, list1, list2, list_id_acc, list_nb_acc, list_streets_network):
    layer_intersected.dataProvider().addAttributes([QgsField("ID_acc", QVariant.Int), QgsField("Nb_acc", QVariant.Int)])
    layer_intersected.updateFields()
    layer_intersected.startEditing()
    request = QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry).setSubsetOfAttributes(['NOMBRE', 'nombre_2'], layer_intersected.fields())
        for feature in layer_intersected.getFeatures(request):
            for i in xrange(0, len(list1), len(list_streets_network)):
                l_streets_t1 = list1[i:i+len(list_streets_network)]
                set_l_streets_t1 = set(list1[i:i+len(list_streets_network)])
                set_l_streets_t2 = set(list2[i:i+len(list_streets_network)])
                l_id = list_id_acc[i:i+len(list_streets_network)]
                l_acc = list_nb_acc[i:i+len(list_streets_network)]
                if feature["NOMBRE"] in set_l_streets_t1 and feature["nombre_2"] in set_l_streets_t2:
                    index1 = l_streets_t1.index(feature["NOMBRE"])
                    feature["ID_acc"] = l_id[index1]
                    feature["Nb_acc"] = l_acc[index1]
                    layer_intersected.updateFeature(feature)
    layer_intersected.commitChanges()
    return layer_intersected


texts: Preserving special chars using osgeo.ogr.Driver to shapefile in Python?

texts:
I am converting KML files to shapefiles and run in to the problem that special characters are not preserved.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os, ogr, osgeo, sys

ogr.UseExceptions()
# ...
shp_driver = ogr.GetDriverByName("ESRI Shapefile")
# String with file name called outpolyext:
outpolys_file = shp_driver.CreateDataSource( outpolyext )
outpolys_layer = outpolys_file.CreateLayer("polys", srs=srs, geom_type=ogr.wkbMultiPolygon)
outpolys_layer.CreateField(ogr.FieldDefn('Name', ogr.OFTString))
# Left out other fields I created above
# ...
kmldr = ogr.GetDriverByName('KML')
datasource = kmldr.Open(kml_file)
# Create new feature "feature" for shapefile
# with new field names:
ftd = outpolys_layer.GetLayerDefn()
feature = ogr.Feature(ftd)
# ...
# Get first layer from KML, loop through each
# KML feature called "f", then:
print(f.items())
feature.SetField("Name", f.GetField('Name'))
feature.SetGeometry(multipolys)
outpolys_layer.CreateFeature(feature)

In the original KML file, I have e.g. the following value "bæk/sti", the special chars in the resulting shapefile end up blank.

In the terminal print(f.items()) results in:

{'Name': 'b\xc3\xa6k/sti', 'Description': ''}

Can anyone see where it goes wrong - already when reading the KML file or later when saving to the shapefile - or somewhere in between?


texts:
I am trying to spatially join point layers of some 12,000 women interviewed in Egypt to Egypt's polygon layers (administrative districts boundaries) using QGIS 2.2 and need some advice. I am a QGIS self-learner.

I have read many QGIS manuals, but I cannot complete the spatial join successfully. I have added the vector layers of both Egypt's women's GPS points and Egypt's polygons (district boundaries). I used the UTF-8 encoding. Whenever I attempt to join attributes by location, the following message appears "Input layers have non-matching CRS". However, when I check the CRS of both layers, I found it the same (EPSG:4326 - WGS 84).


texts: Cannot use transparency slider with hosted tile layer in ArcGIS Online WebAppBuilder

texts:
I created a hosted tile layer from a raster in ArcMap and uploaded this to ArcGIS online.  I added the layer to my WebMap and within that I can use the transparency tool to adjust transparency of the hosted tile layer without any issues.
However, I created a WebApp with this map as the underlying Web Map and when I go to the layer list and click on the three dots at the end of the layer name to adjust transparency, I am just seeing a small circle of dots which I think must mean the request is timing out.
If anyone has any suggestions as to what might be causing the problem I would be very grateful. It is happening to all the hosted tile layers that I have created, not just this one.
I am not sure why the transparency works in the WebMap and not the WebApp - is there a setting somewhere that I need to change?


texts: creating a transparent layer background in QGIS

texts: I  have created a heatmap in QGIS  to cover my point data and I have saved it as geoTiff so that I can import into tilemill, which I managed to do but my background is black .In QGIS however  I want to make the background of the layer transparent not black so that it shows my heatpoints/areas only . I hahve tried playing around with singleband and changing to singleband pseudocolor but this then removes the heat map coverage and leaves single coordinate points?


texts:
I need to read WFS data from one GeoServer and insert it into another GeoServer, using code (C# in this case).
So I have a FeatureCollection, like this (real values masked):
<wfs:FeatureCollection xmlns= ...and so on... https://xxx.xx/xxxx/wfs?service=WFS&version=1.0.0&request=DescribeFeatureType&typeName=xxxx">

What is required to convert this into a WFS insert/update transaction? They seem very similar.
To sum it up: I need to get data from one GeoServer to another periodically, using code. The source I can only access using WFS. The destination runs in a Docker container. I can add to this docker image, maybe installing some tooling for the job, but I was thinking inserting the data with WFS-T was the simplest solution. Don't know for sure, though.


texts: Uninstalling ArcGIS Desktop 9.2 to install ArcGIS Desktop 9.3?

texts:
I had previously installed ArcGIS 9.2. I uninstalled it to install 9.3. but I am getting error:

  ArcGIS Desktop 9.2 exists on your machine. You must uninstall ArcGIS Desktop 9.2 before installing ArcGIS Desktop 9.3

ArcGIS 9.2 does not show up in list of installed programs. I used revo uninstaller to clean the remaining files and registry as well. ArcGIS 9.3 conflict detector does not show any conflicts but still I am not being able to install ArcGIS 9.3. The same package installs in other machines, just not in the machine with 9.2 previously installed. ArcGIS 9.3 license manager is easily installed.

I must use ArcGIS 9.3 as an extension is built on 9.3 and does not work in recent versions. I am using Windows 7 operating system.


texts: I have the same question as this one posted 3 years ago. I use windows 10 and  USB 3.0. I have seen a lot of angry reviews about well known GPS device brands so I was wondering if anyone has been using GPS tool in QGIS and is willing to share info about type of device being used.


texts:
I want to create something in android like this: http://openlayers.org/dev/examples/wmts-getfeatureinfo.html I'm able to load the costum tiles with WMTS but now i need to get data to show in some kind of a popup window when clicked(touched).

I found that you can get extra info with the service getFeatureInfo, for this request i need some extra parameters like "I" and "J". Whitch are the number of pixels from the bottom-left point of the tile to the clicked position on the tile? But I have no idea how I can get these parameters in android.

Is this the right method to do this or is there a better solution?

Any advice would be appreciated.


texts: Least number of polygons to cover vector in QGIS

texts:
I have two layers :

Polygon rectangles layer representing different footprints
Line layer representing roads

I want to know what is the least number of polygons from the polygon layer do I need to cover the roads layer, meaning basically to remove any polygons that covers the same area as other polygons in its layer. I don't care which one will be deleted and which one will remain, as long only 1 remains per for a certain area.
Another image with opacity
Any idea how to achieve this?


texts: Problem with Search Radius (Bandwidth) of Kernel Density Estimate

texts:
My task is to do a Kernel Density estimation with a cell size of 500m and a search radius (bandwidth) of 500m.
The search radius doesn´t seem to work though.

These are the settings:

This is the output I get using 500m Search Radius:

These are the the values when using the default search radius:

Data Frame Properties:

NAD_1983_Statistics_Canada_Lambert

WKID: 3347 Authority: EPSG

Projection: Lambert_Conformal_Conic

False_Easting: 6200000,0

False_Northing: 3000000,0

Central_Meridian: -91,86666666666666

Standard_Parallel_1: 49,0

Standard_Parallel_2: 77,0

Latitude_Of_Origin: 63,390675

Linear Unit: Meter (1,0)

Geographic Coordinate System: GCS_North_American_1983

Angular Unit: Degree (0,0174532925199433)

Prime Meridian: Greenwich (0,0)

Datum: D_North_American_1983

Spheroid: GRS_1980

Semimajor Axis: 6378137,0

Semiminor Axis: 6356752,314140356

Input Layer

Data Type:  Shapefile Feature Class

Shapefile:  C:\Users\lhark\Desktop\SHK\ArcGis_Proj_Canada\shapefiles\baked_goods_stores.shp

Geometry Type:  Point

Coordinates have Z values:  No

Coordinates have measures:  No

Projected Coordinate System:    NAD_1983_Statistics_Canada_Lambert

Projection: Lambert_Conformal_Conic

False_Easting:  6200000,00000000

False_Northing: 3000000,00000000

Central_Meridian:   -91,86666667

Standard_Parallel_1:    49,00000000

Standard_Parallel_2:    77,00000000

Latitude_Of_Origin: 63,39067500

Linear Unit:    Meter

Geographic Coordinate System:   GCS_North_American_1983

Datum:  D_North_American_1983

Prime Meridian:     Greenwich

Angular Unit:   Degree

Inverse Flattening: 298,257222101