texts: How to create a line layer which connects a point layer to a line layer?

texts:
In the map shown below, I have manually created the dummy lines (in green).
Using python, and without using PostGIS or sextante, I would like to automate the creation of these dummy links between a point and a line layer, based on the following criteria:

Depending on the polygon size, I should be able to add an X amount of dummy links between the centroid (point) and its nearest projection(s) onto the roadway layer.
The script so that the largest polygons get the most amount of dummy links and the smallest get the least amount.
The script should connect the newly created dummy link to nearest highway based on a classifications scheme, i.e. connect to "busier" highways first
If a polygon does not have any roads/links in it, the script should allow the centroid of that polygon to connect to roads outside of its boundaries within Y units of distance
The dummy links must be connected to the roadway layer in such a way that this roadway layer could be potentially used for routing including the dummy links.
Therefore, the dummy link layer must have the same attributes as the road layer
Create a new layer with all the new dummy links

Question:
Is this possible using Qgis and Python alone? any pointers as to where to start or can you help in providing a pseudocode to assist me to solve this problem?


texts:
my searchCursor look like this.

arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Ecodis_clip_split" #set workspace to folder with values

eco = arcpy.SearchCursor ("NAME.shp", "", "", "area_ha", "")
#search for the value in attributetable

print(eco)

It is supposed to give a numeric value (area in ha) but instead it bounces back with this:

geoprocessing cursor object object at 0x02934DE0

I am working in PyScripter and I am wondering if that is my value (which should be numeric, precisely 141216) just expressed differently or something that went totally wrong.


texts: regionating vector data for use in Google Earth

texts:
Does anybody know of free/open source tools to regionate vector (kml) data, like the regionate tool in Google Earth Pro (see http://support.google.com/earth/bin/answer.py?hl=en&answer=175007)?

I have only found one, which is Regionator (http://code.google.com/p/regionator/). It is a set of Python scripts and I am not sure how to get it to work so perhaps if anybody has a link to tutorials or the like?.


texts: Highlighting number of overlayed lines in ArcMap?

texts:
Using ArcMap 10.2, I have a feature class composed of several lines representing bus routes for the same transit line. There are different routes according to schedules.
In the image below, each route is represented by a different color. Most route have common sections.

What I need is to highlight the number of overlayed route in different sections, with one color whose lightness varies.

I thought about representing the layer in red line, then apply a 80% alpha, hoping overlaying lines become more opaque. But it doesn't work in ArcMap.

Any suggestion?

I don't have Spatial Analyst extension.


texts: Error in executing Voronoi polygon algorithm

texts:
I am trying to  create Voronoi polygon on point layer of my project. I need to save this new layer as .shp file.

The following is the code I have written in Python console, using QGIS 3.0

import processing
Parameters={ 'INPUT' : 'file:///C:/Test/points.csv?delimiter=,&crs=epsg:4326&xField=Lat&yField=Long', 'BUFFER' : 0, 'OUTPUT' : 'file:///C:/Test/vor.shp' }
processing.run('qgis:voronoipolygons',Parameters)

But I am getting the following error:

Traceback (most recent call last):
File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
exec(code, self.locals)
File "<input>", line 1, in <module>
File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python/plugins\processing\tools\general.py", line 84, in run
return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python/plugins\processing\core\Processing.py", line 190, in runAlgorithm
raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.


texts: Why is result from gp tool not being created on JavaScript API?

texts:
Ok I am really confused. I can see the output of the result from geoprocessing tool. But why am I not able to add it to the map. As you can see the result returns all the attribute values.

This is my application and you will see the job id in the console. If you are not able to view the link I have added for result change the job id with the new one (in the console) that you ran and you will see the results.

What could be the reason? I am using ArcServer 10.


texts:
I have a GeoJSON file of the country of Haiti with the administrative boundaries for communes (i.e. counties). Here is the file:

https://rawgit.com/mattbowlby/bf5243ab097e50eb8864/raw/46818c96bc49629681238d5302d597456227c040/haiti_communes.geo.json

That file has a few properties, for example:

( "properties": {
        "id_com": 112,
        "commune": "Delmas",
        "id_dept": 1,
        "dept": "West Department",
        "shape_leng": 0.242799419022,
        "latitude": 18.5543057016,
        "longitude": -72.297644435,
        "shape_le_1": 0.242799419022,
        "shape_area": 0.00222841371629

Now, I want to be able to add a number of new properties, such as area, population and other things. I have 141 different communes and do not want to enter this information manually. Is there a way to work with a GeoJSON file and add the properties, such as with a spreadsheet or other format? Ideally I would have it in a spreadsheet and add in the new properties as columns.


texts: 64Bit Applications installed with QGIS/osgeo4w

texts:
Is python 64bit and gdal 64bit automatically installed when QGIS 64bit is installed via osgeo4w 64bit?

I currently have a custom 64bit python and gdal setup on my pc. However, this setup causes some issues whenever I want to import OTB, grass or qgis libraries into a py script.


texts: Attribute Assistant - Incorrect Increment Value

texts:
I'm using Attribute Assistant to generate an ID using the GENERATE_ID method.
Instead of incrementing by 1, the value increments by 2 every time.
The log file shows a curious thing:
Field Name: STRUCT_ID was found at index: 2
            Trying: GENERATE_ID
            14 is the existing value and the interval is 1: 5:11:23 PM
            SEQCOUNTER changed to 15: 5:11:23 PM
            SEQCOUNTER changed to 16
            STRUCT_ID set to SWMC_00016
            Finished: GENERATE_ID

Thus it appears that AA recognizes the proper start value, realizes it needs to increment by 1, but then proceeds to increment by 2!
Here are the relevant entries and fields in the two AA tables. There is only one GENERATE_ID rule for this TABLENAME, so that can't be an issue.

Any ideas what is happening?


texts: I try to export a QGIS 3 model (.model3) to Python but the functionality that was available in QGIS 2 doesn't exist anymore. Is there a way to export it ?


texts: How can I expand the zoom levels for a feature type in an .mbtiles file?

texts:
I have a 51 GiB mbtiles file containing vector tiles for the entire world.

I would like a particular type of feature (railways) to be available at more distant zoom levels. At present, this feature is missing from zoom levels wider than 8 (i.e. 1–7).

It seems to me that the necessary information should be already contained within the file (since the closer zoom level information can be used to produce the more distant zoom level information).

Is there some existing tool which can modify the mbtiles file accordingly?


texts: How to rename a field in a shapefile in QGIS?

texts:
I'm attempting to rename a field in a shapefile. I select the layer, enabling editing and then apply changes, and save them. When I open the file again, it's been renamed 'field_to__1'. However this only seems to happen for certain (joined?) fields.

How can I save the field name changes in QGIS?


texts:
I am a beginner with QGIS. I downloaded from different websites administrative levels for a country (0, 1, 2), full package: dbf, shp, shp.xlm, snb, shx, sbx, sbn, prj etc. However, QGIS is not opening the shapefiles (see picture attached)

Actually, only the KMZ files are working.

Could someone help me on this?


texts: How to use "this" Levenshtein Python code into ArcGIS 10

texts:
I need some help with this:

First of all, I wanted to "Compare strings from two fields", and asking before at gis.stackexchange, in this post someone help me with this topic, referring me to a Python code, but, unfortunately, from the deepest ignorance of this kind of issues, I wanted to know if it is possible to run this code within ArcGIS 10 (field calculator).

The code: PyLevenshtein ,

Thanks in advance, and sorry if I wasted your time.


texts: Merging all vector files within directory with QGIS?

texts:
I am unsuccessfully trying to merge all shapefiles in a given directory using QGIS 2.18. I do not want to navigate to the folder using QGIS GUI tools. Instead, I want my script to include all shapefiles in a certain location.
For now, what I am getting, is a copy of one random file only. I suppose the semantics of my qgis:mergevecorlayers algorithm is not appropriate.

import glob, os, processing
path_to_shp = "C:/Temp/Tracks/"
os.chdir(path_to_shp)
for file in glob.glob("*.shp"):
    uri = path_to_shp + fname
    processing.runalg("qgis:mergevectorlayers", uri, "C:/Temp/Tracks/merged.shp")


texts:
I am trying to figure out how to change OpenStreetMap country borders color/thickness.

I've set up a tile server using OSM data, PostgreSQL & Mapnik.

The default mapnik style is used.

There doesn't seem to be a straight forward way to do it globally.

I'm a newbie.


texts: Create polygon after excluding extents from raster data having no-data values

texts:
I want to get extent of raster dataset after excluding the extents containing no-data values in the raster data which I then want to export as a shapefile. Till now I have used the following code:
import rasterio

with rasterio.open(r'C:\Users\Desktop\PRODUCT1\Landsat\T1.tif') as src:
    bounds = src.bounds

from shapely.geometry import box
geom = box(*bounds)

import geopandas as gpd
df = gpd.GeoDataFrame({"id":1,"geometry":[geom]})
df.crs = src.crs
df.to_file(r"C:\Users\Desktop\PRODUCT1\Landsat\T1_boundary.shp")

After running the above code I am getting raster extents containing no-data values.
Can someone please help me out in getting the raster extents excluding the extents having no-data values.


texts: ArcGIS - Model builder (loop through many features with same name)

texts:
I want to use modelbuilder to apply symbology from layer file. I use "Apply Symbology From Layer" with one feature layer and one symbology layer as input.
In my map I often have many features with the same name, and therefore I want to loop through all the features and apply the same symbology, see attached picture.
The different input features are store in each separate geodatabase
But I can't figure out how to make this loop input function. Any tips?

EDIT 1: As response to first answer (new pictures):
ModelBuilder:

To the left: TableOfContents (symbology). To the right:Run log
(The forth loop deletes symbology established in third loop.)

EDIT 2: I tried to make a loop for each lyr-file, but it seems that I am only allowed to only add one iterator. Any tips how I can loop through the map, and link each layer in TOC to its corresponding lyr-file? Now I use "KpSamferdselPunkt" as wildcard in "Iterate Layers" to find the right layer in the first loop.
My idea was to make an iterator to each layer-type, see fig. below:


texts:
I am trying to add a layer to an ESRI ArcGIS Javascript API map using a WFS based on geoserver. The back end is pulling from POSTGIS, if it makes a difference. The WFS appears to work fine in openlayers, but not in the esri javascript api. Looking at the http requests, they are totally different. Specifically the GetFeatures Request. In openlayers  this is the request:

  [URL]:8080/geoserver/wfs?service=WFS&version=1.1.0&request=GetFeature&typename=conryservices:zip_polygons&outputFormat=text/javascript&format_options=callback:loadFeatures&srsname=EPSG:3857&bbox=-8913168.994277833,5390950.7308969125,-8903385.05465733,5400734.670517415,EPSG:3857&_=1469810854811

With a nice json response. ESRI requires me to use their proxy, but even still , the request is a POST and looks like this:

  [URL]/sproxy/proxy.php?[URL]:8080/geoserver/wfs

<?xml version='1.0' encoding='utf-8'?>
<GetFeature
 xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance'
 xmlns:gml='http://www.opengis.net/gml'
 xmlns:ogc='http://www.opengis.net/ogc'
 xmlns:wfs='http://www.opengis.net/wfs'
 xmlns='http://www.opengis.net/wfs'
 xmlns:conryservices='[URL]'
 version='1.1.0' service='WFS' maxFeatures='100'>
 <wfs:Query typeName='conryservices:zip_polygons' srsName='EPSG:4326'>
  <ogc:Filter xmlns:ogc='http://www.opengis.net/ogc'>
   <ogc:BBOX>
    <ogc:PropertyName>conryservices:geom</ogc:PropertyName>
    <gml:Box srsName='EPSG:4326'>
      <gml:coordinates>-122.94620660400169,42.29605553470161 -122.72579339599417,42.395904817819165</gml:coordinates>
    </gml:Box>
   </ogc:BBOX>
  </ogc:Filter>
 </wfs:Query>
</GetFeature>

I'm not sure what I'm missing here, but are these running on the same standard?


texts: ArcScene 3D vrml export / rhino import error

texts:
Arcscene does not give an error when exporting 3D (vrml) but when trying to import into Rhino, an error "" shows up. Importing into AutoCad and 3dsmax also fails. The ArcSene version is 10.5.1. Screenshots of the vrml are attached.
Any ideas why this is happening?


texts: Changing Names in TOC Doesn't Update Legend Properly

texts:
This is related to an ArcMap add-in I coded up this past week.

I'm trying to rename a layer in ArcMap and have the legend refresh properly. The renaming part is working fine, however I cannot get the legend in the layout to update properly.

I call CurrentContents.Refresh() to refresh the name in the table of contents, which works.

If I search through the map surrounds and call .Refresh() on the legend, then do an ActiveView.PartialRefresh() on it, the legend gets all out of whack. This is why I am trying to simply trigger the same event as when the user manually changes the name of a layer.

The API documentation seems to say that MxDocument.UpdateContents() will trigger the events needed to refresh the legend. However, this does not seem to be the case.
In addition to the above methods I have tried ActiveView.ContentsChanged() to no avail.

Edit: I am writing this in C# but can read VB if you feel more comfortable writing in that.

Any suggestions are welcome.

Thanks,

Mat


texts: Creating raster map using Excel file and calculating averages?

texts:
I have thought about using ArcGIS to create a raster map of latitude/longitude (0 to 360 degrees) temperature data (in Excel format). This temperature data reflects the total amount of warming that has occurred per degree Celsius of global warming by location. My data collection process involves me upscaling regional temperature warming into a global mean temperature warming.

What I am attempting to do is transform the provided Excel data into a raster map in ArcGIS. In doing this, I would be able to visually see regional temperature increases per degree Celsius of global warming. The idea is to be able to more flexibly calculate regional temperature averages per degree Celsius of global warming (such as in the NW Pacific, or the North Atlantic), so that I can upscale these averages into a global mean temperature warming.

Would it be possible to calculate these averages in this way? If so, what procedure can I follow?


texts: Editing polyline data, can't snap or modify feature

texts: Editing in Desktop standard 10.3 and I am having issues with snapping polylines. I keep getting the same error below when I try to snap or move a line. Any ideas On why this might be happening?


texts: Changing Excel file to shapefile for use in QGIS?

texts:
I have my data in an Excel sheet.

How do i convert it to CSV and shapefile for use in QGIS?


texts:
QGIS Zurich 3.18.2 Need to pass on a Shapefile to A Nother. There were three Layers each with assigned Label Style (saved) as can be seen
Yes I know I have a lot more work to do on the Label placement/conflicts.
I then successfully merged these three shapefiles (adding a Unique Field Name for each original Layer) and have tried to setup a Rule Based Labels for each of the selected Unique attribute data, selecting then Load Style for that data BUT each time it reverts to a Single Label style

Am I going about it wrong or is it not possible?


texts: Selecting features through ArcPy gives Assertion Error?

texts:
I would like to achieve the following functionality through arcpy:

Querying an existing baselayer on time and create a new layer
Apply symbology from an existing layer and add the new layer to the current session's data frame on top of the basemap

Have given my code below which gives me an Assertion Error when I try to add the layer returned from arcpy.MakeFeature. I am missing something trying to impact the current session. If I run the arcpy.MakeFeatureLayer_management command on the python console it works fine. But I am not able to get it work through the scripts in the toolbox. Please advice me on this.

import arcpy
arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
arcpy.MakeFeatureLayer_management ("base_layer", "new_layer_1","DATE_TIME" = TO_DATE('13-JUL-2010 12:40:00','DD-MON-YYYY HH24:MI:SS'))
arcpy.ApplySymbologyFromLayer_management("out_time",r"C:\WIP\Shape files\KARTHIK.TEST_TABLE.lyr")
arcpy.mapping.AddLayer(df,out_time)
arcpy.RefreshActiveView()


texts: Calculating individual width between buildings using ArcMap?

texts:
I have a building footprint shapefile and I would like to calculate the individual width between each building (i.e., adjacent building and across the street). I have looked at solutions which use the road network but they turned out to be unfeasible as some areas do not have roads.

Is there any way to use the Near tool or other tools to achieve this?

Using the Euclidian distance tool with the building footprint as input, I have generated a raster with a cell size of 1. I am now hoping to find a more accurate and less visual solution to generate data on gaps between buildings.


texts: What event is fired when a feature is edited?

texts:
When the user draws a new feature (line, point, etc.) in ArcMap, there are null attributes.
If the user enters values in these null fields what is the event that is fired?

I want to control the entered values.


texts: Explaining pyproj to_epsg min_confidence parameter?

texts:
Following on from this question How can I get epsg code in pyproj with version==2.1.3, I discovered that sometimes when intialising a pyproj Proj instance with a string ESPG code eventually results in an ambiguity when trying to determine the projection's EPSG code.

from pyproj import Proj
crs = Proj(init='epsg:4544').crs
print(crs.to_epsg())
# prints None

If you adjust themin_confidence parameter, then you can get it out:

from pyproj import Proj
crs = Proj(init='epsg:4544').crs
print(crs.to_epsg(min_confidence=25) # <-- !!
# prints 4544

I feel like this is something I should understand, but I don't. What is the ambiguity here; why does the to_epsg method require this min_confidence parameter to return a valid EPSG code in this case?

Reading the source code is not enlightening regarding the source of the ambiguity:

Parameters
----------
min_confidence: int, optional
  A value between 0-100 where 100 is the most confident. Default is 70.

Returns
-------
int or None: The best matching EPSG code matching the confidence level.

Reading that function, as well as to_authority, I surmise that the ambiguity arises when there are multiple matches, and the confidence score allows pyproj to return the best candidate. I don't particularly understand that either, but it does make sense. What really doesn't make sense is how None can be a better candidate than 4544 in this case, especially when the Proj instance was initialised without an exception or warning.

Notably, using a CRS instance instead, but with the same EPSG string input, seems to be unambiguous:

from pyproj import CRS
crs = CRS('epsg:4544')
print(crs.to_epsg())
# prints "4544"

>>> import pyproj
>>> pyproj.__version__
'2.2.0'


texts:
I'm using Rasterio to read GeoTIFF files from Landsat 8 and calculate NDVI into a new GeoTIFF file.

My code looks like this:

import numpy
import rasterio
import subprocess

with rasterio.drivers(CPL_DEBUG=True):

    # Read raster bands directly to Numpy arrays.
    #
    dsRed = rasterio.open('downloads/LC81970212013122LGN01/LC81970212013122LGN01_B3.TIF')
    bandRed = dsRed.read_band(1)

    dsNir = rasterio.open('downloads/LC81970212013122LGN01/LC81970212013122LGN01_B5.TIF')
    bandNir = dsNir.read_band(1)

    ndvi = numpy.zeros(dsRed.shape, dtype=rasterio.uint16)

    ndvi_upper = bandNir + bandRed
    ndvi_lower = bandNir - bandRed

    ndvi = ndvi_lower / ndvi_upper

    kwargs = dsRed.meta
    kwargs.update(
        dtype=rasterio.uint16,
        count=1,
        compress='lzw')

    with rasterio.open('example-total.tif', 'w', **kwargs) as dst:
        dst.write_band(1, ndvi.astype(rasterio.uint16))

# Dump out gdalinfo's report card and open the image.
info = subprocess.check_output(['gdalinfo', '-stats', 'example-total.tif'])
print(info)

subprocess.call(['open', 'example-total.tif'])

But it isn't really producing the result I was hoping for.

I'm probably getting the data types wrong, I'm not too familiar with Python and NumPy.

The resulting image is all black. If I multiple the result array by e.g. 1000 it becomes visible, but it isn't the correct values.


texts: Can't enable Jsonp in Geoserver 2.7.0 and then publishing my WFS request in Leaflet

texts:
I'm working on localhost:8080 and want to get WFS data in a Leaflet map.
So, I've read in some posts here that I first I have to enable Jsonp in my web.xml file in Geoserver.
So I put

<context-param>
<param-name>ENABLE_JSONP</param-name>
<param-value>true</param-value>
</context-param>

But I can't see that Jsonp is allowed with a getCapabilities :

<ows:Parameter name="outputFormat">
<ows:AllowedValues>
<ows:Value>text/xml; subtype=gml/3.2</ows:Value>
<ows:Value>GML2</ows:Value>
<ows:Value>KML</ows:Value>
<ows:Value>SHAPE-ZIP</ows:Value>
<ows:Value>application/gml+xml; version=3.2</ows:Value>
<ows:Value>application/json</ows:Value>
<ows:Value>application/vnd.google-earth.kml xml</ows:Value>
<ows:Value>application/vnd.google-earth.kml+xml</ows:Value>
<ows:Value>application/x-gpkg</ows:Value>
<ows:Value>csv</ows:Value>
<ows:Value>geopackage</ows:Value>
<ows:Value>geopkg</ows:Value>
<ows:Value>gml3</ows:Value>
<ows:Value>gml32</ows:Value>
<ows:Value>gpkg</ows:Value>
<ows:Value>json</ows:Value>
<ows:Value>text/xml; subtype=gml/2.1.2</ows:Value>
<ows:Value>text/xml; subtype=gml/3.1.1</ows:Value>
</ows:AllowedValues>

Then, I tried to load my WFS data into a GeoJson layer in leaflet with this code :

<!DOCTYPE html>

<meta http-equiv="Access-Control-Allow-Origin" content="*">
<meta charset="utf-8" />
<title>Test WFS</title>

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

And firefox's console return

SyntaxError: missing ; before statement    ows:1:7

Is someone knows why it doesn't work ?


texts: Filling NA values of raster stack using R?

texts:
I am trying to fill NA values of a raster stack with interpNA but the first raster stays empty. I wish to obtain similar results as in substituteNA in a sense that all rasters of the stack are filled.

How can I get around this?

#rasters
set.seed(123)
library(raster)
library(timeSeries)

r1 <- raster(nrows = 1, ncols = 1, res = 0.5, xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, vals = 0.3)
rr <- lapply(1:10, function(i) setValues(r1,seq(1,ncell(r1),1)))
for (i in 2:length(rr)){
  rr[[i]] <- rr[[1]]+i
}

#make and stack with NA
s <- stack(rr)
pb = txtProgressBar(min = 0, max = nlayers(s), initial = 0, style = 3)
for (i in  seq(1,nlayers(s),2)){
  values(s[[i]])[values(s[[i]])%% 2 == 0] <-  NA
  setTxtProgressBar(pb,i)
}
plot(s)

##First raster will be empty:
fun <- function(x) {
  v=as.vector(x)
  z=interpNA(v, type="linear") #substituteNA(v, type="mean")
}
s1<-calc(s, fun)
plot(s1)


texts: SQL geography::EnvelopeAggregate output not accurate

texts:
I am trying to get the geography::EnvelopeAggregate for a set of locations. The geography column values are calculated using the Lat/Long for individual points. When I try to fetch the boundary box(Extent) the box does not cover all the individual points. Is this expected behaviour? Can someone who has expertise in this field provide some inputs. Below is the SQL code and corresponding Map plotted. SQL2014 is being used in this example.

IF EXISTS (SELECT * FROM tempdb.dbo.sysobjects WHERE id = OBJECT_ID(N'tempdb..[#ServicePoints]'))
    BEGIN
        DROP TABLE #Points
    END
CREATE TABLE #Points
(
    PointId nvarchar(max),
    Lat float,
    Long float,
    Location geography
)
insert into #Points(PointId, Lat, Long, Location)
values ('point1', 47.659226, -117.218292, geography::Point(47.659226, -117.218292, 4326))
insert into #Points(PointId, Lat, Long, Location)
values ('point2', 47.659226, -117.218292, geography::Point(47.659226, -117.218292, 4326))
insert into #Points(PointId, Lat, Long, Location)
values ('point3', 47.659226, -117.218292, geography::Point(47.659226, -117.218292, 4326))
insert into #Points(PointId, Lat, Long, Location)
values ('point3', 47.658276, -117.218474, geography::Point(47.658276, -117.218474, 4326))
insert into #Points(PointId, Lat, Long, Location)
values ('point3', 47.658323, -117.21899, geography::Point(47.658323, -117.21899, 4326))
insert into #Points(PointId, Lat, Long, Location)
values ('point3', 47.658341, -117.219487, geography::Point(47.658341, -117.219487, 4326))
insert into #Points(PointId, Lat, Long, Location)
values ('point3', 47.660115, -117.219365, geography::Point(47.660115, -117.219365, 4326))
insert into #Points(PointId, Lat, Long, Location)
values ('point3', 47.658478, -117.216265, geography::Point(47.658478, -117.216265, 4326))

select * from #Points
--select geography::Point(47.658323, -117.21899, 4326)
Select E1.ID,
                   E1.SPEXtents.STPointN(3).Lat as BottomLeft_Lat,
                   E1.SPEXtents.STPointN(3).Long as BottomLeft_Long,
                   E1.SPEXtents.STPointN(1).Lat as RightTop_Lat,
                   E1.SPEXtents.STPointN(1).Long as RightTop_Long,
                   E1.SPExtents.STPointN(2).Lat as BottomRight_Lat,
                   E1.SPExtents.STPointN(2).Long as BottomRight_Long,
                   E1.SPExtents.STPointN(4).Lat as TopRight_Lat,
                   E1.SPExtents.STPointN(4).Long as TopRight_Long

        from
        (SELECT 1 as ID,
                geography::EnvelopeAggregate(sp.Location) AS SPExtents
        FROM #Points sp
        ) E1;

And the representation of these points in map is below. As can be seen 2 points are lying outside the boundary.

Or Should I be plotting a circle in which case all points will fall within the circle as shown below.


texts: Sharing/exporting map created in QGIS for different people in Company?

texts:
I created a QGIS map with several different layers.
I want to know a light way method to share it locally with staff. If I can be pointed to the different export formats that may work for me in Qgis. I can experiment and see what may work
The map has to be able to display certain attribute features when a certain layer is clicked.
I exported it as a kml file. But it is very large and slow.


texts: Query Layer to access data on Sharepoint

texts: I want to connect to a Sharepoint list from ArcMap.  Without knowing very much about these things I thought that it would be possible to read this non-spatial data to create lookup tables in my file geodatabase.  In searching for an answer, I can't seem to get away from ArcGIS for Sharepoint, which isn't what I'm looking for at all.  I'm using ArcGIS 10.0


texts:
I've been using insert, update and select cursors found in arcpy.da and so far it's met my needs.

Today, I found another situation where I really just want to throw a SQL statement at the database, from my python/arcpy code, and get the resultset. Specifically, I want to get the Id of a polygon that contains a point; very easy to do with a SQL statement in my SQL tool.

There have been other examples over the last few weeks so this question is not about how do perform the point / polygon intersection. It's about running SQL statements from within my ArcPy code. In other parts of my code, ArcPy uses ArcSDE connection files to access the ArcSDE feature classes and it would be great if I could use the existing data access infrastructure and just send a SQL statement.


texts: Creating polygons (squares) from point using PostGIS?

texts: I want to create a regular grid of polygons/squares from this point knowing that between two successive points 100 meters, so each square would have 100*100 m2 in its area.


texts: GeoServer Maven Quickstart fails during build process

texts:
I am trying to build GeoServer using Maven. I followed the instructions from the GeoServer  website Instruction from GeoServer site. When i get to the build stage it fails. I have tried to search for cases where this happened but not been able to find any solutions.
In the error reads Invalid escape character at line 1...... attached is a screen shot of the error.

I am using Windows 10, apache-maven3.6.0 bin. Any ideas on how to solve this?


texts: ArcMap trace the creation/modification/deletion events by windows user for all changes (not only the latest change)

texts:
I want to check whether there is a way to trace all the creation, modification and deletion events in the editing session. What I aware now is how to get the latest updated information of an object in the attribute table of a layer.

I am using ArcMap 10.2.2.


texts: Finding Nearest Feature Methods VB ArcObjects

texts:
I have:

(a) point feature class
(b) polyline feature class

I am trying to loop through the (a) point feature class to find the nearest feature from the (b) polyline feature class.

Just as a start to more complex command, I'd like to find and select the nearest (b) polyline, and create a message box that tells me its ObjectID, just so I know that the method for finding the nearest line is working.

I've come across two options, both of which I'm having issues with.

(1) The IIndexQuery.NearestFeature method, which I could use in combination with the IFeature.GetFeature method to select the polyline feature of interest.

I've successfully entered the parameters (I think), but the NearestFeature method returns a "-1" for each of the features FIDs. Maybe this is because my features have an 'ObjectID' field, but no 'FID' field, though, I tried the same code with a polyline feature that did have an FID field, and it still gave me -1's across the board. Could anybody speak to this?

(2) The IHitTest.HitTest method, which returns a segment index of the nearest feature.

The problem here is that I'm looking for a value in a field for the whole line, not just a segment of the line. Is there a way to return the whole feature from the segment index?

Any other suggestions on finding the nearest line from a given point?

I couldn't figure out how to get the coloring right on the code below...

Private Sub NHD_Prep_Click()

'Set general variables

Dim pMxDoc As IMxDocument
Set pMxDoc = ThisDocument
Dim pMap As IMap
Set pMap = pMxDoc.FocusMap

'Set variables for the line layer of interest

Dim pNhdLayer As ILayer
Set pNhdLayer = pMap.Layer(0)
Dim pNhdFLayer As IFeatureLayer
Set pNhdFLayer = pNhdLayer
Dim pNhdFClass As IFeatureClass
Set pNhdFClass = pNhdFLayer.FeatureClass
Dim pNhdFeatureIndex As IFeatureIndex
Set pNhdFeatureIndex = New FeatureIndex
Set pNhdFeatureIndex.FeatureClass = pNhdFClass
Dim pNhdIndexQuery As IIndexQuery2
Set pNhdIndexQuery = pNhdFeatureIndex

'Set variables for the point layer of interest

Dim pLayer As ILayer
Set pLayer = pMap.Layer(1)
Dim pFLayer As IFeatureLayer
Set pFLayer = pLayer
Dim pFClass As IFeatureClass
Set pFClass = pFLayer.FeatureClass

'Set number variables

Dim nearReachFID As Long
Dim distance As Double
Dim strReach As String
Dim dblmeas As Double

'Check to make sure the layer of interest is on top

If pMxDoc.FocusMap.Layer(0).Name <> "NHD Flowline" Then
    MsgBox "Please place NHD Flowline on the top spot of the TOC..."
    Exit Sub
End If

'Start Cursor on point layer features

Dim pFCursor As IFeatureCursor
Set pFCursor = pFClass.Search(Nothing, True)
Dim pFeature As IFeature
Set pFeature = pFCursor.NextFeature

Do Until pFeature Is Nothing

    'Use nearest point method to find the line FID and distance to it

    Dim pPointFeature As IPoint
    Set pPointFeature = pFeature.Shape

    pNhdIndexQuery.NearestFeature pPointFeature, nearReachFID, distance

    MsgBox "This is the " & pFeature.Value(pFeature.Fields.FindField("Name")) & " feature. " _
    & "The distance from this point to the nearest line, " & nearReachFID & ", is " & distance

    Set pFeature = pFCursor.NextFeature
Loop

End Sub


texts: Displaying bounding box array on leaflet map

texts:
I have a set of bounding box arrays which I want to display over leaflet map in my application. The arrays look something like this with the bounding box co-ordinates in it:

Array[4]
0:
-2.362060546875
1:
51.96245837645124
2:
7.62162297964096
3:
51.96245837645124

The number of arrays can increase depending on the searching done by user. I want to have a functionality where it goes through each array and displays the bounding box on map. Any suggestions?

Here is how I am fetching the data and storing bounding box in separate array

    function map(obj){

               var o = obj.data.hits.hits;
               var b=[];
                for(var i in o){
                    o[i]._source.id = o[i]._source.compendium_id;
                    try{
                    b.push(o[i]._source.metadata.o2r.spatial.union.geojson.bbox);
                    throw error()
                  }
                  catch (g){
                    console.error("missing spatial");
                    }
                }
                $log.debug('mconfiging result: %o', o,b);
                return obj;
            }
here is my map creation

angular.extend($scope, {
                center: {
                    lat: 51.505,
                    lng: 10.09,
                    zoom: 3
                },
                controls: {
                  scale:true,
                  draw: {}

                },

                layers: {
                    baselayers: {
                        mapbox_light: {
                            name: 'Mapbox Streets',
                            url: 'https://api.mapbox.com/styles/v1/mapbox/streets-v10/tiles/256/{z}/{x}/{y}?access_token=pk.eyJ1IjoicmVoYW5zNTE2IiwiYSI6ImNpeWxjcWNkODAwNGwzM3FxamR6a2gxOXkifQ.PuUfs90MyfmVGYVqx0AoUw',
                            type: 'xyz',
                            layerOptions: {
                                apikey: 'pk.eyJ1IjoicmVoYW5zNTE2IiwiYSI6ImNpeWxjcWNkODAwNGwzM3FxamR6a2gxOXkifQ.PuUfs90MyfmVGYVqx0AoUw',
                                mapid: 'mapbox.streets',
                                format: '@2x.png'
                            },
                            layerParams: {
                                showOnSelector: false
                            }
                        }
                    },
                    overlays: {
                        draw: {
                            name: 'draw',
                            type: 'group',
                            visible: true,
                            layerParams: {
                                showOnSelector: false
                            }
                        }

                    }
                }

           });

           leafletData.getMap().then(function(map) {

               leafletData.getLayers().then(function(baselayers) {
                 var drawnItems = baselayers.overlays.draw;
                  map.on('draw:created', function (e) {

                    var layer = e.layer;

                    drawnItems.addLayer(layer);

                    console.log(JSON.stringify(layer.toGeoJSON()));
                  });
               });
           });


texts: Why does geotools DataStoreFinder return null jdbc dataStore

texts:
I want to create a JDBC DataStore with these parameters:

Map<String,Object> params = new HashMap();
params.put("dbtype","postgis");
params.put("host","localhost");
params.put("port",5432);
params.put("schema","public");
params.put("database","mydb");
params.put("user","postgres");
params.put("passwd","postgres");

DataStore store = DataStoreFinder.getDataStore(params);
System.out.println(store);// prints null

Iterator it = DataStoreFinder.getAvailableDataStores();
while(it.hasNext()){
  System.out.println(it.next());
}

the output:

null
org.geotools.data.shapefile.ShapefileDirectorFactory@5e316c74
org.geotools.data.shapefile.ShapefileDataStoreFactory@7ee955a8

my libraries (pom.xml):

<dependency>
  <groupId>org.geotools</groupId>
  <artifactId>gt-main</artifactId>
  <version>19-SNAPSHOT</version>
</dependency>
<dependency>
  <groupId>org.geotools</groupId>
  <artifactId>gt-jdbc</artifactId>
  <version>19-SNAPSHOT</version>
</dependency>
<dependency>
  <groupId>org.geotools</groupId>
  <artifactId>gt-data</artifactId>
  <version>19-SNAPSHOT</version>
</dependency>
<dependency>
  <groupId>org.geotools</groupId>
  <artifactId>gt-shapefile</artifactId>
  <version>19-SNAPSHOT</version>
</dependency>

Why does it create a null datastore?


texts: ArcGIS Pro select-by-location using intersect not working as anticipated

texts:
I'm using ArcGIS Pro and I have a polygon layer that needs to be selected if it intersects with another polygon layer in "red", as seen in the picture. However when I use the select-by-location tool and use intersect it selects the entire layer and not just the polygons that intersect.
I've tried using repair geometry, ensured my CRS is matching, and ensured they are not multi-polygons but I haven't figured out what is causing this. Any suggestions?


texts: Save an 8 band raster in one TIFF rasterio

texts:
I'm trying to write an 8 band image in just one TIFF using rasterio.
I actually opened an 8 band raster, normalize it for further processing and tried to re-write the normalized bands: I used
with rasterio.open('CONCAT_NORM.tif','w',**profile) as outds:
     for b in range(b_array.shape [-1] ):
         outds.write(b_array[:,:,b].astype(rasterio.float32))

with b-array is the normalized array read from the original image in shape (3214, 5478, 8)
However this error appears:
ValueError: Source shape (3214, 5478) is inconsistent with given indexes 8

I tried to transpose the array's columns from (0,1,2) to (2,0,1) but the result was totally changed.


texts: Exporting merged layers so coordinates are in the same CRS using QGIS

texts:
I have two different data point sets, one with Easting/Northings, one with Latitude/Longitude. I can display them on one map and also merge the layers in QGIS 3 (3.4.4) just fine. However, when I export the merged layer as CSV-file, the coordinates are still either East/North or Lat/Long.
How do I get all of them in the same format (East/North or all Lat/Long)?


texts:
I am a beginner in QGIS (v3). I have a point layer with a lot of overlapping and I made diagrams for each point. I would like to displace them as the grid, just like the usual point displacement renderer.

Do you think it is possible? If so how? I couldn't find any answers in the documentation.


texts: Icons of different layers cover each other in ArcMap?

texts:
I have converted 6 KML files (including a set of points) into 6 layers by ArcMap 10. For points in each layer, I have considered a special icon with a specific color. I also decreased the size of icons as much as possible. But the icons of different layers cover each other and the picture is not desirable enough to be inserted in an article.

Is there any way to convert the icons into transparent ones?


texts: How to use GDALGridInverseDistanceToAPower in a python script

texts:
I'm trying to use GDAL in a python script to do an inverse distance interpolation. Within the C documentation API I found GDALGridInverseDistanceToAPower function but I don't know how/if I can call that from Python.

Does somebody know how to perform this analysis or call this GDAL function in python?


texts: OpenLayers with Google, Bing & MapQuest - Does it violate the terms of service?

texts:
to my understanding, you are not allowed to use the Google Maps API to display Google and Bing maps (Layer Toggle) (From the terms of service). However, I'm unclear if I can have an OpenLayers map that allows the users to switch between Google, Bing and MapQuest. I want to let the user be able to choose the best Satellite/Aerial data available for their location. I'm really unclear if this is a Terms Of Service violation.

Any ideas?

Thanks,

Nick,


texts: Time-slider in ArcGIS skipping once exported to video

texts:
The set-up: I am using ArcGIS 10.3. I have exactly 528 raster maps in a mosaic dataset. Each one corresponds to the same geography but at different months (from January 1969 to December 2012). I would like to use time-slider to move through these images and export to a video.

The time-slider animates nicely through the rasters inside ArcGIS, but once I try to export it to an .avi video it skips the large majority of images and I end up with only about every sixth raster being displayed.

Potentially relevant: the time-slider is operating on a string field that is the concatenated year and month (YYYYMM format). I created this field using the field calculator. It looks like this:

I picked this format because it fits one of the options in the Layer Properties dialogue. Note that in this dialogue I do not set the layer time extent because if I have ArcGIS calculate it directly, it decides the time starts at 1/1/0001 for some reason (this is highlighted in the next picture)

So in the time-slider properties dialogue, I change 1/1/0001 to 1/1/1969 and it seems to work fine. From here, the animation works in ArcGIS.

Using the default settings (below) to export to video causes the skipping. When I tested the default settings on a smaller subset (24 months) it worked fine. It's hard to test settings on the whole thing since the video usually takes a long time to export (and often doesn't work if I open any windows over the top of it).


texts: How to set label font color using data defined property for PyQGIS 3?

texts:
Set label font color using data defined property in PyQGIS 2 is well established, e.g. in another thread

However, PyQGIS 3 completely changed the way this works with no well documented examples. So, how to setup labels whose colors are based on data defined property. Here is my initial try, the label text show up, but the label's font color is not changing as i specified to be dependent upon the 'color' field values:

prop=QgsProperty()
prop.setField("Color")
pc=QgsPropertyCollection('mycol')
pc.setProperty(0,prop)
pal_layer=QgsPalLayerSettings()
pal_layer.setDataDefinedProperties(pc)
pal_layer.fieldName="value"
pal_layer.enabled=True
labeler=QgsVectorLayerSimpleLabeling(pal_layer)
layer.setLabeling(labeler)
layer.setLabelsEnabled(True)


texts: Convert a GeoJSON file coming from Overpass in shapefile using ogr2ogr in command line

texts:
I need to convert a GeoJSON file coming from an Overpass query in a ESRI Shapefile using ogr2ogr in command line (I'm on Ubuntu 15.10 .. )

This is my query edited in Overpass Turbo

[out:json][timeout:250];
{{geocodeArea:Italia}}->.searchArea;
(
  node["amenity"="hospital"]["emergency"="yes"](area.searchArea);
);
out body;
>;
out skel qt;

and how is converted in Overpass XML (using Export functionality in Overpass Turbo .. )

<osm-script output="json" output-config="" timeout="250">
  <id-query into="searchArea" ref="3600365331" type="area"/>
  <union into="_">
    <query into="_" type="node">
      <has-kv k="amenity" modv="" v="hospital"/>
      <has-kv k="emergency" modv="" v="yes"/>
      <area-query from="searchArea" into="_" ref=""/>
    </query>
  </union>
  <print e="" from="_" geometry="skeleton" limit="" mode="body" n="" order="id" s="" w=""/>
  <recurse from="_" into="_" type="down"/>
  <print e="" from="_" geometry="skeleton" limit="" mode="skeleton" n="" order="quadtile" s="" w=""/>
</osm-script>

I put this in a file named ps-Italy-from-OSM.txt and then I'm download the data using this command:

wget -e use_proxy=yes -e http_proxy=<proxy>:port> --post-file=ps-Italy-from-OSM.txt http://overpass-api.de/api/interpreter --output-document=ps-Italy-OSM.json

Al works fine but, when I try to convert the json data in an ESRI Shapefile using this ogr2ogr command from command line

ogr2ogr -nlt POINT -skipfailures ps-Italy-OSM-WGS84.shp ps-Italy-OSM.json OGRGeoJSON

I obtain this error ...

ERROR 4: Failed to read GeoJSON data
FAILURE:
Unable to open datasource `ps-Italy-OSM.json' with the following drivers.
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> SDTS
  -> TIGER
  -> S57
  -> DGN
  -> VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> NAS
  -> GML
  -> GPX
  -> LIBKML
  -> KML
  -> GeoJSON
  -> Interlis 1
  -> Interlis 2
  -> GMT
  -> GPKG
  -> SQLite
  -> DODS
  -> ODBC
  -> WAsP
  -> PGeo
  -> MSSQLSpatial
  -> OGDI
  -> PostgreSQL
  -> MySQL
  -> PCIDSK
  -> OpenFileGDB
  -> XPlane
  -> AVCBin
  -> AVCE00
  -> DXF
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> VFK
  -> PGDump
  -> OSM
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> PDS
  -> WFS
  -> HTF
  -> AeronavFAA
  -> Geomedia
  -> EDIGEO
  -> GFT
  -> GME
  -> SVG
  -> CouchDB
  -> Idrisi
  -> ARCGEN
  -> SEGUKOOA
  -> SEGY
  -> XLS
  -> ODS
  -> XLSX
  -> ElasticSearch
  -> PDF
  -> Walk
  -> CartoDB
  -> SXF


texts: How can models be integrated into Workflow Manager?

texts:
I have a ModelBuilder model that determines whether a field in a table contains null or duplicate values. The outputs of the model are two integer parameters, one with the count of nulls and one with the count of duplicates.

I have a workflow in Workflow Manager that I would like to integrate the model into in order to implement logic branching.

Is there a way in Workflow Manager to:

Run the model as a step in a workflow
Examine the output of the model and branch into 1 of 2 logical workflow branches.

For 1) I assume that I must run the JTXSteps.ExecuteGPTool custom step, is that correct?

For 2) (How) can I set the return code in the model? Must I refactor the model so that it only has one output parameter, or should I actually cause an error to occur for the return code to be non-zero? How can I use the return code to specify the logical branch to step into?


texts:
We've recently started working with ArcMap, having used MapInfo previously. I've created a shapefile in ArcCatalog and want to update some columns using SQL. In MapInfo I'd write the following command in the MapBasic window.

UPDATE Clearance SET GTID = "GT" & FID

I can't find where to do this in ArcMap - any pointers?


texts:
I have calculated NDVI (vegetation index) in ENVI 5.2 and I saved it in TIFF format. Now I am trying to open it in ArcMap, but its not opening properly. The image is total black and the values are showing from 9 to -999.

How do I open the NDVI file properly in ArcMap?


texts: I have some photos. They have GPS information. I want to add them with a layer like markers to OpenLayers map. How can I do that?


texts: Point Map Automation with Data Driven Pages

texts:
I need to create a map book for a large point dataset. The end user wants all points on each page numbered 1 through N and have each number correspond to a table in the layout showing one field from the attribute table. The number of points in each extent ranges from 1 to 300+ in some cases. Ideally this would be an automated process so that when I update data in the future, I can simply export and deliver a new set of maps.

So far, the only option I have found is using the Production Mapping toolbar which I don't have access to. I could get a brief trial but that would end and leave me out to dry for future map iterations. Are there any tools or processes I'm missing that can be used for this?


texts: Unable to access Geoserver on Ubuntu Server

texts:
I installed the complete opengeo-suite on Ubuntu Server and it displays that it is working fine here,

opengeo-geoserver:
  Installed: 2.5
  Candidate: 2.5
  Version table:
 *** 2.5 0
        500 http://apt.opengeo.org/ubuntu/ lucid/main amd64 Packages
        100 /var/lib/dpkg/status

But when I try to run the url (http://serverurl/geoserver/web) i am not able to get the admin page.


texts: Calculate Quickest Route Between All Features

texts:
I am looking for a way to create a line feature class that is the quickest path (as the bird flies) between all features in a feature class. I feel like there has to be an ArcGIS Tool that can do this, but my Google-fu isn't working. I think I could program something that uses the proximity analysis then manually create line vertices via the x and y of the closest point, but again there has to be some ready made tool for this isn't there? I have the latest ArcGIS Pro with all the extension.
Here is an example of what I'd like to do.
Points:

Line connecting the dots:


texts:
I want to use the command gdal_polygonize image.png bla in a Java for-loop, where image.png and bla changes in every iteration. How can I do this?

I searched for it on google and stackexchange and I only found a similiar question for Python (but it is not quite answered I think).

Any suggestions? Maybe a Batch script which I could run in every iteration?


texts: NDVI time series chart showing multiple values for the same point in Google Earth Engine

texts:
I tried to generate a time series chart of NDVI using Sentinel-2 data for multiple points. But the chart has vertical plots with two values for same point and same date.

Here is the code:

var points = ee.FeatureCollection([
ee.Feature(geometry2, {label: 'm1'}),
ee.Feature(geometry3, {label: 'm2'}),
ee.Feature(geometry4, {label: 'm3'}),
ee.Feature(geometry5, {label: 'm4'}),
ee.Feature(geometry6, {label: 'm5'}),
ee.Feature(geometry7, {label: 'm6'}),
ee.Feature(geometry8, {label: 'm7'}),
ee.Feature(geometry9, {label: 'm8'}),
ee.Feature(geometry10, {label: 'm9'}),
ee.Feature(geometry11, {label: 'm10'})
]);

Map.addLayer(points);

// NDVI calculation

var addNDVI = function(image) {
  return image.addBands(image.normalizedDifference(['B8', 'B4']));

};

  var with_ndvi = puli.map(addNDVI);

  var rgb_vis = {min: 0, max: 5000, bands: ['B4', 'B3', 'B2']};
  Map.addLayer(puli, rgb_vis, 'RGB');
  Map.addLayer(with_ndvi, {bands: 'nd', min: 0, max: 1}, 'NDVI');

   // Time Series Analysis for NDVI

var TimeSeriesndvi = ui.Chart.image.seriesByRegion(with_ndvi,points,ee.Reducer.mean(),'nd',200, 'system:time_start', 'label')
.setChartType('ScatterChart')
        .setOptions({
          title: 'NDVI variations',
          vAxis: {title: 'NDVI'},
          lineWidth: 1,
          pointSize: 4,
          series: {
          // 0: {color: 'FF0000'},
          // 1: {color: '00FF00'},
          //   2: {color: '0000FF'},
}});
print(TimeSeriesndvi);


texts: Plotting latitude and longitude coordinates (from Bing Maps) in ArcGIS

texts:
I have an excel spreadsheet containing latitude and longitude values for a set of businesses.  I want to plot the businesses in ArcMap.  I was told the latitude and longitude values were obtained from Bing Maps, which I believe uses the coordinate reference system WGS84 (I found this info on various forums, so it may be wrong!).

However when I plot the data using WGS84 (using the 'display xy data' function), the points do not appear to plot correctly. The points look to be systematically displaced west and slightly north.

Indeed, when I manually geocoded the food businesses (using the postcodes of the food businesses) this systematic displacement was evident (see image).  Does anyone know what might be going wrong?


texts: How to get polygon coordinates for Google Map from .poly files (preferably in nodejs)?


texts: OL application is at different production server than GeoServer

texts:
My OL ASP.Net web application is hosted at 172.24.100.16 production server (windows server 2008)whereas Geoserver is hosted at 172.24.100.15 production server (Linux, UBUNTU). When i access page from outside the production environment, then Layers does not show over google map in page.

I need suggestion to fix this problem. Also, i need explanation why it happen ?


texts: Convert ST_Union polygon to PNG file in PostGIS?

texts:
I am working with Postgres 9.6, PostGIS 2.3 and GDAL 1.11.5. I would like to output the union of some polygons as a single PNG (or other image) file, with the polygons in black on a white background. (Not tiles - I just want an image.)

I know how to query the union of the polygons:

SELECT ST_Union(geom) FROM holdings WHERE town='Macclesfield';

But I don't know how to convert this to a PNG file. I've tried:

SELECT ST_AsPNG(ST_Union(geom)) FROM holdings WHERE town='Macclesfield';

But this errors:

ERROR:  function st_aspng(geometry) does not exist
LINE 1: select ST_AsPNG(st_union(geom)) from holdings where town='Ma...
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

I think maybe I need to convert the union into a raster layer first?


texts: How to read MapInfo text type using GDAL?

texts:
I'm trying to read MapInfo files using GDAL instead of MITAB library - https://github.com/mapgears/mitab.
With the MITAB C API I can call mitab_c_get_type() to identify the feature object type, and if it's a text type, get the text string by calling mitab_c_get_text().
Is there any way I can read this data by calling exported functions from GDAL? If this is a wrong approach, how would you handle reading MapInfo text type?


texts: geopandas to_postgis: ValueError: geom (geometry(MULTILINESTRING,4326)) not a string

texts:
I have a GeoPandas GeoDataFrame with a column named geom having MultiLineString values as such:
print(gdf.to_markdown()) # using package 'tabulate' to display a markdown table
|    | uid   | key   | install_year | update_date   | geom                                             |
|---:|------:|------:|:-------------|:--------------|:-------------------------------------------------|
|  0 | 57247 |  9247 | 1976-07-29   | 2018-01-21    | MULTILINESTRING (( Human readable coordinates )) |
|  1 | 57248 |  9248 |              | 2018-01-21    | MULTILINESTRING (( Human readable coordinates )) |
|  2 | 57249 |  9249 | 1991-12-24   | 2016-01-31    | MULTILINESTRING (( Human readable coordinates )) |
|  3 | 57250 |  9250 | 1927-02-15   | 2015-01-21    | MULTILINESTRING (( Human readable coordinates )) |
|  4 | 57251 |  9251 | 1991-12-23   | 2015-01-21    | MULTILINESTRING (( Human readable coordinates )) |
|  5 | 57262 |  9262 |              | 2015-01-21    | MULTILINESTRING (( Human readable coordinates )) |
|  6 | 57263 |  9263 |              | 2014-01-21    | MULTILINESTRING (( Human readable coordinates )) |

(tabulate pypi web page)
Now, I'd like to take benefit of the new .to_postgis() method in order to overwrite an existing PostGIS table (using an SQLAlchemy engine instance).
But when I do this:
import psycopg2
from sqlalchemy import create_engine
connector = psycopg2.connect(**database_parameters_dictionary)
engine = create_engine('postgresql+psycopg2://', creator=connector)
with engine.connect() as connection:
    gdf.to_postgis(
        "psql_table_name",
        connection,
        schema="public",
        if_exists="replace",
        index=True,
        dtype={'geom': Geometry(geometry_type='MultiLineString', srid=4326)}
    )

I encounter this error:
ValueError: geom (geometry(MULTILINESTRING,4326)) not a string

What does this mean, I don't get it (heck, geometries should not be geometries but strings?! This sounds strange for a postgis operation...), and how can I fix this error?
I'm on Python 3.6.9 / Ubuntu 18.04 using GeoPandas 0.8.1, sqlalchemy 1.3.20, and psycopg2 2.8.6 (dt dec pq3 ext lo64)


texts: How to make auto zoom when GeoJSON file is added on Leaflet map

texts:
I want it to auto zoom / zoom to center when any GeoJSON file is loaded on a Leaflet map. Help please.

$.getJSON("{{asset("storage/upload/shapefiles/districts/atd_district_level.geojson")}}",function(data){
                // add GeoJSON layer to the map once the file is loaded
                L.geoJson(data, {

                    style: {color: "#999", weight: 2, fillColor: "#00ad79", fillOpacity: .6},

                    onEachFeature: function (feature, layer) {
                        layer.bindPopup("<strong>" + feature.properties['DISTRICT'])
                    }

                }).addTo(map);
            });

I am new with Leaflet maps.


texts: Understanding difference between Coordinate System and Projection?

texts:
Can somebody please explain what is the difference between the Coordinate system (WGS 84 for example) and a Projection (Universal Transverse Mercator for example)?

What is the difference between a projected coordinate system and projected CRS


texts:
Whenever I try to open a kea file in TuiView, it throws me an error stating "NoneType object has no attribute GetOverviewCount".

What should I do to fix it?


texts:
I want to use two different GDAL versions on my Ubuntu 16.04. I installed GDAL via the UbuntuGIS ppa and via anaconda. My GDAL for python works fine, but when I want to use gdal_translate via the shell to convert a .jp2 to .tif the following error occurs:

  ERROR 4: `/media/lastro/.../L2A_T32ULB_20170430T103021_AOT_10m.jp2'
  not recognized as a supported file format

When I remove anaconda and have only gdal from the ubuntuGIS ppa installed then it works, but when I then install anaconda with the gdal binding the same error occurs.

Does someone know a workaround for this problem?


texts: Finding the distance between points along a river (Linear referencing)

texts:
I am new to gis and am struck. I have a network of rivers which is made up of small linestrings (polylines). In a different shapefile, i have locations of several stations on a river. I need to find distance along the river from each station till the end of river.

The problem is that I don't have one complete linestring of river. And there is no common attribute I can use to union or merge the linestrings. Also, some stations are on tributaries.

So the bottom line is I need to make a large network of rivers including tributaries and perform liner referencing. But I don't know how to do this. Any help is deeply appreciated. Please explain as clearly as possible.

Please comment if you need more information.


texts: Removing negative values from DEM using ArcGIS for Desktop?

texts: i'm using (ArcGIS 10.1) and i want to remove a negative values from DTM in order to keep just a ground topographic surface.


texts: Calculate distance in km to nearest points (given in lat/long) using ArcGIS DEsktop and/or R?

texts:
I have two point datasets in ArcGIS, both of which are given in WGS84 lat/lon co-ordinates and the points are spread across the whole world. I would like to find the nearest point in Dataset A to each point in Dataset B, and get the distance between them in kilometres.

This seems like a perfect use of the Near tool, but that gives me results in the co-ordinate system of the input points: that is, decimal degrees. I know I could re-project the data, but I gather (from this question) that it is difficult (if not impossible) to find a projection that will give accurate distances all over the world.

The answers to that question suggest using the Haversine formula to calculate distances using the latitude-longitude co-ordinates directly. Is there a way to do this and get a result in km using ArcGIS? If not, what is the best way to approach this?


texts: How to input an array in R script running inside QGIS

texts:
I am using QGIS 2.18.21 and running an R script which needs various files as inputs. The number of files to be inputted is also entered by user. I am thinking of taking the input of all file paths in an array. I am aware that in QGIS, we must specify the type of input at the beginning of script using ##. I have tried the following:

## file_arr = array
## file_array = list

but I got strange errors. The script didn't even run at all. For example, with a script as simple as the following:

##arr=array
print(arr)

I got the error:

and on clicking the View mesage log, I got the following:

Now, first of all I have no idea why is this a Python error since I am running a R script,
and secondly there is no proper documentation of how to specify the input of different R data types in QGIS.
Please help me out to take an array as input.


texts: Spatial Overlapping - Filter by % of Coverage

texts:
I am working with 2 Polygon Shapefiles: parcels & wetlands
Essentially what I want to accomplish is to filter the parcels dataset, where I only keep parcels that contain less than 50% Wetlands. For example the yellow parcels marked blue would be kept while the others would be filtered out, as the wetlands layer covers too much of the parcel

I've researched and tried various Spatial Joins, but none offer an option for %'s, and the gpd.overlay option only works on a DataFrame Level, not on individual polygons.
I was thinking some sort of comparison of Shape Areas, where:
filtered_parcels = filtered_parcels[(wetlands['shape_area'] <= parcels['shape_area'] * 0.5)]
However, I don't know how to obtain the Shape Area of the Wetlands contained only within a certain parcel.


texts: Does anyone have any references for this? (Is it OK to ask about it here?) Specifically the ImageControl API. I'm tasked with writing a program using it.


texts:
Working in R Studio I'm executing the lines:

setwd("D:/Perfil_IBGE/4. IBGE")
pasta = "1. Shapefiles Brasil"
shapefile_brasil = "BRUFE250GC_SIR"
shp <- rgdal::readOGR(dsn = pasta, layer = shapefile_brasil)

the result is: "no features found"

I work with Windows 10 and I checked about GDAL in my machine and I have this library installed.

What should I do?


texts: Android arcgis how to use only a layer from a service containing 3 layers?

texts:
I have this service:

http://*************/ArcGIS/rest/services/mod_streets/MapServer

Layers:

Numbers (0)
StreetsGraph (1)
Streets (2)

And I'm trying to use only Numbers (note /0)

    map = (MapView) findViewById(R.id.map);
    ArcGISDynamicMapServiceLayer layerStreets = new ArcGISDynamicMapServiceLayer("http://*************/ArcGIS/rest/services/mod_streets/MapServer/0");
    map.addLayer(layerStreets);

This is not working.

But this works (removing /0 but all layers are displayed into the map)

    map = (MapView) findViewById(R.id.map);
    ArcGISDynamicMapServiceLayer layerStreets = new ArcGISDynamicMapServiceLayer("http://*************/ArcGIS/rest/services/mod_streets/MapServer");
    map.addLayer(layerStreets);

I want to see only Numbers in my Android map...Any ideas?


texts: Calculate the exact area of polygons that fall within another polygon layer

texts:
I have two polygon layers. The first one contains the different “local districts” within a city, and the second one is the “greenery in different land use types”.

An important distinction to make, is that, it is not guaranteed that a particular greenery area, will always fit in a given district (i.e. a greenery area can be shared between 2+ adjacent districts).

When I join the “greenery in land use types” layer to the “districts” one, and I try ask to calculate the sum of my attributes, I get the greenery area sum for each district, I get double counting/summing on the "shared" areas.

Is there any way to consider the whole amount of greenery area that falls in each district, but to ignore the section of greenery that falls in another district?


texts:
I'm working on processing some geo-spatial raster data, JOG-A raster to be exact, and know the imagery is supposed to be to a scale of 1:250000. However I would like to calculate this in code but am not coming up with the correct values.

Here are the bounds of a single image and other values needed for calculations. (coordinates are in degrees, EPSG:4326)

Meters Per Degree: 111319
Image Size: x: 1536 y: 1536
DPI: 90
Upper Left:   Longitude: -89.5431 Latitude: 33.6529
Upper Right:  Longitude: -88.8578 Latitude: 33.6529
Bottom Right: Longitude: -88.8578 Latitude: 33.1351
Bottom Left:  Longitude: -89.5431 Latitude: 33.1351

I thought I could say

(degLatA - degLatB) * meterPerDeg / imageSizeY * dpi / 0.0254
33.6529deg - 33.1351deg = 0.5177deg
$1 deg * 111319m/deg = 57634m
$2 m / 1536px = 37.52m/px
$3 m/px * (90px/in * 1in/0.0254m) = 132953

which gives a scale of 1:132953 ... not even close!
The units all cancel great, but as you can see from the calculation, the value is not even close to 250k.

Can anyone explain where I'm wrong?

EDIT:
The math for calculation is all correct, I just needed the correct DPI which for CADRG imagery is 169 according to the SPEC documents. (buried on page 10 if you want to verify or read for yourself).


texts: Improving Sentinel 1 Interferometric DSM Accuracy - Removing Surface Anomalies

texts:
I have been working through the ASF's "Recipes" for the creation of a DSM from Sentinel 1 (A and B) imagery.

I am able to produce a DSM, and am relatively happy with the results in terms of the surface produced - it appears to follow the trends of my known datasets (hills where there should be hills and valleys where there are valleys).

I have however noticed artifacts or static in some areas of the surface produced. Please see the image below showing what I mean - basically areas of spikes and troughs in random places on the surface. These areas tend to be consitent amongst different periods too (for example I produced a 2018 surface and the same issues appeared in my 2016 surface).

The image above shows an example of a cross section through an area. As can be seen, the elevation drops from around 300m to less than 100 and then back up to 400m. All within around 400m in the horizontal.

At moment I have been using the default values in the ESA SNAP Toolsets.
My question I guess can be broken into two parts. 1) What created this "static" in the image - is it topography (large hills creating interference) or a setting / value in SNAP I should be tweaking? 2) If it is a setting, what should I be looking at? The process of generating the DSM is a long one, and having to repeat each step is time consuming, so a general idea of where I should focus my attention would help.


texts:
I'm trying to create a map image with a coordinate grid. I have not been able to get a grid to appear that displays according to my CSR.
I have used this a couple of guides: https://geoafrikana.com/how-add-grids-and-grid-frame-map-qgis/ , https://www.qgistutorials.com/en/docs/making_a_map.html.
I have not had any luck when going to add a grid. When the interval is selected for "map units" there is no display, even after playing with different interval numbers for X and Y (i.e. 2, 10, 30, 100). I get a display when using "centimetres" or "millimetres" but that does not appear accurate to the CSR.
My CSR is GDA94.
What is wrong here? Here's an image for extra reference.


texts: De-highlight (change opacity of) other features that don't share same values

texts:
I have this Arcgis Online webmap that is a base for my testing in Arcgis Javascript. I have not done anything with the javascript file, other than the basic setup, as I am new to this. I was wondering if it was possible to de-highlight all features that dont share a color with it when hovering over them?

The data in the picture below consists of the same layer copied twice in a WebMap. The first layer filters out some points that becomes squares, and the other layer becomes circles. Both layers origin from the same feature layer, and the colors are determined by a string attribute of 3 characters. So if I hover over any of the blue points, my ideal map should look like the second picture, where all other points have an opacity of 80%. If no points are hovered, all should be visible (0% opacity).

main.js:
require([
    "esri/WebMap",
    "esri/views/MapView",
    "esri/widgets/Popup"
], function (
    WebMap,
    MapView,
    Popup
) {

    var map = new WebMap({
        portalItem: {
            id: "xxxxx"
        }
    });

    var view = new MapView({
        container: "viewDiv",
        map: map,
        popup: {
            dockOptions: {
                     buttonEnabled: false
                 }
             },
        ui: {
            components: ["zoom", "compass", "attribution"]
        }
    });

});


texts: Raster distance for statistical analysis

texts:
I am trying to some kind of hedonic analysis. I need to calculate driving forces for new urban areas. In order to do this, i have a raster with new urban area's and one with existing urban areas (30x30 meter), 6 different road types, and the population growth in the smallest administrative area. My hypothesis is that existing urban areas and population growth are the most important drivers for new urban areas.

However, to do the analysis in Stata i need the distances from each cell to each other feature (existing urban and roads) and in each cell the population growth corresponding to that administrative area. Then I'll do a regression analysis in Stata in which the dependent variable is urban or not (1 or 0). The distances and population growth per cell are then the independent variables.

My question is, how can i calculate these distances? and have a resulting table with this? I might need to change the resolution, since i have 254 million cells (entire Bangladesh).

I just did the distances using the Generate Near Table-tool, however that gave me distances from the clustered-new-urban-cells-in-polygons to the line-features (roads). And I don't have a dependent variable. So there I made a thinking error.


texts: How to Fix Inconsistent Data Conversions in SFCGAL Expressions


texts: Batch raster operations to get Maximum Precipitation in Catchment

texts:
I have information of maximum precipitation per grid cell for a certain area over 55 years, each year is a different matrix with the corresponding precipitation values. I know that I need to do the following steps:

Convert the precipitation points into a raster
clip the precipitation raster with the catchment area
Get the maximum value of precipitation for 1 year from the clipped raster.

However, I don't want to do this for each year, given that I have a total of 5 catchments so I will have to do this procedure 5*55=275 times.

Is there a way to program this operation in qgis/Grass and get a vector with the maximum precipitation per year for a given catchment?


texts: ArcGIS 10.1 Converting a model tool to a script tool

texts:
I have recently created a model tool for a class project and my instructor would like it to be converted to a script tool for practice purposes. I have exported the script from my model and loaded it into a script within a toolbox but there are several things I am unsure about. I have defined the parameters as they were in the model but when I run the script tool it give me multiple errors. I know that the model doesn't pass model only tools to the script (i.e. Select Data) so I'm sure my problem stems from that but I'm not sure what exactly I need to add to fix the problem.

I will attach a copy of my script, tool parameters, tool properties, the original model and the errors I get. I am an amateur with Python so bare with me if I'm asking too much.

NOTE: the '*' represents the actual values but can't be shared

# ---------------------------------------------------------------------------
# CROB1.py
# Created on: 2013-10-08 11:45:20.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: CROB1 <ArcSDE_Connection_File_Location> <Select_Your_AP_County> <Select_Your_AP> <Your_Counties> <Your_AP>
# Description:
# This tool creates the layers for management plan map for your AP boundary. The user must still format the Layout View before exporting.
# ---------------------------------------------------------------------------

# Set the necessary product code
import arceditor
# Import arcpy module
import arcpy

arcpy.env.overwriteOutput = True

# Script arguments
ArcSDE_Connection_File_Location = arcpy.GetParameterAsText(0)
Select_Your_AP_County = arcpy.GetParameterAsText(1)
Select_Your_AP = arcpy.GetParameterAsText(2)

# Local variables:
FDEP_SDE_Connection_sde = ArcSDE_Connection_File_Location
Counties = FDEP_SDE_Connection_sde
Aquatic_Preserves = FDEP_SDE_Connection_sde

# Process: Create ArcSDE Connection File

arcpy.CreateArcSDEConnectionFile_management(ArcSDE_Connection_File_Location, "FDEP_SDE_Connection", "******", "****", "", "DATABASE_AUTH", "****", "*****", "SAVE_USERNAME", "SDE.DEFAULT", "SAVE_VERSION")

# Process: Select Data
arcpy.SelectData_management(FDEP_SDE_Connection_sde, "BASE.FL_SHORELINE_AREAS")

# Process: Select
arcpy.Select_analysis(Counties, Your_Counties, "\"NAME\" = '%Select Your AP County%' ")

# Process: Select Data (2)
arcpy.SelectData_management(FDEP_SDE_Connection_sde, "DEP.AQUATIC_PRESERVES")

# Process: Select (2)
arcpy.Select_analysis(Aquatic_Preserves, Your_AP, "\"LONG_NAME\" = '%Select Your AP%'")`


texts: ArcGIS Geoprocessing Script runs fine in Desktop but crashes as Geoprocessing Service?

texts:
I have a small and simple geoprocessing script that I am trying to expose as a geoprocessing service in ArcGIS Server 10 SP4.  The script is below.  Pardon the hard-coding...

With the hard-coded variables, this script runs perfectly fine in desktop.  However, once published as gp service, it crashes in a bad way, like I can't even trap the error.  The exact message in the GIS Server log - "Container process 6436 has crashed on machine westchamp24-pc."

A couple of potential causes I can eliminate right off the bat:

Permissions:  The entire D: drive is wide open to the 'Everyone' user
Licensing:  ArcGIS Server is licensed at Enterprise level

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute CostDistance
outCostDistance = CostDistance("origin",SOURCE,"#","D:/RasterStuff/backlinkRaster.tif")

# Execute CostPath
outCostPath = CostPath("destination", outCostDistance,"D:/RasterStuff/backlinkRaster.tif")

# Convert Result to Polyline
arcpy.RasterToPolyline_conversion(outCostPath, "D:/RasterStuff/costPath.shp")
featSet = arcpy.FeatureSet()
featSet.load("in_memory/costPath")

arcpy.SetParameter(0,featSet)

except:
    # do stuff with error
    print
finally:
    print


texts: Create GeoJSON polyline with elevation from a shapefile

texts:
Here is my situation: I have volunteers using smartphones to collect bike trails. These GPX files are a good start, but they need editing due to canopy cover, etc. I export these to shapefiles and edit to match the actual trails where I can see them with the latest ortho imagery, get rid of spikes, other issues. However, the edited shapefile have lost the elevation (and time) data from the original GPX file. I want to export this shapefile as a GeoJSON to show it in a web map. The resulting geometry should have z values like this:

"geometry": {
  "type": "LineString",
    "coordinates": [
      [
        -83.43572,
        39.1932,
        272.9
      ],

Is there a way to take the polyline shapefile, add elevation information, and export that to GeoJSON. Maybe first to a point, then to a GeoJSON linestring? If there is a geojson editor that may work as well.


texts: Create a shapefile from dataframe in R, keeping attribute table

texts:
I am trying to create a shapefile from a .csv file in R. I can successfully create a working shapefile with this code, but when I try to view the attribute table for it in arcMap 10.4, the table only contains columns for objectID, latitude, and longitude. I need the info in the other columns for symbolizing the data correctly.

Is there a way to specify which columns of the data to keep in the shapefile?

(I know how to make a data frame into a shapefile within arcMap, but I would like to automate this step in R if possible.)

library(maptools)
library(rgdal)
library(sp)

#data
site <- c("a","b","c","d")
prop_c <- c(0.88,0.48,0.15,0.47)
prop_b <- c(0.17,0.18,0.09,0.08)
minus_c <- 1-prop_c
minus_b <- 1-prop_b
lat <- c(44.22,38.38,33.35,43.48)
long <- c(-124.45, -123.70, -124.40, -124.05)
MyData <- cbind.data.frame(site, prop_c, prop_b, minus_c, minus_b, lat, long)

#convert data to shapefile

WGScoor<-  MyData #data to convert
coordinates(WGScoor)=~long+lat #column names of the lat long cols
proj4string(WGScoor)<- CRS("++proj=longlat +datum=WGS84") # set coordinate system to WGS
WGScoor.df <- SpatialPointsDataFrame(WGScoor, data.frame(id=1:length(WGScoor)))
LLcoor<-spTransform(WGScoor.df,CRS("+proj=longlat"))
LLcoor.df=SpatialPointsDataFrame(LLcoor, data.frame(id=1:length(LLcoor)))
writeOGR(LLcoor.df, dsn=getwd(),layer="MyShapefile",driver="ESRI Shapefile")
#successfully creates shapefile that opens in arcGIS,
#but shapefile does not include other columns in original data frame


texts: Using intersect/clip for polyline to polyline features resulting in Warning 000117 ArcMap?

texts:
I am trying to measure the length of a polyline from its origin to its intersection with another polyline. The application is to find out the distance from a meteorological station to the coastline. I created polylines with origin in the met station outwards every 30 degrees. When I use clip/intersect I get 'Warning 000117: Warning empty output generated'. I have tried the online help for related topics and still have not solved it (exporting to a shp, to a new gdb, re-starting ArcMap, creating a new project, resetting extents, checking al layers have same coordinates).

Files I am trying to use: ArcMapFolder

https://drive.google.com/folderview?id=14paLxlLnTMGhzTB35pRjrWd0WQkKNLiK

I have tried also the solutions suggested in Intersect (analysis) of polygons generating empty output


texts: What annoys you about your GIS Environment?

texts:
I have a question regarding your favourite GIS environment.  Over the past few decades, there have been the usual major players, MapInfo, ESRI, Intergraph (Grass?), and we all love and hate aspects of the major players of GIS.

The purpose of my question is to provide some forum (community wiki?) where we can place suggestions that better our working lives.  So please be specific.

Please try and be constructive, give descriptions of process, data, expected outcomes and scenarios.

  Constructively, What are your primary issues with your GIS package of
  choice?

The Seed for the Question
Looking at the Hadoop project by Apache and some of the impressive work completed there, which started the thinking cogs.  Since GIS is such a computationally intensive profession, how has this technology not being incorporated in our favourite GIS software!  There are some fantastic solutions that could enhance a GIS professionals life.

The other purpose is to look at the possibility of taking some of these ideas and potentially incorporate them in FOSS4G projects (QGIS? GeoServer? New Project???), increase usability and functionality and give the big guys a run for their money ;)


texts:
I've got 2 layers being rendered separately with this:

var layerSource = {
        user_name: 'jonmrich',
        type: 'cartodb',
        sublayers: [{
            sql: "SELECT * FROM " + carto_table + " WHERE type = 'red_blank'",
            cartocss: '#location_fusion_table {marker-file: url(https://s3.amazonaws.com/com.cartodb.users-assets.production/production/jonmrich/assets/20150203194453red_pin.png);marker-width:25;marker-allow-overlap: true;}'
        }, {
            sql: "SELECT * FROM " + carto_table + " WHERE type = 'small_blue'",
            cartocss: '#location_fusion_table {marker-fill:#1d5492;marker-width:5;marker-line-width: 0.5;marker-line-color:#fff;marker-allow-overlap: true;}'
        }]
    };

and later:

cartodb.createLayer(map_object, layerSource, options)
        .addTo(map_object)

Here's what ends up happening:

What I'd like to have is where the red pins ALWAYS show on top of the blue circles. If I turn overlap to false, I lose a lot of detail, but the blue circles aren't on top of my red pins. Is there some sort of "z-index" or "bring to front" that I can employ here?


texts: QGIS 3.x: Show only labels / TilesXYZ TileLayer of labels only

texts:
I'm using QGIS 3.6 and the TilesXYZ extension (deprecated for 3.8+) to export tilelayers of my map. I'd like to export a base layer without labels and a label layer to show on top of features.

The first part is possible, using a one-click label toggle button:

action = QAction(QIcon(""), "Turn labels" + "\n" + "ON/OFF", iface.mainWindow())
action.setCheckable(True)
iface.addToolBarIcon(action)

def label_control():
    for layer in QgsProject.instance().mapLayers().values():
        if layer.type() == QgsMapLayer.VectorLayer:
            if action.isChecked() == True:
                layer.setLabelsEnabled(True)
            else:
                layer.setLabelsEnabled(False)
        layer.triggerRepaint()

action.triggered.connect(label_control)
# Uncomment line below if you want to remove the icon yourself,
# otherwise it will be removed automatically when you restart QGIS
# iface.removeToolBarIcon(action)

Is there a similar way to hide all my features but show labels?


texts: OpenLayers-plugin error after installing update

texts:
After installing an update for the OpenLayers-plugin today, this plugin does not work any longer. Reinstalling did not function neither with QGIS 2.12.0-Lyon nor with QGIS 2.14.1 ESSEN. Following error message occurred:

"Konnte Erweiterung openlayers_plugin nicht laden aufgrund eines Fehlers beim Aufruf der initGui() Methode

UnicodeEncodeError: 'ascii' codec can't encode character u'\xdc' in position 11: ordinal not in range(128)
Traceback (most recent call last):
  File "C:/PROGRA~1/QGISLY~1/apps/qgis/./python\qgis\utils.py", line 315, in startPlugin
    plugins[packageName].initGui()
  File "C:/Users/Stalze/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py", line 115, in initGui
    print action.text()
UnicodeEncodeError: 'ascii' codec can't encode character u'\xdc' in position 11: ordinal not in range(128)

Python-Version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS-Version: 2.12.0-Lyon Lyon, cd9d645

Python-Pfad:
C:/Users/Stalze/.qgis2/python/plugins\processing
C:/PROGRA~1/QGISLY~1/apps/qgis/./python
C:/Users/Stalze/.qgis2/python
C:/Users/Stalze/.qgis2/python/plugins
C:/PROGRA~1/QGISLY~1/apps/qgis/./python/plugins
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg
c:\osgeo4~1\apps\python27\lib\site-packages\python_dateutil-2.1-py2.7.egg
c:\osgeo4~1\apps\python27\lib\site-packages\six-1.3.0-py2.7.egg
C:\PROGRA~1\QGISLY~1\bin\python27.zip
C:\PROGRA~1\QGISLY~1\apps\Python27\DLLs
C:\PROGRA~1\QGISLY~1\apps\Python27\lib
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\plat-win
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGISLY~1\bin
C:\PROGRA~1\QGISLY~1\apps\Python27
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\PIL
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/Stalze/.qgis2//python
C:\Users\Stalze\.qgis2\python\plugins\DigitizingTools\tools
C:\Users\Stalze\.qgis2\python\plugins\mmqgis/forms
C:\Users\Stalze\.qgis2\python\plugins\SelectTools/tools
C:/Users/Stalze/.qgis2/python/plugins\shpsync\xlrd-0.9.4-py2.7.egg
C:/Users/Stalze/.qgis2/python/plugins\shpsync\xlwt-1.0.0-py2.7.egg
.
C:\PROGRA~1\QGISLY~1\apps\qgis\python\plugins\fTools\tools"

As I work with the Bing aerial layer quite often, what could I do to make this plugin run again?


texts: Draw horizontal lines on points shapefile using ArcGIS

texts:
I have one shapefile of 200 points in which I would like to draw one horizontal line of 2Km on each of those points.

How can I do this?


texts: In the attempt to install Crayfish plugin from QGIS Plugin Repository, I often encounter such problem as "Crayfish depends on a platform specifi complied library which was not found. Would you like to automatically download and install one from the developer's website?" By clicking on Yes, a window with a message "Download and installation successful" will be seen. But by clicking "Ok", another message saying "couldn't load plugin crayfish due to error when calling its initGui() method" would be displayed.


texts: Getting HiRes Tiles via WMTS results in Incorrect size of download block?

texts:
To use maps in prints we'd like to use HiRes Tiles offered by a server (basemap.at) providing WMTS tiles sized 512x512.

The example on https://www.gdal.org/frmt_wmts.html works well for standard tiles sized 256x256 but fails for the HiRes (512x512) resolution resulting in the following error:

ERROR 1: GDALWMS: Incorrect size 512 x 512 of downloaded block, expected 256 x 256, max 256 x 256.
ERROR 1: GDALWMS: ReadBlockFromFile (https://maps1.wien.gv.at/basemap/bmaphidpi/normal/google3857/19/182218/285835.jpeg) failed.
ERROR 1: GDAL_WMS>, band 1: IReadBlock failed at X offset 0, Y offset 0
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 0

Are we missing any (undocumented) tags in the WMTS XML?

wmts.xml:

<GDAL_WMTS>
<GetCapabilitiesUrl>https://www.basemap.at/wmts/1.0.0/WMTSCapabilities.xml</GetCapabilitiesUrl>
  <Layer>bmaphidpi</Layer>
  <Style>normal</Style>
  <TileMatrixSet>google3857</TileMatrixSet>
  <ZoomLevel>19</ZoomLevel>
  <Format>image/jpeg</Format>
  <DataWindow>
    <UpperLeftX>1810893.1</UpperLeftX>
    <UpperLeftY>6109237.0</UpperLeftY>
    <LowerRightX>1811214.4</LowerRightX>
    <LowerRightY>6109038.2</LowerRightY>
  </DataWindow>
  <BandsCount>4</BandsCount>
  <Cache />
  <UnsafeSSL>true</UnsafeSSL>
  <ZeroBlockHttpCodes>204,404</ZeroBlockHttpCodes>  <ZeroBlockOnServerException>true</ZeroBlockOnServerException>
</GDAL_WMTS>

Python script:

from osgeo import gdal
ds = gdal.Open('wmts.xml')
gdal.GetDriverByName('JPEG').CreateCopy('map_z19.jpg', ds)

GDAL version 2.2.2, Python 3.6


texts: Erreur gdal_polygonize.py package ARfun with R

texts:
I'm a french student, in advance I'm sorry for my english. I'm installed the package APfun on R to transform Raster to Polygons with Appolygonize function. But I have a error :

crownsPoly <- raster(paste(getwd(), "/Output_data/crownpoly.tif", sep = ""))

   APpolygonize(crownsPoly,
+              readToMemory = TRUE,
+              outFile = paste(getwd(), "/Output_data/segmentation.shp", sep = ""),
+              OSGeoPath = "C:\\OSGeo4W64",
+              connectivity = 4)

Error in APpolygonize(corwnsPoly, readToMemory = TRUE, outFile = paste(getwd(),  :
  Could not find required file 'C:\OSGeo4W64/bin\gdal_polygonize.py'

This function need to OSGeo4W (for GDAl), I have realised this installation with https://trac.osgeo.org/osgeo4w/ like hep package advised.

I find gdal_polygonize.py within the folder C:\OSGeo4W64\apps\Python37\Scripts. But the function research this file within C:\OSGeo4W64/bin\gdal_polygonize.pyso it is nowhere to be found.

Why OSGeo4W installation put this file within bad localisation ? When I put this file in the good folder (where the function search) it's a other file which is absent...

I realised a new installation of OSGeo4W6 and this error was yet present...

When I search gdal_polygonize.py in r, I have nothing :

> Sys.which("gdal_polygonize.py")
gdal_polygonize.py
                ""

Do it have other manipulation between OSGeo4W64 and the launch of the function APpolygonize ?

I use Windows, Rstudio with "R version 3.6.1 (2019-07-05)" and I have Python 3.7.3 in my computer.


texts: Convert a mask image to polygon shapefile using GRASS

texts:
I used image classification technique in GRASS GIS v7 and generated the image mask below in which white area has NULL value and red area contains 1.

I am going to clip those red color areas in a LiDAR file so I need to convert those to a polygon shapefile. I know that I need to use a combination of image analysis techniques to achieve this but don't know whether GRASS has similar functionality to generate polygon shapefile from a mask.


texts:
I need to split large contour polyine datasets every 100m to be used for tiling. It currently has each (25cm) contour as a single (very long) line that extends beyond the tile and impacts labeling etc.

Is it possible to break it up using MapBasic to 100m segments?

I can see lots of advice for a manual process (https://groups.google.com/forum/#!topic/mapinfo-l/-T0SX4lqMXk) but this is not practical in such large datasets.

Note
The data is multiple mapinfo extended files of 8 to 12 GB each. So the process has to work in MI 16, 64 bit.

The Use-Case
The 25cm contour lines in SSA are very long and this impacts on performance of the SSA server. It also results in poor labelling of corner objects in the map window (typically only labelled in top LHS corner tile).

UPDATE
We ran the process on a 1.5GB Tab file in 32 bit and imported it into SSA. It fixes the issues we faced but the issue is with the mapinfo extended files and the dll error when the toolbox is opened 64-bit.

" MapBasic applications that call a 32 bit DLL or other external code will need a 64 bit version of the DLL. This sort of application will require modification to run in the 64 bit version of MapInfo Pro. "

We've tried to bypass the DLL altogether, but there's function declarations that are needed

We've also tried to access the dll's to extract what I need - unsuccessfully.

What can we do?


texts:
I have numerous HDF5 files(.nc) which I need to batch process warp using one of the gdal utilities gdalwarp. When I tried to warp the files an error occurred:

INPUT:

gdalwarp -geoloc -te 109.975 3.475 135.025 25.025 HDF5:"@file"://geophysical_data/chlor_a %out_path%\@fname.tif"

RESULT:

ERROR 1: Unable to compute a GEOLOC_ARRAY based transformation between pixel/lin
e
and georeferenced coordinates for HDF5:A2015045060500.L2_LAC_OC.nc://geophysical
_data/chlor_a.

Update1:

Just to make it clear, do you mean in lat.vrt, lon.vrt and chlor.vrt I should remove the GCP Id's and MDI key and insert this section:

   <metadata domain="GEOLOCATION">
     <mdi key="X_DATASET">oc-long.vrt</mdi>
     <mdi key="X_BAND">1</mdi>
     <mdi key="Y_DATASET">oc-lat.vrt</mdi>
     <mdi key="Y_BAND">1</mdi>
     <mdi key="PIXEL_OFFSET">0</mdi>
     <mdi key="LINE_OFFSET">0</mdi>
     <mdi key="PIXEL_STEP">1</mdi>
     <mdi key="LINE_STEP">1</mdi>
   </metadata>

in between this section?

<VRTDataset rasterXSize="1116" rasterYSize="1610">

  ###### metadata section here #######

 <VRTRasterBand dataType="Float32" band="1">
    <SimpleSource>
      <SourceFilename relativeToVRT="0">HDF5:A2015194044000.L2_LAC.SeAHABS.nc://geophysical_data/chlor_a</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1116" RasterYSize="1610" DataType="Float32" BlockXSize="373" BlockYSize="17" />
      <SrcRect xOff="0" yOff="0" xSize="1116" ySize="1610" />
      <DstRect xOff="0" yOff="0" xSize="1116" ySize="1610" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>


texts: Placing elevation numbers on contours with uphill orientation and position in QGIS?

texts:
I want to orient my elevation labels uphill.

Is that possible?

I'm currently using the labeling tool from QGIS (not that one in layer settings). I want the labels curved and uphill oriented (and I would like to define how regularly the numbers are displayed). It's all for a print-map.


texts: Creating two-colored map of region in one color, and neighbors in another color?

texts:
As they say, a picture is worth more than a thousand words:

Turkey in Wikipedia

What I need is to create maps for an educational system, and though we can create them manually in a lot of environments, my requirement is to make maps dynamically based on different criteria.

A simple starting point is to create a two-colored map to show the region of interest in an attractive color, among its neighbors which are rendered in gray color. This image from Wikiepdia is another example.

I have the WKT database of countries, and cities, and other geographical data. How can I create images from them this way?

I've tried SQL Server's spatial results pane, yet it renders regions in a multicolored map, with random color given to any region. Not useful at all.
I also searched the web to find some online WKT rendering pages, but they also don't give me the flexibility of programming, and also their background map is Google Map with no access to make it a flat simple map.


texts: Selecting all pixels from Pixel Inspector?

texts: While using the Pixel Inspector, whenever I click "copy all values" it only copies values that you can see on the screen. Anyone know a solution? I need all the values to count the number of times 1, 2, and 0 appear on the graph.


texts: How do you enable HTTP GET for submitJob for a GPService?

texts:
I published a GPService to ArcGIS Server 10.2--the service is a copy of an existing service with a different name. I receive no errors during publishing, however I get 404 errors when attempting to submit a job via HTTP GET.  HTTP Post works fine--job is successfully submitted and it completes successfully.

I did not set up the original service, so I'm thinking I missed a checkbox or some configuration during the publishing process.  Is there a way to enable HTTP GET for a GPService after it has been published?  If not, how do you enable during publishing?


texts: I nead to find the nearest places based on latitude/longitude points

texts:
I have a java aplication backed by a mysql database which contains events with locations (latitude/longitude), the location data is saved as POINT data type.

I nead to find the nearest events in some radius. I am already able to do this but the results are innacurate.

I'm also using hibernate spatial to make the query.

I'm using the folowing libraries :
1.) com.vividsolutions.jts  ... for creating geometry types
2.) Geotools org.geotools ... for calculating distance between points

First I make a circle Geometry with a center placed in a current location this circle is given some radius. For example 500m if I want to find the events in the distance of 500m. The problem is that the Geometry construct takes neads the radius in degrees. If i'm not mistaken the distance in degrees changes based on a point on earth.

then I query the database with the within hibernate spatial function, but the results are unnacurate.

So I dont know how to calculate the distance in degrees based on location on earth. I think there should already be some library that is doing this, but for now I havent find an example. Can you advise me some code?

My Code:

//calculate the shape where the events should be
Geometry circle = createCircle(currentLongitude,currentLatitude, radiusInKm);
//retrieve the events from database that are within the calculated shape
List<si.arctur.model.Event> seznam = find(circle.toString());

    //the methods
private static Geometry createCircle(double x, double y, final double RADIUS) {

          double latitude = 110.54;
          double longitude = 111.320;

          //TODO beter accuracy
          double radiusInDegreesLatitude = RADIUS / latitude;
          double radiusInDegreesLongitude = RADIUS / Math.abs((longitude * Math.cos(latitude)));

          GeometricShapeFactory shapeFactory = new GeometricShapeFactory();
          shapeFactory.setNumPoints(1000);
          shapeFactory.setHeight(radiusInDegreesLatitude * 2);
          shapeFactory.setWidth(radiusInDegreesLongitude * 2);
          shapeFactory.setCentre(new Coordinate(x, y));//there are your coordinates
          Polygon shape = shapeFactory.createCircle();
          return shape;
    }

    private List find(String wktFilter) {
        Geometry filter = wktToGeometry(wktFilter);
        Query query = em.createQuery("select e from Event e where within(e.location, :filter) = true", si.arctur.model.Event.class);
        query.setParameter("filter", filter);
        return query.getResultList();
    }

    private Geometry wktToGeometry(String wktPoint) {
        WKTReader fromText = new WKTReader();
        Geometry geom = null;
        try {
            geom = fromText.read(wktPoint);
        } catch (ParseException e) {
            throw new RuntimeException("Not a WKT string:" + wktPoint);
        }
        return geom;
    }


texts: Classifying coordinates on land/on water using QGIS?

texts:
I am a new user of QGIS. I'm looking for a solution for my problem. I have my own dataset providing the latitude and longitude of different unique ID's. Now, I would like to classify those ID's that are showing on the rivers as 'on water' and those that are on land as 'on land'. There are thousands of data points, so manual classification is not possible.

Is there any way to load in some data of rivers, borders and/or coastlines in order to make classification like 'on land' or 'on river' possible?

I found a dataset providing the routes of waterways. Although this data does not show the true width of rivers, it might be useful as an example: http://nationaalgeoregister.nl/geonetwork/srv/dut/catalog.search#/metadata/2f58c248-b229-48af-a838-3279f7f0173c