texts: Is there an xsd schema or UML diagram for GDF and if not then why?

texts:
I am a GIS newbie exploring the GDF Geographic Data Files format.
The question "GDF — Can anyone suggest some good reading?" was a great starting point (despite dead links that have lots of dead copies on the archive), as was this from openstreetmap and this from Safe, makers of FME.
I feel it would be really useful to browse a visual UML diagram of the format (perhaps by opening the xsd schema), more or less alongside all the other reading. The best I could find so far:

This question took me to lots of great examples for other
formats but not GDF.
This page on github has the description "Information model and GML schema for GDF" but does not seem to be what I'm after. From my understanding of xsd schemas in general, I was under the impression I would be able to find a single .xsd file for the GDF format, like I was for KML, and that page does not seem to have that, unless I am failing miserably to spot it?

So is there an xsd schema for GDF and if not then why? The only reason I can think of would be that it's because it's controlled by an ISO standard - but intuitively, this doesn't stop a UML diagram, perhaps for an old version, from being publicly available?


texts: pgRouting gives fragmented results when trying to calculate route

texts:
I have a road network and I used pgr_createTopology and pgr_nodeNetwork to convert it for use in pgRouting library. However the results are fragmented and there is no continuous path. Furthermore I have to give source and target values instead of ID in the below query. here 747 and 311 are source and target values from table rather than IDs.

pgr_dijkstra('
    SELECT id, source::int, target::int, length::double precision AS cost FROM guru_noded',
    747, 311, false, false);

I am new to this.


texts: Qspatialite plugin problems with queries using inner join

texts:
I'm trying to create a layer in QGIS from a sqlite database using a query from 2 tables: one has spatial data and the other doesn't.

When I build the query using Qspatialite using the 'Advance SQL Editor' it throws an error saying that the first column in the query has an ambiguous column name. Copying and pasting from an MSAccess query, it finds the matching records without a mutter.

Here are segments of the two queries:

MSAccess version (working)

SELECT Inspections.ID, Inspections.InspectionDate, Inspections.InspectionBy, Inspections.PathID, Inspections.Section, Inspections.GISLink, Inspections.Observations, Inspections.Notes, Inspections.ActionPriority, Sections.Geometry AS Geometry
FROM Inspections INNER JOIN Sections ON Inspections.GISLink = Sections.GISLink;

Qspatialite version (not working)

SELECT 'Inspections'.'InspectionDate',
'Inspections'.'InspectionBy',
'Inspections'.'PathID',
'Inspections'.'Section',
'Inspections'.'GISLink',
'Inspections'.'Observations',
'Inspections'.'Notes',
'Inspections'.'ActionPriority',
'Sections'.'GISLink',
'Sections'.'Geometry' AS Geometry
FROM 'Inspections', 'Sections'
inner join 'Inspections'
on 'Inspections'.'GISLink' = 'Sections'.'GISLink'

Could anyone give me an idea of what wrong's here - when I remove the INNER JOIN the query runs, but obviously gives me the duplicate values you'd expect!

Is this a problem with QSpatialite in the way it constructs the queries with '' all over the place?

I've just tested both queries in QGIS DB manager and same result the copy-paste version works, but the 'Advance SQL editor' version throws the same error...

Many thanks - and hoping that I can save someone else frustrating hours of SQL - I was assuming the differences in MSAccess would prevent the query from running in Qspatialite and as a last ditch I did the copy-paste. Doh!

I'm using QGIS 1.8.0 and QSpatialite 6.0.4


texts: Error publishing feature service in GISserver 10.2

texts:
I installed ArcGIS 10.2 for Server just now. I started to create a map using the ArcMap application of ArcGIS for Desktop.

I created My map with just: 1 shape file + 1 gdb file...

there were not any error when I push the anlyze button(I mean I have checked the enterprise db for my gdb file to be published correctly).

After all, when I want to publish my map as a feature service(I mean a service map with the 'feature access' capability enabled) I encounter an error like this(the feature service definition is created with no errors but just not published and running):

I searched all the web but there where no solution for the error!!!
actually there where no useful log in GP results window as you see below:


texts:
I'm currently using QGIS 1.8.0

I would like to "merge" a raster data set on precipitation with a shapefile of administrative districts.

Basically my the aim is to figure out the precipitation in different administrative areas. The data I have on precipitation is formatted in a lat/lon 0.5 grid. Each administrative district however intersects several of these 0.5x0.5 squares. How can I connect each administrative district to multiple corresponding 0.5x0.5 squares and then weigh the value of the squares that are only partially contained in the administrative district/ split up the squares that are divided between several districts?

Thanks!


texts:
I have saved NaturalEarth data set in my desktop. Then I created a view as below.

        <script>
            var viewer = new Cesium.Viewer('cesiumContainer', {
                imageryProvider: Cesium.createTileMapServiceImageryProvider({
                    url: 'C:\Users\name\Desktop\NaturalEarthII'
                }),
                baseLayerPicker: false,
                geocoder: false
            });
        </script>

But the viewer does not show the map. What I want is to use Natural Earth data set offline. But if I saved the data set in 'Assets/Textures/NaturalEarthII' in cesium folder in the project it works fine with the below code.

 var viewer = new Cesium.Viewer('cesiumContainer', {
                imageryProvider: Cesium.createTileMapServiceImageryProvider({
                    url: Cesium.buildModuleUrl('Assets/Textures/NaturalEarthII')
                }),
                baseLayerPicker: false,
                geocoder: false
            });

What's the wrong I have made?


texts: Computing shortest distance to line segment from set of points using QGIS?

texts:
I have a shapefile that describes the urban boundaries of a city and I want to calculate for a set of points the nearest distance to the boundary. Given the number of points is very large (millions) I decided to first test my planned procedure to make sure it is producing the correct output. This is what I did:

Created a grid of test points using 'Create grid'
Converted the boundary to a set of points using 'Convert lines to points'
Identified the shortest path to the boundary using 'Distance to nearest hub (lines to hub)' with meters as the measurement unit

I have encountered a problem in that for a number of test points, particularly those for which the shortest path to the boundary is in a north-south direction, the nearest point on the boundary identified by 'Distance to nearest hub' is incorrect (ie. there are other points on the boundary that are closer).

As an example, for the selected test point, the closest point on the boundary identified by 'Distance to nearest hub' (ie. the northern most point) is 15037.985 meters away, but in fact there is a closer point (the western most point) that is only 14679.774 meters away. I have also checked these distances (and a number of others) against Google Maps and they seem to be fairly accurate (within 10-20 meters).

Not sure what is going on here. I am relatively new to QGIS.  I am using QGIS version 3.2.3-Bonn.


texts: How do I create a new vertex in the middle of an edge in postGIS?

texts:
I have a postGIS database populated with OSM data with the osm2po import tool. I have a number of fixed points that I want to route between using pgrouting. In order to do this accurately, I want to create new nodes (vertices) at nearest points on nearest edges to the actual coordinates to be used as sources and targets for routing. I also want to replace the original edge with the two edge fragments split from the original and have one of their endpoints be the new node and the other endpoint be the original endpoints of the original edge. Since the points are fixed, creating a fixed number of these new nodes permanently to the database should not be an issue.

What would be the best way to achieve this? It looks like the ST_NewEdgesSplit function [1] might be able to do this but it requires a parameter "varchar atopology", which I assume is a topology schema name. I'm not sure I have one, since I've constucted the databse with osm2po.


texts: Why does QGIS convert IntegerLists into something strange and how do I process these in PyQGIS?

texts:
I have GeoJSON data which has an integer list and I'd like to query that list in QGIS via GUI (select by expression, visualization) and probably later with PyQGIS.

The problem is that QGIS recognizes the field as an "IntergerList" and returns it as a "QString". That wouldn't be much of the problem but QGIS somehow transforms the original value ('[13, 35]') to something very strange ('(2:13,35)'). So QGIS somehow recognizes a List, but instead of returning a String that could be easily converted back to an array via string_to_array you have to process the data to be able to get the values, e.g.:

string_to_array(substr("long_list", strpos("long_list", ':')+1, -1))

So why is QGIS doing this, what are the advantages?

And is there an easy way to process the data in PyQGIS?

Here is the example data:

{
  "type": "FeatureCollection",
  "features": [
    {
      "geometry": {
        "type": "Point",
        "coordinates": [
          12.0,
          52.0
        ]
      },
      "properties": {
        "long_list": [
          13,
          35
        ]
      },
      "type": "Feature"
    }
  ]
}


texts:
I do have a point layer (data from a GeoJSON) and a polygon layer (also GeoJSON) in Leaflet. Both layers are created the same way and I made both layers "non-interactive". Unfortunately the point layer still reacts on the click-event on the point. The polygon layer is working just fine. Is this a bug or a feature?

Both layers are defined like this:

var layer = L.geoJson(json_data, {
            interactive: false
}).addTo(map);

My goal is to set the points to "non-clickable" aka non-interactive as well.
A working example is given here: https://codepen.io/butcher/pen/dBpWRP


texts: How to calculate azimuth in QGIS field calculator?

texts:
Is there any way to calculate the azimuth of two points in the QGIS field calculator? The coordinates are stored in an attribute table.

x1 | y1 | x2 | y2


texts: Why the Select By Location analysis calculates distance in a wrong way?

texts:
I want to select cities, which are within a distance of 200 km from my pre-selected city. However my result depends on whether I'm creating a buffer and use it to clip the cities (good result) or use the Select by Location analysis and supply it with the same data. My ArcGIS version is 10.2.1.

I figured out that the problem is that the analysis doesn't consider the data in a sphere (geographical) way, but in a planar (geometry) way. What I mean by that is, it calculates distance as in XY coordinates, so that the supposed buffer representing the distance has a perfect disk shape, which in Y direction is OK, however in X the results are tottaly wrong), meanwhile my data is geographical one and it uses deegrees (the data is in WGS84).

I made an experiment and created a buffer in PostGIS on the data imported as geometry, not geography and the result corresponds perfectly to the distance Select By Location uses. I've also created a buffer and used it to clip the cities and the result is correct.

I'm attaching the following files:

cities (a layer with some alaskan cities)
sel_city (a layer with a pre-selected city for the buffer)
city_buff200 (correct buffer create in ArcMap)
buff200_exp (buffer created in postgis, when data type was set as geometry and not geography) [this is a representation of the wrong distance buffer used by select by location]
cit_within200 (a layer with cities extracted using an incorrect distance representation)
cities_buff200_inter2 (correct result)

Data link: http://student.agh.edu.pl/~wysek/stack_data.rar

The code which I used to create the incorrect result is following:

cities_lyr = arcpy.mapping.Layer(cities)
arcpy.SelectLayerByLocation_management(cities_lyr, "WITHIN_A_DISTANCE", selected_city, "200 kilometers", "NEW_SELECTION")
arcpy.CopyFeatures_management(cities_lyr, "cit_withi200")


texts:
Does anyone have a download link for the NDVI colour scheme to use in ArcGIS for Desktop?

I am looking for MODIS mod13 data and cannot find an appropriate colour scheme in the programme.


texts: Interpolating LineString by adding new Point in middle between each two points

texts: I want to add interpolated points in a LineString using shapely. The interpolated points should be in the middle of every two Points present in the LineString. The shapely interpolate function adds the Points for the specified distance.


texts: How can I automatically deselect a mapped layer in Google Earth Engine?

texts:
I am currently mapping land cover changes between two years, e.g. 2001 and 2010.
The code I am writing adds the Land Cover layers to the map for both years, so a map of 2001 and one of 2010. Then, I add a layer containing exclusively the Land Cover Changes for a specific LC, e.g. cropland. I would like my code, when the LC change layer loads, to automatically deselect the 2010 LC layer so that the LCC layer 'rests' on the 2001 map rather than the 2010 one.
// --- SELECT STUDY AREA ---

var studyArea = ee.Geometry.Rectangle([-99.55454883653354,19.126233922295793, -98.93794117051792,19.612392241681285]);
Map.addLayer(studyArea, {}, 'studyArea');
Map.centerObject(studyArea);

var baseYear = 2001;
var studyYear = 2010;

// --- GET IMAGES ---

var landCoverStudy = ee.ImageCollection('MODIS/006/MCD12Q1')
                  .filter(ee.Filter.date(studyYear + '-01-01',  studyYear + '-12-31'))
                  .select('LC_Type1').reduce('mode').clip(studyArea);

var landCoverBase = ee.ImageCollection('MODIS/006/MCD12Q1')
                  .filter(ee.Filter.date(baseYear + '-01-01',  baseYear + '-12-31'))
                  .select('LC_Type1').reduce('mode').clip(studyArea);

var modisStyling = {
    min: 1,
    max: 17,
    palette: ['05450a','086a10','54a708','78d203','9900','c6b044',
              'dcd159','dade48','fbff13','b6ff05','27ff87','c24f44',
              'a5a5a5','ff6d4c','69fff8','f9ffa4','1c0dff']
};

Map.addLayer(landCoverBase, modisStyling, 'Land Cover (' + baseYear +')');
Map.addLayer(landCoverStudy, modisStyling, 'Land Cover (' + studyYear +')');

// --- CALCULATE DIFFERENCES ---

var selectedLandCover = 12; // 12 is cropland

var slcBaseYear = landCoverBase.remap([selectedLandCover], [1], null, 'LC_Type1_mode');
var slcStudyYear = landCoverStudy.remap([selectedLandCover], [1], null, 'LC_Type1_mode');

var remapedLcBaseYear = landCoverBase.remap([selectedLandCover],[2], 0, 'LC_Type1_mode');
var remapedLcStudyYear = landCoverStudy.remap([selectedLandCover],[1], 0, 'LC_Type1_mode');

var lcLoss = remapedLcBaseYear.subtract(remapedLcStudyYear).remap([2],[1],null);

Map.addLayer(lcLoss, {palette:['Black']}, 'Land Cover Loss');

In the screenshot below I have deselected the 2010 layer, but I had to do it manually.

Instead I would like to be able to code a command that deselects the layer without manual intervention.


texts: How to find if a lat/lon point is on a street?

texts: Are there any services to determine whether or not a give lat/lon point is on a street?  I am working on an app that uses user entered lat/lon points but it only really makes sense if its on a street.  I want to generate a bunch of points for testing but the points need to be on a street.  Any input is appreciated!


texts: I am looking to programatically generate an image from a WKT string (or any other common format, such as WKB or SQL Spatial geography). Are there any existing libraries that allow for this? I am working in .Net, but I am open to anything. The idea is that I could generate dynamic tile layers with the images instead of having to load up my Silverlight UI with heavy vector layers. Any ideas?


texts:
I have a SpatialLineDataFrame with 20 observations of 4 variables (in WGS84). I want to create waypoints every 10 meters. Every waypoint will has the x y coordinates. I tried to do it manually like:

point_first=SpatialPoints(cbind(x,y))

for the first point, but there are many points I want to do, so it is not possible to do it that way. I think I need a loop but I don't know how to do it in order to have the coordinates for every waypoint.

Is there any command that I should use?


texts: QGIS live GPS data via internal GPS (Durabook, Windows 10)

texts:
I have a Tablet (Darebook with Windows 10) on which I installed QGIS (v. 2.8)  So I tried to do the "live GPS tracking" with internal GPS (Sierra wireless sensor).

Next, from the menu of QGIS I went to "GPS information" and I clicked on the "Autodiscover" but I always get the message "unable to connect to the GPS." So I tried with the other options: internal, serial device, gpsd, but always the same message.

What am I doing wrong? What steps do I need to be able to do live tracking?


texts: ApplySymbologyFromLayer_management Overwrites Graduated-Color Class Break

texts:
My goal is to generate any number of MXD files with custom class-break-values using a consistent, pre-defined style.

For example, in ArcMap (10.2.2 and 10.3.1), I have an MXD with a graduated-color layer that has a single class-break-value symbol in it.  I want to update this layer's symbology from ArcPy using a .lyr file without affecting that class-break-value.

I can do this very easily from ArcMap: TOC > Layer Properties > Symbology > Import > Select the .lyr file > "Just the symbols" option > OK

Is there a way to mimick this from ArcPy?

First Attempt

When I try the following code (using the same .lyr file as above), the TOC layer's style is updated correctly, but its class-break-value is overwritten by the class-break-value from the .lyr file:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
dataFrame = arcpy.mapping.ListDataFrames(mxd)[0]
sourceLayer = arcpy.mapping.Layer(pathToLyrFile)
destLayer = arcpy.mapping.ListLayers(mxd, existingLayerName)[0]
arcpy.ApplySymbologyFromLayer_management(destLayer, sourceLayer)
arcpy.RefreshTOC()

Second Attempt

UpdateLayer has a symbol-only option, but it also overwrites the class-break-value:

arcpy.mapping.UpdateLayer(dataFrame, destLayer, sourceLayer, True)

Third Attempt

Originally, I tried importing the .lyr file into a blank MXD file after modifying the class-break-value, but when I change the value this way, the class-break's symbol is lost (and we are back to where we started):

import arcpy
mxd = arcpy.mapping.MapDocument(mxdPath)
dataFrame = arcpy.mapping.ListDataFrames(mxd)[0]
layer = arcpy.mapping.Layer(pathToLyrFile)
layer.symbology.classBreakValues = [0, 50]
arcpy.mapping.AddLayer(dataFrame, layer, 'TOP')
mxd.save()


texts: Saving styles for all layers using PyQGIS - SLD, QML, databases

texts:
I wrote a script (for QGIS 3) that saves styles for all visible vector layers in QML and SLD files:

import processing
self=qgis.utils
layers = self.iface.mapCanvas().layers()

for layer in layers:
    layerType = layer.type()
    if layerType == QgsMapLayer.VectorLayer:
        name = layer.name()
        pathqml = 'C:/Style/QML/'+str(name)+'.qml'
        pathsld = 'C:/Style/SLD/'+str(name)+'.sld'
        layer.saveNamedStyle(pathqml)
        layer.saveSldStyle(pathsld)

How can I do this for the layers from SQLite, PostgreSQL/PostGIS ?


texts: gdalwarp producing empty rasters when called from python script, but not via command line

texts:
I'm trying to aggregate some raster data (all .tif) up to a coarser resolution (from .05 degree to .25 degree) using gdalwarp in python, but the command is not working. Instead of getting an output with a wide range of values, all values of the output are 0. The resolution and pixel depth/type are correct, but the values are not.

Here is the documentation to the gdalwarp command: http://www.gdal.org/gdalwarp.html

I have two input files which I wish to aggregate up to .25 degree resolution, producing several outputs:

'NDVI_raster': The first input is a 16-bit signed raster representing NDVI, with values ranging from -10,000 to 10,000 and nodata values of -15,000.
'nodata_mask': The second is a NoData mask, 32-bit float, where 1 = values of "good" data, and 0 = NoData.

With 'NDVI_raster' as input, I want to produce 7 different outputs, each representing a different statistic. I do this by calling gdalwarp 7 times, each time setting the resampling method (-r) to one of the following: average, mode, max, min, median, q1, q2. I'll call the outputs NDVI_ave.tif, NDVI_mode.tif, etc. Right now, I am using GDAL 1.10.1, which only allows average and mode, so I'm testing with just these two stats now.

Using 'nodata_mask' as input, I want to ultimately produce a QAL (quality assurance layer). To do this, I use gdalwarp, with the resampling mode set to 'average' to aggregate up to .25 degree. This results in each pixel representing the ratio of good pixels/total pixels from the input. Let's call the output QAL.

Here is what's in my code (using mode as example for the first input):

os.system('gdalwarp -tr .25 .25 -r mode -srcnodata -15000 %s %s' % (NDVI_raster, NDVI_mode))

And for the QA layer:

os.system('gdalwarp -tr .25 .25 -r average -srcnodata -15000 %s %s' % (nodata_mask, QAL))

The results are rasters with the correct resolution, projection, and pixel depth, but the pixel values are all 0.

Anyone familiar with python/gdal know what's going on?

When calling the gdalwarp command from the command line (linux), I get the desired result. When using os.system to call gdalwarp from python, I get empty rasters. So maybe something is wrong with my gdal/python bindings?

Instead of calling the command via os.system, I used subprocess. The tool via this method also appears to run smoothly, but the result is the same: a raster full of 0's.

I tried putting the gdalwarp call in a bash shell script and calling that shell script from python, but the result is a bunch of -1s instead of 0s. Weirdly enough, I had tested it before and am pretty sure it worked, but the test got deleted off my server and now I can't recreate it for some reason.

Putting the gdalwarp command in a bash shell script and then calling that shell script from the command line gives me the desired result. Calling the same shell script from python, though, doesn't. It's looking like something is wrong with python, but what and how do I fix it?


texts: Creating dropdown option when setting up python script into arctoolbox?

texts:
I created a script that I will import to a arctoolbox. On the display name section, I have asked the user to choose a workspace, select the shapefile, field name, and an option to make a pdf, jpeg or none.

For the last option, instead of making the user write whether they want pdf, jpeg, or none, is there an option where I can create a dropdown menu with the 3 selection?


texts: "Copy" Image coordinates to another image that is nd.array

texts:
I have image with three bands. I have script predicts the first band pixels values by using the two  other bands and in the end generates an image.The generated image is at the beginning pandas table, that transformed into nd.array and then displayed and saves as tif image using imageio.
My problem is that during this processing I lose the coordinates so the result image needs to be georeferencing.
During the process I use reshape in order to get one table of all the pixels and their 3 values and I believe this is where I lose the coordinates but I don't know how can I keep them , if they should become a column in the new table? or they stores somehow different?
My script:
#open the raster I have download before
img=rasterio.open("img_new.tif")
show(img,0)

#create pandas df where each pixel is a row, the column are the bands
#probably here i'm losing the coordinates
df_all=pd.DataFrame(array.reshape([3,-1]).T)
df_all

#use random forest regressor to predict the first band by bands 2,3:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
rf=RandomForestRegressor()
y_train=y_train.values.ravel()
rf.fit(X_train,y_train)
rf_pred=rf.predict(X_test)
rf_pred

#apply rpediction for all the data
pred_all=rf.predict(data)

#.....creating new df table with the prediction value for all....
df_join=df_all.merge(df,how='left',left_on='index', right_on='index')

#convert back to image that doesn't have the coordinates:

#convert to numpy
rf_array=df_join['Prediction'].values
rf_array

#reshape
rf_array=rf_array.reshape(869,1202)

plt.imshow(rf_array)

I thought maybe if I could tell python to choose some pixels from the first image, save their value and paste it later somehow in the generates result image it could work .
My end goal: to "copy" the coordinates from the first image to the result image , so when I open the two images in QGIS they will overlap.
Edit: just to clarify: the result image is constructed from numpy.ndarray that I save with imageio as  image.
Edit2: I have c found the img.transform and could get it from the original image:
img.transform
>>>Affine(10.0, 0.0, 208810.0,
       0.0, -10.0, 7583530.0)

but now I don't know how to get this coordinates to be pasted on the result image.
Edit3: definition of X, y, df, array:
#Definition of df
#df is the pandas dataframe, constructued from the original tiff:
img=rasterio.open("image_original.tif")

#array
#shape
array=img.read()

#create pandas df
df=pd.DataFrame(array.reshape([3,-1]).T)
df.columns=['band1','band2','band3']
df=df.reset_index()
df

#define X and y, y is the predicted values (I wanted to rpedict y using columns X #with Random Forest)
X = df.iloc[:, 2:]
y = df.iloc[:,1:2]


texts:
I have a shapefile that contains three lines. I want to add a specific numeric id to the field "line_id" and a string id to the field "border_id".
Basically, for the first line, I want to add the value '1' to "line_id" and 'North' to field "border_id". For the second line I want to add the value '2' to "line_id" and 'South' to field "border_id".

Here is my code, I think I don't get how the changeAttributeValues command works on pyQGIS.

lyr.startEditing()
provider = lyr.dataProvider()
provider.addAttributes([QgsField("line_id",QVariant.Double),(QgsField("border_id",QVariant.String))])
lyr.updateFields()
lyr.commitChanges()


texts: How to -explode only one subfeature from only one feature?

texts:
I have a topojson which treats French Guiana as part of the "France" feature (therefore also, part of the Western Europe layer). To match my dataset, I want to split this out from France, while keeping Corsica as part of France:

On browsing the Mapshaper docs, -explode looks like my best bet. But there are a few problems:

It explodes out every multi-part feature in a layer, and there doesn't seem to be any way to limit it to just one feature or field value: target= is layers only. This isn't a huge problem for me in this particular example, as for this data the other features on this layer aren't multi part and I could temporarily move it onto a dummy layer, but I'll also need to do this with a variant where there are other multi-part features on the layer and the dummy layer workaround is rather clumsy (I'd like to preserve the order of the features in the layers if possible), so if there is a way to specify the target feature that would be useful.
It explodes all the parts of the feature, and there are two problems here:

The obvious way to recombine them, -dissolve or -dissolve2, in my experience seems to require the features to be touching or overlapping.
If I can get around that, another difficulty with -dissolve is that I would need to specify Corsica but not French Guiana using field data, but after exploding they get identical field data to the original feature (France). This technique using MapShaper's internal feature ID is a possibility here, not tried it in this context yet

So:

Is there a simpler, cleaner workflow for this general case than moving "France" to a dummy layer, exploding the dummy layer, then trying to dissolve Corsica back into France?
If no, how can I reliably dissolve Corsica back into France when they don't touch or overlap?

Please, no puns like "Can Mapshaper do this? Of Cors-i-can!"...


texts: How can I fix misshaped GeoJSON polygons?

texts:
In 2018 Ghana changed its administrative divisions (from 11 to 16 regions, and from 216 to 260 districts). Unfortunately, I have not been able to find the GeoJSON (or equivalent) data for the new divisions. I managed to scrape map data from a website I will not name, by looking into AJAX requests. Except the coordinates of the regions were using the Web Mercator system (ESPG:3857).
I coded a small script to convert them into Lon/Lat (ESPG:4326) and build a GeoJSON file, and here is the result :

Some polygons appear to overlap each other... I believe this may be caused by the massive lake in the center of the country (and also the smaller one in the south east).
I thought maybe my conversion was off so I tried generating my GeoJSON without converting the coordinates, and the result is the same.
What is wrong with my data? And how can I fix it?
Is there a script I can use to remove points that are outliers?


texts: How can I add points to a LineString in shapely?

texts:
I have two Polygons which intersect along one edge as shown below.

I need to add the red vertices at each end of the intersection to the blue polygon. I thought I could do this with by calling union on the Polygon and the LineString produced by intersection of the two polygons but it appears not.

Is there a way of doing what I need to do?


texts: Free Digital Elevation (Terrain) Model of Europe or Germany

texts: Is there any free (or not expensive) source of Digital elevation model of Europe or Germany with the resolution of 10 meters or better and vertical accuracy of +-5m or better?


texts: Different pixel values in ArcGIS Desktop and R?

texts:
I am using the NAFD _NACP open source data set available at

https://daac.ornl.gov/NACP/guides/NAFD-NEX_Forest_Disturbance.html

As noted in the metadata, a raster for 2010 has 6 classes (Data gaps, Open Water, Forest Cover, no Forest Cover etc). When I open this raster in Arcmap, I see that only these 6 values are present. Also the attribute table has only these values.

However, when I open this raster in R (reading it using raster() in raster package), the min and max values are 0 and 255.

Why is there this difference between ArcGIS Desktop and R?

Also, if I have to reclassify this raster to keep only pixel values=3 (Forest Cover), should I recalssify everything below and above 3 as NA using recalssify(nafd 2010 raster,c(0,3,NA,3,255,NA). I tried this and when I export the result, the raster becomes a float 32 bit unsigned type (the original is 8 bit unsigned).

Why is this happening?

I have tried to specify in the writeRaster() that the dataType=INT1U, but this does not do anything.


texts: QgsRasterFileWriter.writeRaster() is crashing python program

texts:
I am doing ndvi analysis on a farm image. After ndvi I need to convert the image into Singleband pseudocolor image using following python code.

def calculate():

    file_name =  CUR_DIR + "\sampleodm.tif"
    out_file = CUR_DIR + "\sample_output.tif"
    out_file_ndvi = CUR_DIR + "\out_ndvi.tif"
    try:

        file_info = QFileInfo(file_name)
        rasterName = file_info.baseName()
        raster = QgsRasterLayer(file_name, rasterName)

        if not raster.isValid(): print 'Invalid raster'

        ir = QgsRasterCalculatorEntry()
        r = QgsRasterCalculatorEntry()

        print ir, r

        ir.raster = raster
        r.raster = raster

        ir.bandNumber = 4
        r.bandNumber = 3

        ir.ref = rasterName + "@4"
        r.ref = rasterName + "@3"

        references = (ir.ref, r.ref, ir.ref, r.ref)
        exp = "1.0 * (%s - %s) / (%s + %s)" % references
        print exp
        print 'here2'

        output = out_file

        e = raster.extent()
        w = raster.width()
        h = raster.height()
        entries = [ir, r]

        ndvi = QgsRasterCalculator(exp, output, "GTiff", e, w, h, entries)

        ndvi_res = ndvi.processCalculation()
        print ndvi_res

        file_info_out = QFileInfo(output)
        rasterName_out = file_info_out.baseName()
        lyr = QgsRasterLayer(output, rasterName_out)
        if lyr.isValid(): print 'valid raster'

        min = 0.0341278
        max = 0.56844
        range = max-min

        algorithm = QgsContrastEnhancement.StretchToMinimumMaximum
        limits = QgsRaster.ContrastEnhancementMinMax
        lyr.setContrastEnhancement(algorithm, limits)

        i = []
        i.append(QgsColorRampShader.ColorRampItem(0, QColor(0,0,0,0), 'NODATA'))
        i.append(QgsColorRampShader.ColorRampItem(0.0341278, QColor(120,69,25,255), 'Lowest Biomass'))
        i.append(QgsColorRampShader.ColorRampItem(0.168, QColor(255,178,74,255), 'Lower Biomass'))
        i.append(QgsColorRampShader.ColorRampItem(0.244, QColor(255,237,166,255), 'Low Biomass'))
        i.append(QgsColorRampShader.ColorRampItem(0.301, QColor(173,232,94,255), 'Moderate Biomass'))
        i.append(QgsColorRampShader.ColorRampItem(0.435, QColor(135,181,64,255), 'High Biomass'))
        i.append(QgsColorRampShader.ColorRampItem(0.494, QColor(3,156,0,255), 'Higher Biomass'))
        i.append(QgsColorRampShader.ColorRampItem(0.56844, QColor(1,100,0,255), 'Highest Biomass'))

        c = QgsColorRampShader()
        c.setColorRampType(QgsColorRampShader.INTERPOLATED)
        c.setColorRampItemList(i)

        s = QgsRasterShader()
        s.setRasterShaderFunction(c)
        print lyr.rasterType()
        ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), lyr.rasterType(), s)
        lyr.setRenderer(ps)

        QgsMapLayerRegistry.instance().addMapLayer(lyr)

        extent = lyr.extent()
        width, height = lyr.width(), lyr.height()
        print width, height
        renderer = lyr.renderer()
        print renderer.type()
        provider = lyr.dataProvider()
        crs = lyr.crs()

        pipe = QgsRasterPipe()
        pipe.set(provider.clone())
        pipe.set(renderer.clone())

        file_writer = QgsRasterFileWriter(out_file_ndvi)
        print file_writer
        print 'write'
        # upto here everything is working fine and program is generating black and white image.
        error = file_writer.writeRaster(pipe, width, height, provider.extent(), crs)
        # print error
        print 'Completed NDVI analysis!'

    except Exception as e:
        print 'Exception: ', e

Consider every module is imported in program.
line error = file_writer.writeRaster(pipe, width, height, provider.extent(), crs) is causing Python.exe has stopped error.

I tried this program on QGIS Desktop application's python console and it is saving the singleband pseudocolor image perfectly.

How do I rectify the error?


texts: Styling coordinates in QGIS's Map Composer

texts: Is it possible to achieve the following style for coordinates in a map grid in Map Composer?


texts:
I have downloaded landsat image from this site and now Im wondering if i can get NDVI index from it.

In raster calculator I have only 3 ranges, shouldnt be four?

Or am I doing something wrong?

If it helps its Landsat 8 sensor: OLI


texts: Performing spatial join on existing feature class?

texts:
I have feature classes containing points (Events) and polygons (Sectors). I want to add a new field to the Events FC to store a field from the Sectors FC.

The Events FC would be the Target Features and the Sectors would be the Join Features, using the terminology of a Spatial Join (Toolbox > Analysis > Overlay > Spatial Join). I need to get information for points from the polygons in which the events occurred.

I could perform a Spatial Join but that results in a third FC. I prefer to add a new field to the existing Events FC. How do I perform a spatial join and add a field from the resulting join to an existing FC?


texts: Problem with loading OSM planet into PostgreSQL

texts:
I try download OSM planet into PostgreSQL database.

I do all steps from instruction http://switch2osm.org/serving-tiles/manually-building-a-tile-server-12-04/

I also modify string "kernel.shmmax"

Now I try download:

osm2pgsql --slim -d gis -C 800 --number-processes 4 planet-latest.osm.pbf

But I have error:

Reading in file: planet-latest.osm.pbf
Processing: Node(2041920k 189.9k/s) Way(199297k 3.39k/s) Relation(47930 14.42/s)osm2pgsql: PolygonBuilder.cpp:261: geos::geomgraph::EdgeRing* geos::operation::overlay::PolygonBuilder::findShell(std::vector<geos::operation::overlay::MinimalEdgeRing*>*): Assertion `shellCount <= 1' failed.

Emergency stop(made a memory dump)

How can I solve this problem?

UPDATE

It is important to know, that the problem is in data. I use planet-latest.osm.pbf     10-Oct-2013 13:29   22G

When I download rostov.osm.bz2 there are no problem.

Which version of libgeos is usefull for me?

UPDATE 2

Well. It is interesting to know, that I try to download planet-131002.osm.pbf file. I verify md5. It is equal.

UPDATE 3

I instal the latest version of libgeos. It is 3.3.3-1. Then I try download planet-131002.osm.pbf, and then planet-130614.osm.pb

But I have the same problem.


texts: How to control output of Jiffle Processes in GeoTools so to return float 32bit instead o double 64bit

texts:
While developing custom map algebra in GeoTools, I adopted Jiffle.
(see : Creating custom Operation for encapsulating some map algebra using GeoTools?)

This allowed very easy implementation of logic but two major issued raised:

Jiffle as stated here does all calculations using double values
Dumping to GeoTIFF the output ++can take a lot of time++: more than one would expect.

This is how I dump a GridGoverage2D from a Jiffle process.

[.. jiffleProcess defined ouside this snippet ..]
Map<String, Object> output = jiffleProcess.execute(inputs, null);
GridCoverage2D ndviCoverage = (GridCoverage2D) output.get(JiffleProcess.OUT_RESULT);
File outputNDVIGeoTIFF = new File("NDVI.tif");
GridCoverageWriter writerNDVI = geotiffFormat.getWriter(outputNDVIGeoTIFF);
writerNDVI.write((GridCoverage) ndviCoverage, null);

So the questions are:

Can file saving time be shortened in any way? (I mean from GeoTools / Jiffle perspective, not hardware one)

.. but more important:

Is there any way to have Jiffle output GridCoverage2D in float to reduce the bits required to store data ?


texts: Shifting on the X-axis while reprojecting Gauss-Krüger zone 3 to UTM zone 33N

texts:
I'm having trouble converting geodata from Gauss-Krüger Zone 3 to UTM 33N.
See attached image.

I don't know why, but there seems to be a shift on the X-axis. I found a similar  question here but am not sure how to do this in QGIS:
https://sourceforge.net/p/saga-gis/discussion/790705/thread/a824273a/?limit=25#567a

The left map is from here: https://gadm.org/download_country_v3.html

The right map is from here: https://opendata.dwd.de/climate_environment/CDC/grids_germany/multi_annual/precipitation/


texts: Forcing ArcMap Add-in ComboBox to refresh using ArcObjects?

texts:
I am having problems with an add-in combobox that doesn't refresh.  When the currently selected item is removed from items collection, the name remains visible in the text area.

It should show the combobox's caption, but that doesn't happen until the mouse is moved over it.

Is there a way to force it to refresh programmatically ?


texts: Using ArcGIS Online for public-facing site, while keeping source datasets private?

texts:
I would like to build a public-facing map using the ArcGIS Server JavaScript API, which will display datasets stored in ArcGIS Online.

ArcGIS Online seems to work on the premise that all users have a Named User account, and will be signing in to their Named User account before accessing any secure content from ArcGIS Online.

This works well for an organisation, but is less suitable for a public-facing map which doesn't require users to log in.

If the map is to be accessible without requiring the users to log in, then the datasets must be shared with Everyone on ArcGIS Online, which means that any ArcGIS Online user has access to the source datasets.

Is there an option for creating a public-facing map which doesn't require users to log in to ArcGIS Online, without requiring that the source datasets be shared with Everyone in ArcGIS Online?


texts: Dropdown list and map.fitBounds in Leaflet

texts:
In my Lealfet application I have a dropdown list which uses data from a GeoJSON to directly populate the list.
HTML

Seleziona il Foglio

  <select id="fogli-dropdown" name="foglio">
  </select>

JAVASCRIPT
let dropdown = $('#fogli-dropdown');
dropdown.empty();
dropdown.append('<option selected="true" disabled>Seleziona il Foglio</option>');
dropdown.prop('selectedIndex', 0);
var fogli = new L.geoJson(fogli, {
    onEachFeature: function (feature, layer) {
        dropdown.append($('<option></option>').attr('value', layer.feature.properties.fg).text(layer.feature.properties.fg));
    }
});

According to the select item in the list (fogli) I would like that my map fits the current extent of the selected item.


texts: Trying to extract arc information from a Suburb shapefile

texts:
I downloaded some shapefiles from the Australian Bureau of Statistics (ABS). It's basically showing a bunch of suburbs when I view the shapefile through a GIS tool. The GIS tool that I am using at present is the GRASS 6.4.2 tool.

So basically the shapefile has the below statistics

Number of points:       0               Number of areas:      15231
Number of lines:        0               Number of islands:    7003
Number of boundaries:   31362           Number of faces:      0
Number of centroids:    15220           Number of kernels:    0

I'm trying to extract from this shape file the arcs. That is the polylines of the boundaries. This is because I need to find for each arc its accompanying adjacent boundary.

For example arc1 may have boundary with Cat id - 135 on its left and Cat ID = 125 on its right.

Any clues as to how to do this would be appreciated.


texts: Transforming azimuthal to wgs projection using Python?

texts:
I want to plot a circle around a given lat/long with a certain radius on the Earth. The shapefile of the earth has WGS84 coordinate system.

First, I projected the circle to azimuthal equidistant map to apply the real distance on the map and transform it to WGS84. Then, I plotted is as follows:

(with lat=34.2255804,
lon=139.294774527, radius=13006 km)

aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84', lat_0=lat, long_0=lon)

wgs = pyproj.Proj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

circle = transform(functools.partial(pyproj.transform, aeqd, wgs, Disk(0, 0, radius*1000))

But, when I compared the result with google-radius, the result seemed reverse of my plot as follow:

Am I missing something when I project or transform the circles? Can anyone explain why it is happened?


texts:
I have some data that encodes locations in certain countries, such as the US and China, that I would like to map. The data includes location name in Latin, extended Latin and other scripts (Japanese, Thai etc), location latitude, location longitude and some categories.

Ultimately I would like to be able to show this data on a map with administrative boundaries (roughly state/province level) and/or conurbation overlay. Likely output would be colour PNG images for use on the web.

I thought about using R (my main tool for data manipulation) to do a trial run with some of the US data but immediately stepped into a welter of jargon concerning shapefiles, attributes and the problem of trying to make a decision at least 3 different GIS packages for R. There are so many tools and moving parts out there that I am having difficulty seeing the forest for the trees.

Q. Can anybody provide a brief step-by-step summary (or a minimal example) of how I should go about mapping such data using R and some qualitative feedback on which R tools would be the most suited to me needs?

An example of my data, expressed in R code, would look something like this:

locname <- ("Shanghai Stock Exchange")
loclocalname <- c("上海证券交易所")
lat <- c(31.236259)
long <- c(121.511299)

mydf <- data.frame(locname, loclocalname, lat, long)
mydf$latlong <- paste(mydf$lat, mydf$long, sep = ",")
str(mydf)


texts: Creating multiple hyperlinks in attribute table on relative paths

texts:
I got the task to create a layer with multiple hyperlinks within the attribute table. One link should go to the PDF document, that describes the feature further. The other link should go to one, or multiple pictures of the same feature. I would like to have only one field for all pictures.

Is it possible to create one hyperlink to multiple pictures? I will have to create hyperlinks for around 500 features, so the Option to create hyperlinks manually is not suitable.

My second problem is the fact, that I have to use relative paths for the hyperlinks.

I'm using ArcGIS 10.6 and QGIS 3.4 so if there is a solution within QGIS with a product, that works in ArcGIS as well, I would be happy with that too.


texts:
I am trying to extract raster values within multiple polygons. It would be nice if I could choose between say the maximum or an average raster value within the polygon.

For example, I have created depth grids as my raster to give the flood depth. Now I would like to calculate the highest/average raster value within each building footprint to estimate flood risk.

I have tried zonal statistics as recommended in a few similar posts, but the results are calculating uniform for every polygon instead of each individual building. I also have done polygon to point and extracted the rater value to the point, but in some cases the depth grid does not cross the center of the building footprint where the point is placed.


texts: Calling `Rscript` from QGIS Python Console / QGIS Plugin not working

texts:
I want to run a RScript that plots some data which is fetched by a QGIS Plugin. However it returns a non zero-exist status and breaks (when running with subprocess.check_call) or just skips the command when calling subprocess.call / subprocess.Popen.

Things I've done so far:

check_call:

script_filename = "/myRscript/that/plots/some/stuff.R"
csvpath = "/path/of/csvdata.csv"
output_dir = "/output/directory"
p1 = subprocess.check_call(['/Library/Frameworks/R.framework/Versions/3.5/Resources/Rscript', script_filename, "-f", csvpath, "-o", output_dir], shell=False)

returns following Error:

>>> subprocess.CalledProcessError: Command '['/Library/Frameworks/R.framework/Versions/3.5/Resources/Rscript', '(...)createProfile.R', '-f', '(...)profile_data.csv', '-o', '(...)4 MASTERARBEIT']' returned non-zero exit status 1.

Popen:

p1 = subprocess.Popen(['/Library/Frameworks/R.framework/Versions/3.5/Resources/Rscript', script_filename, "-f", csvpath, "-o", output_dir], shell=False)
p1.wait()
>>> 1

However, the calls work fine when I try to call them from a standard python session outside of QGIS. What is even more irritating is that a command like subprocess.check_call(['open', output_dir], shell=False)works inside QGIS.

BTW, I'm on Mac OS X High Sierra and QGIS 3.2.0 Bonn / Python 3.6.


texts: QGIS TimeManager with multiple features at one time?

texts:
Since the data I have to deal with is really large, it is not practical to separate every feature into one line in the CSV file.

Is there any way to animate them as a whole layer?

I have approximately 26000 features in one layer for each time interval, and 300 time intervals that i need to implement the animation.

EDIT:

I want to display links and animate their colour changes based on the magnitudes.

Basically, all links should show up at all time, but since their magnitudes change with time, and I colour coded the style, the colour of the links will change as time passes.

My previous data format is: magnitude; time; coordinates in WKT format. one example is: 4; 2015-01-01; (642044 4862161, 642319 4862144). However, in this form, if I need to include all the data, it requires 7,000,000 lines which exceeds the maximum row limit for excel. So now I want to know if there is any chances to animate them as a whole layer, or some ways to solve this problem.

Also if I input all the data, the program probably becomes really slow.


texts: Can't see the Fields I'm trying to add in pyQgis?

texts:
On Qgis 2.14.11 I successfully made buffered envelopes from my "roads" Qgsvectorlayer, and they load fine on the screen.

I am now trying to feed a QgsFields object to a vector layer shapefile and I don't see that reflected in the data table.

Nor have I been able to copy the attributes from every feature in "Roads" to each new feature I put into "layers"

I suspect that my code is lacking...

 from PyQt4.QtCore import *
    import os

working_dir = '/example'
os.chdir(working_dir)
project = QgsProject.instance()

roads = iface.addVectorLayer("/example/ALL_zones_updated_dissolved.shp", "ZDOT_Roads", "ogr")
if not roads:
    print("Layer failed to load!")

temp_path = 'temp'
os.chdir(temp_path)

fields = QgsFields(roads.fields())

crs = QgsCoordinateReferenceSystem(4523, QgsCoordinateReferenceSystem.PostgisCrsId)
writer = QgsVectorFileWriter("multi_gon.shp", "utf-8", fields, QGis.WKBMultiPolygon, crs, "ESRI Shapefile")
if writer.hasError() != QgsVectorFileWriter.NoError:
    print "Error when creating shapefile: ",  w.errorMessage()

layer = iface.addVectorLayer("/examples/temp/multi_gon.shp", "Roadbounds", "ogr")
if not layer:
    print("Layer failed to load!")

for feature in roads.getFeatures():

    inAttr   = feature.attributes()
    inFields = roads.fields()
    inGeom = feature.geometry()
    bf_inGeom = inGeom.buffer(500,1).boundingBox()
    gPolygon = QgsGeometry.fromRect(bf_inGeom)

    feat = QgsFeature()
    feat.setGeometry(gPolygon)
    feat.setFields(QgsFields(roads.fields()))
    feat.setAttributes(inAttr)
    layer.dataProvider().addFeatures([feat])


texts: calculate average commute time within 15 miles of city center

texts:
I am trying to figure out how to calculate the average commute time to work within 15 miles of the city center (and then map it in a form of continuous symbols (red to green gradient) with a legend, scale and everything. That would be called “average commute times to work).

Could you help me for the first part of the problem? Thank you very much


texts: How to calculate weighted mean coordinates within polygons?

texts:
I'm using QGIS 1.8 on OSX 10.8.3.

I have a polygon shapefile of all of the US state legislative districts.

I also have a layer of points that are centroids of US Census Tracts, with a field for the population of the tract.

How do I generate a layer of points that are the weighted mean of each census tract point within each legislative district?

My hunch is to play with spatial joins and use fTools to do mean coordinates, but I'm not sure how to calculate the mean for each polygon instead of the whole layer itself.

Any thoughts?


texts: Anybody know any WMS clients for Android?

texts: I am interested in finding a project (with source) for a WMS client for Android.


texts: Torque heatmap on CartoDB with time control

texts:
I am interested in create a Torque Heatmap in which you can see the events occured every year in one step, for twelve years in total. Every event has a aggregation or ponderation value (they are the number of times that something occurs in a place, 771 places surveyed). I tried two ways:
- Unique 771 places with ponderation value per year. I used 11 steps on cartodb_id so i have all 771 places shown every year (step) but ponderation value is missing (I used Torque Aggregation Function as "count (value)", "sum (value)" and "value" itself, not working).
- Adding a date field like "12/31/year" and replicating each place as many times as its ponderation value for this year, i used again 11 steps on "date field" and Torque Aggregation Function as "count (cartodb_id)". Now i can see intensity in aggregated points, but no correspondence in dates because there are not the same number of events in every year.

Can anybody help me?


texts: ENVI Error running FLAASH Atmosphere Correction

texts:
I've made a Radiometric calibration on a Landsat-8 image.
Next, I've tried to make a FLAASH Atmosphere Correction from the Toolbox but I get a warning message and the process stops with errors.

The warning indicates that for western hemisphere the longitude should be negative and in fact, the Landsat-8 image I am processing has a positive longitude.

The warning message is:

How can I fix this problem?

The Input parameters I use for the FLAASH Atmosphere Correction are:

For the Multispectral settings I've used:

and...

and the finnal error message:


texts: Understanding arcpy.CalculateField_management with Join?

texts:
I am using arcpy.CalculateField_management to fill a column field with another column values.

I have two shapefiles, fishnet.shp and new_run.shp. First I joined both using arcpy.AddJoin_management as follows:

arcpy.MakeFeatureLayer_management("fishnet.shp","fishnet_layer")
arcpy.AddJoin_management("fishnet_layer","FID","new_run_grid.shp","FID_fishne")

Now I want to assign fishnet.shp's column rungrid the value of the joined table's column final_run value.

How to use  arcpy.CalculateField_management here?

I tried the following code but it gave error:

arcpy.CalculateField_management("fishnet_layer","rungrid","!new_run.final_run","PYTHON")

I am getting the following error:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000728: Field rungrid does not exist within table
Failed to execute (CalculateField).

Further

I checked the names of the fields of fishnet_layer and these two names are there fishnet.run_grid and new_run.final_run. After giving the following command:

arcpy.CalculateField_management("fishnet_layer",'fishnet_f.run_grid','new_run_grid.run_grid1',"PYTHON")

I am getting the following error:

ExecuteError: ERROR 000539: Error running expression: new_run_grid.run_grid1
Traceback (most recent call last):
File "<expression>", line 1, in <module>
NameError: name 'new_run_grid' is not defined
Failed to execute (CalculateField).


texts: How to import or incorporate formated Excel tab into QGIS print composer

texts:
I would like to know if there's a way to import formated Excel Tab into QGIS print composer.

Below is an attachment of formatted Excel tab example I has been asked to display on qgis map composer.

Regards!


texts: How to show shape stored in attribute table from .geodatabase in QGIS?

texts:
I have a geodatabase (filename ends with .geodatabase) and when I drag this into QGIS, I see all the tables and can browse them.

On some of those tables there attributes that look like they are the shapefile encoded or in a binary format or something. See this example

In the table properties, source fields, this field is defined as Type:QString TypeName: Binary.

Assuming I am correctly understanding this and it is the shapefile. How can I get QGIS to see this as a layer?


texts: Seeking Python script for creating .mxd files?

texts:
I am new to both ArcGIS and Python. My requirement is to automate the below MANUAL process:

Creating a layer in ArcGIS for Desktop. To put it in another words, creating an ArcMap document (.mxd).
Publishing the created ArcMap document (in Step 1) as a service to ArcGIS Server.

Currently we are doing this manually. I have come across scripts to automate step 2 using Python.

How can I automate step 1 and step 2?


texts: Mapproxy transformation from Web Mercator to WGS84

texts:
I want to transform tiles defined in Web Mercator (EPSG:3857) to WGS84 (EPSG:4326) using Mapproxy. Unfortunately the resulting WGS84 map gets more and more distorted to the poles. Here is my configuration

services:
  demo:
  tms:
    use_grid_names: true

layers:
  - name: topo
    title: Topographic map from ArcGIS Online
    sources: [topo_cache]

caches:
  topo_cache_source:
    grids: [GLOBAL_WEBMERCATOR]
    disable_storage: true
    sources: [topo_tiles]
    format: image/png
  topo_cache:
    grids: [GLOBAL_GEODETIC]
    sources: [topo_cache_source]
    format: image/jpeg

sources:
  topo_tiles:
    type: tile
    url: http://services.arcgisonline.com/arcgis/rest/services/World_Topo_Map/MapServer/tile/%(z)s/%(y)s/%(x)s
    grid: GLOBAL_WEBMERCATOR

Around the equator it looks fine but it gets more and more distorted closer to the poles:

As I zoom in, the distortion gets smaller and smaller. In Leaflet zoom level 5 and up it is hardly noticeable. At max zoom it matches up perfectly.

Does anyone have an idea about what is going on? I'm not sure if I have done something wrong in the configuration or if there is a bug in Mapproxy. I have tried using Open Street Map as source and get exactly the same result. I have also tried defining a custom grid with specified bounding box but nothing I do gets rid of the distortion.


texts: I installed PySAL but i can't find the toolbox for ArcGIS.


texts: Why do my out fields change what map points are displayed on my map?

texts:
I have a feature layer within my map control. When I do not list out fields on the feature layer control within the LayerCollection, I see the points, and they appear as expected as I navigate around the map.

BUT, when I add a single outfield that exists, the map points no longer appear! I know the layer can display them, because when I filter the layer for a specific value within the outfield attribute, The point appears.

when it appears though, navigating around the map will change if the point is displayed. for example, if I zoom in, where the point no longer is in the map, and pan to its position, it no longer is there. When I zoom out or in, it appears, but if the extent does not display the point when zooming, it disappears.

What the heck is going on?


texts: Properly setting Extents in Arcpy Using path

texts:
I am trying to set the extents of my processing using Arcpy, but when my processing is completed, there has been no change to the extents (it results in an extent that did not go as large as I desire).

Included is the code I have been using for the extents:

arcpy.env.extent = 'C:\data\myData.tif'

What I am using as an input would bea .tif file, so I am not sure if that would matter or not.

After using this, I am attempting a Euclidean Distance analysis, where I have:

arcpy.gp.EucDistance_sa(input, euclidean_output, "", "200", Output_direction_raster)

As I said earlier, it just returns the same output without the extent change. Would anyone be able to give me some tips?


texts: Returning 1 or 0 in column based on point location using QGIS?

texts:
Within QGIS - I have two layers. One is a point layer, the other is a polygon layer.

I'd like to create a column in the point's attribute table using the field calculator where 1 indicates the point is within polygon and 0 indicates the point is not within polygon.

Any ideas?


texts: Allow new tiles to be loaded while dragging map

texts:
At least by default (for example with the OL3 examples), new titles are not loaded while the map is being dragged. A blank area is shown where tiles have not loaded yet. When the dragging stops, the tiles are loaded and the blank areas fill in.

How can I configure Openlayers3 to allow new tiles to be loaded while dragging the map?


texts:
I've got 69 .asc files to make a nice colored and shadowed elevation map or levels contours, and maybe a 3D view.
They all open and get shadows. One by one…
So I would like to create a VRT ("menu:raster:miscellaneous:Build virtual raster" or with toolbox:GDAL. When I open the window to select my files and hit "ok", both VRT window and select box close and it does nothing.
When I try to merge several .asc file I get the same result.
What went wrong? are .asc not suitable for VRT? A Python crash? A badly configured QGIS?
Or just out of memory? (each asc file has 1 million points, it made a lot, but with just 4 tiles it does the same "nothing")
MacOS High Sierra (10.13), QGIS 3.10 Coruña, GRASS 7.8, GDAL/OGR 2.4.1, Python 2.7
The installation of GRASS was tricky (config of path to GRASS)
When I  try to create contour, it runs but log says that the result layer was not correctly generated.


texts: Missing v.in.ogr module of GRASS in QGIS Processing toolbox?

texts:
I am a new user to QGIS/GRASS and I am having a problem utilizing GRASS.

I have several files that I generated in QGIS that I want to import into GRASS so I can use some of the Grass modules.

However, I cannot find the v.in.ogr module in the processing toolbox selection.

All I have is v.in.dxf and v.in.wfs.

Is there a module/plugin/extension that contains the v.in.ogr module that I need to download?


texts:
Here, the Python for create a list of feature types but....

import sys
from pyfme import *
try:
    log=FMELogfile("shapeToTAB.log")

    log.log("Opening reader",FME_INFORM)
    reader=FMEReader("GEODATABASE_FILE")
    reader.open("C:\depot\export\TE-PMB.gdb")

    GeomList = []

    schemaFeature=FMEFeature()
    log.log("Copying schema features")
    while reader.readSchema(schemaFeature):
        GeomList.append(schemaFeature.getFeatureType())
        log.logFeature(schemaFeature)

    reader.close()
    print GeomList

except FMEException, err:
    print "FMEException: %s" % err
    sys.exit(1)

I would like return a list of feature type of my geodatabase.
Is it possible ?


texts:
I'm finishing my last degree on informatics but I'm having trouble getting Tiles/Layers.
I have a DB on MongoDB - with this Schema (GeoJSON)

db.stations.insert({
    "name": "station 1",
     "geometry": {
"coordinates": [
            50,60
       ],
         "type": "Point"
     }
 })

I'm trying to integrate with Geoserver to get WMS. I got the Mongo Plugin and I can use it to add a New Store.
After the Store. It shows empty. What should I do next? Or what am I doing wrong?


texts: TIFF, Grid, FGDB Grid, etc. - what performs better in ArcScene animation?

texts:
I am creating a flyover animation in ArcScene and I would like my surface raster to be rendered in full resolution (an exported function raster dataset - a blend of shaded relief and high resolution imagery)

I've tried to export several clips with the function raster dataset itself but the uncompressed AVI output was quite jerky. I then converted the function raster dataset into TIFF and the output video did not improve.

Will Converting to FGDB Grid or any other grid make a difference in the final animation rendering.  It plays smoothly in ArcScene itself; it's only the exported uncompressed video that seems to be quite jerky.  Any other recommendations?

What I meant by the above statement "I would like my surface raster to be rendered in full resolution" was not a question but a statement;  this raster is set and must remain set to preserve detail to HIGH (Quality Enhancement for Raster Images: HIGH)

After reading the first question (thanks @GIS-Jonathan) I realized that I might want to include my machine specs:

Intel Core i7-980X Extreme Edition 3.33Ghz CPU
OCZ DDR3 PC3-6GB RAM (12GB total)
Asus P6X58D Premium Motherboard
WD VelociRaptor 450 GB, SATA 6 Gb/s, 32 MB Cache, 10,000 RPM Hard Drive
NVIDIA Quadro 4000 w/2.0GB, Dual Link DVI Video Card
Windows 7 Professional 64bit


texts: How to georeference using Latitude and longitude in QGIS?

texts:
I'm trying to georeference a picture using Latitude and longitude.
The points are these:

top-left:  (116.30950927734372,39.970805680527725)
top-right:  (116.35620117187497,39.970805680527725)
bottom-left (116.35620117187497,39.93606595478707)
bottom-right (116.30950927734372,39.93606595478707)

When I tried to calculate the area of a polygon,the result is very small(something like e-7).I realized to my horror that the QGIS had treated the Unit as meter.What should I do?

(QGIS1.8)


texts: How does ee.Reducer.median() differ from just ee.ImageCollection.median() - Issues with a chart?

texts:
I am trying to create a chart that displays an image series from an Image Collection and there are two things that I do not understand:

In the first part of the code that you will see below, I do not get what the difference is between .median() and using a reducer: ee.Reducer.median() while plotting the chart? Does using it twice have a difference?
My chart is plotting more than 12 images when there are only 12 images in the imageCollection.

Link to code: https://code.earthengine.google.com/6f163a02975f01d6f8728feaadf98272

     var clipper = function(image){
     return image.clip(WG);
    };

    // Load the MODIS data that is being processed

    var modis = ee.ImageCollection('MODIS/006/MOD09GQ').map(clipper);
    var months = ee.List.sequence(1,12);

    // Group by month, and then reduce within groups by mean();
    // the result is an ImageCollection with one image for each
    // month.

    //ee.ImageCollection.fromImages - essentially includes a list of images

    var byMonth = ee.ImageCollection.fromImages(
      months.map(function(m) {
        return modis.filter(ee.Filter.calendarRange(m, m, 'month'))
                    .filterBounds(WG)
                    .median() // Unable to run this chunk of code without this function
                    .set('month', m);
    }));

    print("Monthly Images Median",byMonth);

    print(ui.Chart.image.seriesByRegion({
    imageCollection: byMonth,
    regions: WG,
    reducer: ee.Reducer.median(),
    band: 'sur_refl_b01',
    scale: 250,
    xProperty: 'system:index',
    seriesProperty: 'label'
    }));

The chart looks like this:


texts: FME remove polygon holes of under size ### sqft and chop the ends off

texts:
Quetion #1
I have a wokspace that takes the latest version of my data and creates buffers and then merges those to one polygon.
When 2 buffer polygons are near each other and close to create a hole. I would like to find a reasonable area size and remove only the holes. If there are small polygons for one reason or another I need to keep them.

Another question not exactly related mentions the donuthole extractor.  Do I need to extract these and then delete them?
EDIT:
I used the solution with a tester and WAS successfull.
The part I was missing was I had to put all of this "behind" my first disolver.
It wasn't an actual donut until the two parts were together. I suppose I could get fancy and put some other locgic in front to make it that way.
The other problem I was having was it is using sq meters not hectares or acres. :)
On to the second problem which I have not had any luck with yet.

Question #2
I like the rounded version not the butt of bufferer. However would like to limit the ends (after the merge) to being only a specified distance from the end of the actual line.

OK I figured out I need to  use the parameter created when I set the AreaCalculator. But I still don't see proper results.

I would like to chop off the end flat at maybe 20' offset from the actual end of the line.
Can I get a direction or pointer on accomplishing this?


texts: Portion of PDF missing in QGIS Export to PDF

texts:
I'm trying to produce a custom sized PDF (841 x 3567mm). The preview in QGIS composer looks fine but when the Export to PDF is completed I'm missing the bottom portion of the PDF image.
The composition is comprised of vector and raster layers (ECW) and some of the vector have a transparency setting. I've trying turning off the transparency to see if it is having an affect but still getting the same output, so it doesn't seem to be related.
See image attached.


texts:
I have a PDF file (10mb) which contains a NEATLINE polygon in metadata.

$ gdalinfo I5601_GeoPDF.pdf

  Driver: PDF/Geospatial PDF Files: I5601_GeoPDF.pdf Size is 5906, 3543
  Coordinate System is: PROJCS["UTM Zone 56, Southern Hemisphere",
      GEOGCS["GDA94", ... GeoTransform =
    128363.7573358897, 42.324133621968, 0.8865642156959999
    6471730.899290663, 0.8865642156959999, -42.324133621968 Metadata:   CREATION_DATE=D:20080515105609+10   CREATOR=BAE SYSTEMS MapMaker.exe
  2.9.830.19227   NEATLINE=POLYGON ((191125.203093351592543 6340529.201165121048689,188350.332458852295531 6473000.141725543886423,378308.662319099530578 6476979.201503321528435,381083.532953590154648 6344508.260942890308797,191125.203093351592543 6340529.201165121048689)) ...

Now if I extract the NEATLINE data into a .cvs and use gdalwarp to cut it:

$ for i in *.pdf; do gdalinfo $i|grep NEATLINE | awk -F= '{print "id,WKT\n1,\"" $2 "\""}' > $i.csv; gdalwarp -crop_to_cutline -cutline $i.csv -overwrite $i $i.tif; done

I get a geotiff that isn't clipped correctly:
picture of results

I've tried this on US 250k maps and it works. I'm wondering whether anybody can tell me if I'm doing something wrong, or is the metadata in the PDF wrong ?


texts:
I'm trying to geocode about 250k addresses and after.  Several crashes later I split the file into several smaller files.  I even created a test file of 4 known addresses to check if it's functioning properly.
First run of the test file was near instant.  Followed by some successful and some unsuccessful runs of files with somewhere between 3000 and 50000 addresses.
I noticed it seems to be getting much slower.  So a reran the test file of 4 addresses and it took about 5 hours.  The last file I ran was about 3000 addresses and took about 30 hours.
I've tried uninstalling(which was it's own challenge as there's no standard uninstall function btw) and reinstalling the software to get the initial speed but that didn't work.
Any suggestions are welcome.


texts: Leaflet on mouse wheel event change zoom by delta

texts:
I am trying to find a way to change zoom by delta when using mouse wheel event (scroll) in Leaflet.

E.g. current map zoom is on 0 and if "scrolled" into map, zoom should be 5.

Also, zoom change should take zoom location in account.

E.g. map is 1000px * 1000px, then if mouse wheel is done on [x:1000, y:0], then it should be zoomed in that part (in right upper corner) and not random or in center.

Attempts / Results

adding mousewheel listener on map div - zoom delta working by calling zoomIn(), zoom position random, can not get info on coordinates since it is not Leaflet event
adding mousewheel listener on tile div - not able to do anything with it
adding zoomstart listener on map - not able to do anything  with it


texts:
My question is very simple. I have a point feature class with attributes joined to the attributes of another point feature class. I would like to, if I select a feature in the target feature class, have the feature class in the joined dataset also selected. Is there any easy way to do this?

I know a table relate will select the related entry but this is not the same for a join. Also I would rather not use a relate since it is critical for me to see the attributes of both feature classes side by side.


texts:
I have a lot of coordinates in attribute table like this one:

I'm looking the way how to find and select on active layer duplicated points using Python in QGIS ? I needn't delete them but only find and show (or select).


texts: Layers and styles for all different geometric objects

texts:
I want to change the graphical style of my geoserver map.

There are like 400 different geometric objects and 20 different zoom-levels.

Do I have to create 8000 (400x20) different layers/styles and upload all of them to my geoserver?

We have a university project called OHDM (Open Historical Data Map):

OHDM on GitHub
OHDM Page

If you check out the layer previews from the geoserver of this project, you can see, that they tried to create 228 different layers by hand.

Our goal is to change the style to something like OSM looks like.


texts: Appending Feature Class From File Geodatabase to SQL Server Express Geodatabase?

texts:
I have a set of file geodatabases with the same schema as my SQL Server Express geodatabase.  Using this script, I am trying to append the feature classes of each FGDB into the corresponding feature class in the SQL Server Express GDB.

import arcpy, os

inputGDB = arcpy.GetParameterAsText(0)
sqlDB = arcpy.GetParameterAsText(1)
arcpy.env.workspace = sqlDB

sqlFCs = arcpy.ListFeatureClasses()
for fc in sqlFCs:
    sourceFC = str(fc)[14:]
    if arcpy.Exists(os.path.join(inputGDB, sourceFC)):
        arcpy.AddMessage("Appending " + sourceFC + " to " + fc)
        fcPath = sqlDB + "\\" + fc
        arcpy.Append_management(os.path.join(inputGDB,sourceFC), fcPath)

When I attempt this with my python script, I get the following error:

ERROR 000732
Target Dataset: Dataset Database Servers\WOS2UA52417ZT_SQLEXPRESS.gds\SCAT_Data (VERSION:dbo.DEFAULT)\SCAT_Data.DBO.trackjsonpnts does not exist or is not supported

It is the same if I try to simply use the Append tool from ArcToolbox.  Is there a way to append feature classes in an SQL database?


texts:
I installed qgis in debian linux and I want to run a standalone app in python. So I import qgis.core and qgis.gui in this app.
But it gives me an error:

  ImportError: libqgis_core.so.2.9.0: cannot open shared object file: No such file or directory

I installed the last version using:

deb     http://qgis.org/debian jessie main
deb-src http://qgis.org/debian jessie main

My PYTHONPATH is defined to: /usr/share/qgis/python/plugins. And I search libqgis_core and I found the libqgis._core2.8.0. So why the app is searching the 2.9.0 version?


texts: What are key things to test in ArcMap 10.1, service pack 1?

texts:
I have been asked to test service pack 1 for ArcMap 10.1 before it is provided to other users at my university. Are there any key things I should do for any SP test, or should I try to look at things that have specifically been addressed in this SP?

My plan was to create a new .mxd file and run a few basic geoprocessing commands. Realistically I'm not going to able to check everything, but are there other key things I should do?


texts: Adding new features to selection list using PyQGIS?

texts:
The code written by Joseph answering to the question Features to stay selected by default in QGIS?

layer = qgis.utils.iface.activeLayer()
# Use active layer

def select(featureAdded):
    layer.setSelectedFeatures([featureAdded])

layer.featureAdded.connect(select)
# Connect "featureAdded" event to "select" function

works, but as soon as another feature is created, a new selection is made, resulting in only the new feature been selected.

I would like to create many features that stay selected, and that are added to the features already selected. I am not good enough with scripting to create a working code.


texts:
I am currently trying to create a KMZ from a GeoTIFF by using gdal_translate. This Works fine and I am accomplishing this by using by using this command

gdal_translate.exe -of KMLSUPEROVERLAY c:\in.tif c:\out.kmz -co FORMAT=JPEG -co ALTITUDE=50.0 (Meters) and add -co ALTITUDEMODE=absolute

However, I can't seem to figure out how to add a fixed altitude relative to the ground raster than the absolute value when brought into Google Earth.

I understand that the KML file should have the elements

  <altitude>2744.0</altitude>
  <altitudeMode>relativeToGround</altitudeMode>

but cannot seem to figure out how to incorporate this into the code.


texts: Collecting and merging Service Analysis Line features using ArcGIS Pro ModelBuilder

texts:
I try to collect and merge Service Analysis Line features after iteration. I name my analysis layers with the iteration values, but it doesn't help. The output feature datasets of the analysis named with consequent numbering (ServiceArea1, ServiceArea2, etc.) and the line feature classes like (SALines1, SALines2, etc.). The line features represents the right output (I tested). When I try to collect data, all the Collected Values call "Lines". I tried to write out the lines to feature class for every run by selecting the child data and feature class to feature class tool, but it didn't work. The result was not the output of the analysis result (SALines1, SALines2) but the same line features for every analysis result.
It seems that I miss a step that refers to the saved line feature classes in the created feature datasets. With this model I'm going to generate hundreds of Service Area layers and need to collect and merge the line outputs either as table or feature class and afterwards generate different statistics. Adding the results manually is out of the question.

Can anybody help with this?

[


texts: Automate the writing records from csv to shapefile using pyshp module

texts:
I have the following data in a csv file('file.csv')

Date_Time,Latitude,Longitude,Wind_Speed,Status
09-10-2013/1130,13,93,50-60 gusting to 70,Observed
09-10-2013/2330,13.501,92.501,70-80 gusting to 90,Observed

I am able to convert the above data into a point shapefile using the pyshp module with the following code
import csv, shapefile
point = shapefile.Writer(shapefile.POINT)
with open('file.csv', 'rb') as csvfile:
     csvreader = csv.DictReader(csvfile)
     header = csvreader.fieldnames
     for field in header:
         points.field(field) ## I could define the fields on the go in a loop
         for row in csvreader:
            parts.append([(float(row['Longitude'])), float(row['Latitude'])])
            points.point((float(row['Longitude'])),float(row['Latitude']))
            points.record(row['Date_Time'], row['Latitude'], row['Longitude'], row['Wind_Speed'], row['Status'])
            # for writing records each record has to be specified
points.save('points')

I would like to automate the record part i.e., point.record(row['Longitude'], row['Latitude'], row['Date_Time'],... so if the no.of fields are more writing each field name in "point.record(...)" can be automated.


texts: Joining data of two different datasets from varying distances and integrating into one dataset in ArcGIS

texts:
I have a pointfile of some 200 electrical boxes throughout the county. While the electrical boxes are "mapped", they are done so by using township range section data, so there actual location is not known. They are simply at the centroid of a section or subsection. My task is to verify their real world location and move them from the cetroid of the section or subsection to that real location.

I have a 2nd pointfile of electrical boxes that have been GPS located by workers from a utility company. This pointfile contains some 70 features that I know their location is correct. However, only 40 of the points in this file correspond to the 200 in the above file.

Both files contain a unique identifying number for all of the electrical boxes allowing for a 1:1 join of the files to be completed as the GPS points need the attribute data that is present on the county dat. However, since a tabular join like that doesn't actually change any spatial data, it isn't enough to call the job finished. I need to get the attribute data from the 200 to the matching 40, while eliminating those 40 points from the 200. I then need merge the 40 GPS'd points into the county dataset, so all my points are unified in one single file.

I cannot do a spatial join on these two files, because the distances between the GPS verified and the 200 from the county are not uniform due to the TSR mapping of the county data. In other words, joining them based off distance X wouldn't work because the closest points will not necessarily match in real life.

What I really need is to be able to select attributes across two different datasets. I tried merging the two and then running a selection of where UNIQUE_NO = UNIQUE_NO1, but while it was verified it returned no results.

Any idea how to move forward on this? Working in ArcMap 10.2


texts: Is it possible to add a simple image or graphic to the map canvas?

texts: Is it possible to add a simple image or graphic to the map canvas in QGIS (e.g. a symbol marker)?


texts: Using fme to create shapefile from several kmz files, labeling points in shapefile with name of kmz?

texts: I am trying to create a shapefile of point data. I have lots of kmz files each with a single point, and each named specifically. I'd like to end up with one shapefile of points all with a name attribute that is named for the original kmz file. Sort of a reverse fanout.


texts: Raster layer statistics QGIS 3.4 doesn't give information about pixel count

texts:
I created a plugin in QGIS 2 AND I used the qgis:rasterlayerstatistics , I could get the 'COUNT' information : pixels number, but when I migrate my plugin to version 3.4 I found out that this information (COUNT) is not available anymore .

is there any other way to have raster statistics about pixels number ?


texts: Import polygon shapefile with SPCC Sections

texts:
I have a data set with records containing Township, Range and Section information but no X,Y coordinates.  I would like to bring in the data to ArcGIS with each section as a polygon or a x,y point in the center of each polygon.  Is this possible?

   #  name  date    time        T.N R.W S.E
    1 801   11/16/061100    A   64  9   33SE
    2 801   1/7/05  1122    A   63  9   7SW
    3 801   6/19/00 1352    A   63  8   34NW
    4 801   9/25/00 1455    A   63  8   26SW

What I am ultimately hoping to do is to make a polyline shapefile with the lines moving from one observation to the next but was thinking of starting with points and then working to lines.


texts: Calculation of intermediate points along shortest path not working correctly

texts:
I am trying to implement the calculation of points along the shortest path using the algorithm shown here:
LatLng getGreatCirclePoint(double lat1, double lon1, double lat2, double lon2) {
  double d = 2*asin(sqrt(pow((sin((lat1-lat2)/2)),2) + cos(lat1)*cos(lat2)*pow((sin((lon1-lon2)/2)),2)));
  double f = 0.5;
  double A = sin((1-f)*d)/sin(d);
  double B = sin(f*d)/sin(d);
  double x = A*cos(lat1)*cos(lon1) +  B*cos(lat2)*cos(lon2);
  double y = A*cos(lat1)*sin(lon1) +  B*cos(lat2)*sin(lon2);
  double z = A*sin(lat1) + B*sin(lat2);
  double lat=atan2(z,sqrt(pow(x,2)+pow(y,2)));
  double lon=atan2(y,x);
  LatLng point = LatLng(lat,lon);
  return point;
}

Unfortunately this does not return the correct point but rather this here:

Where did I go wrong in my implementation?


texts: Access string output of QGIS 3 processing script in graphical modeler

texts:
In the good ould days of QGIS 2 I created a processing script, that takes as input a table and a table field (table has only one record) and returns the string attribute value of the table field:

##mygroup=group
##myname=name
##input=table
##a=field input
##text=output string

formeln=processing.getObject(input)
for fo in formeln.getFeatures():
    text = fo[a]

Embedding this script in graphical modeler in QGIS 2 and using the output e.g. in a Field calculator algorithm lets me chose its output string as input for (here) the formula:

So far so good, now on to QGIS 3. At least I succeeded (?) in rewriting this script, my result:

from PyQt5.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                   QgsFeatureSink,
                   QgsProcessingException,
                   QgsProcessingAlgorithm,
                   QgsProcessingParameterFeatureSource,
                   QgsProcessingParameterField,
                   QgsProcessingParameterString,
                   QgsProcessingParameterFeatureSink,
                   )
import processing

class ExampleProcessingAlgorithm(QgsProcessingAlgorithm):
    FORMULA_COLLECTION = 'FORMULA_COLLECTION'
    FORMULA_ID = 'FORMULA_ID'
    OUTPUT_STRING = 'OUTPUT'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return ExampleProcessingAlgorithm()

    def name(self):
        return 'output_table_field_value_as_string'

    def displayName(self):
        return self.tr('Table field value as string')

    def group(self):
        return self.tr('giswg')

    def groupId(self):
        return 'giswg'

    def shortHelpString(self):
        return self.tr("This returns the value of the specified field of a table. Implies, that the table has exacly one record.")

    def initAlgorithm(self, config=None):
        #input parameter
        self.addParameter(QgsProcessingParameterFeatureSource(self.FORMULA_COLLECTION, self.tr('The formula collection'), [QgsProcessing.TypeVector]))
        self.addParameter(QgsProcessingParameterField(self.FORMULA_ID, None, self.tr('The formula identificator'), self.FORMULA_COLLECTION, QgsProcessingParameterField.String))
        #output parameter
        self.addParameter(QgsProcessingParameterString(self.OUTPUT_STRING, self.tr('Output string'), ''))

    def processAlgorithm(self, parameters, context, feedback):
        formula_collection = self.parameterAsSource(parameters, self.FORMULA_COLLECTION, context)
        formula_id = self.parameterAsString(parameters, self.FORMULA_ID, context)

        out_string = self.parameterAsString(parameters, self.OUTPUT_STRING, context)

        #this portion is taken from the old QGIS 2 script
        for fo in formula_collection.getFeatures():
            out_string = fo[formula_id]

        return {self.OUTPUT_STRING: out_string}

Now I would like to do the same in QGIS 3, but the 'Formula' field refuses to let me chose the output of my script as input:

Any suggestions what might be wrong with my script or what else I'm probably missing?


texts: Splitting polyline with ArcGIS Pro snapping environment?

texts:
I would like to take advantage of the 64-bit environment and speedier screen refreshing for some editing we do.

Although I personally don't find the ribbon interface very intuitive, the ArcGIS Pro editing environment seems quite nice. Except I found a problem....

When I set my snapping to only Intersection and Vertex I am unable to SPLIT my polylines at intersections and vertices. Instead, in most situations, the line is split somewhere near the intersection or vertex but not exactly at intersection/vertex.

The problem occurs when snapping options are set to either pixels or units. I tested it on 3 different machines, using different spatial references, and in both ArcGIS Pro 1.0.0 and 1.0.2.

I expect the line to be spit exactly where it snaps at vertices or intersections but I am simply not able to.

a result of line split when snapped to vertices:

a result of a line split at intersection:


texts: Conversion of an OSM Multipolygon to a Shapefile Polygon Using QGIS

texts:
I am quite new to the topic of handling geodata, so I face my first big challenge now:

I downloaded an osm file containing information about any kind of data related to power generators, transformer stations etc. (http://enipedia.tudelft.nl/wiki/Extracting_Power_Data_from_OpenStreetMap). I successfully loaded this file to QGIS, where I started using the save as option to simply save it as a shapefile. This worked fine for points and lines, but not for multipolygons and multilines. As I think I know now it´s because shapefiles do not support these types (correct?). So I tried an advice from our forum to use the geometry tool multi to onepart. But now when I try this I get the error

  len() should return >= 0 See log for more details

Has anybody faced the same problem? My overall goal here is just to use this dataset as a shapefile.

Update: Maybe also an relevant information is that the properties of the layer display "Number of objects in this layer: -1"


texts: Xamarin Forms equivalent of InfoTemplate

texts:
In the past, I've used InfoTemplates to display layer attributes using the ArcGIS API for JavaScript.

Now, I'm trying to make a Xamarin Forms app using the ArcGIS SDK for .NET. I'm stuck trying to figure out a way to get similar functionality. I want to be able to display a callout for each layer when a user touches the screen. I found this, but the map I'm consuming doesn't have any popups.

What I really want to do is to get the layer's attributes and use them to display a callout (or multiple callouts - one for each layer).

Is this possible with this sdk?

Is there a better way to accomplish what I'm trying to do?

I'm not looking for someone to write the code for me, just a nudge in the right direction or a link to something I can follow.

** Edit **

Ok, I'm able to get the attributes for a layer by doing the following:

 private async void MyMapView_GeoViewTapped(object sender, GeoViewInputEventArgs e)
    {
        try
        {
            var queryParams = new QueryParameters();

            queryParams.WhereClause = "1 = 1";
            var results = await _featureTable.QueryFeaturesAsync(queryParams);

            foreach(var result in results)
            {
                var attst = result.Attributes;
            }
        }
        catch (Exception ex)
        {
            DisplayAlert("Ugh!", ex.Message, "Close");
        }
    }

Now, how can I get the layer(s) that contain a given point where the user touched (e.Location.X and e.Location.Y)?


texts: Limiting returned records for an interaction in OpenLayers 3

texts:
I have the following script that returns a list of all the customer names at that point when it is clicked on. I am looking to limit this to only show the first (or just one) record.Could this be achieved through a modification of the following?

   function customer(browserEvent) {
    var coordinate = browserEvent.coordinate;
    var pixel = map.getPixelFromCoordinate(coordinate);
    var el = document.getElementById('customer, single');
    el.innerHTML = '';
    map.forEachFeatureAtPixel(pixel, function(feature) {
      el.innerHTML += feature.get('Proj_Client') + '<br>';
    });
  }
  map.on('pointerdown', customer);


texts: How could I turn my data into a heatmap / intensity map?

texts:
I have database that tracks a large number of entities. Each entity has a 'coverage area' or 'catchment area' - a radius in miles from a specified UK postcode, eg

Leyton office     4.3m    E11 4LL
Stratford office  5.2m    E15 6ZZ
East Ham office   2.1m    E4 8QQ

I want to turn this into a geographic intensity map - a map of a region which is shaded or coloured to indicate the number of entities that cover each point on the map. Segmentation into small regions is also acceptable.

This is probably a one-off project, so it's unlikely my organization will provide money to buy software. That means zero or at least low-cost solutions are preferable.


texts: How to get geoJSON from spring REST Controller?

texts:
I am developing a GIS application with java(Spring-4.1.5 + Hibernate-4.3.8) and OpenLayers. For this project I use GeoTools-13RC, HibernateSptial-4.3, jts-1.13 and jackson-2.5.
In this project, I have a layer in client side and in server, I save the features of this layer in a class. I defined the class below:

@Entity
@Table(name="MyPoint")
public class MyPoint{
    @id
    @GeneratedValue(strategy = GenerationType.IDENTITY)
    private Long Id;

    @Column
    private String name;

    @Type(type = "org.hibernate.spatial.GeometryType")
    @Column(name = "the_geom")
    private Point geometry;

    /*
    *
    Getter and Setter
    *
    */
}

In start up of application, I need to init the layer in client side. for this, I need return from server side a json string to client for this layer. I don't want to use ST_AsGeoJson or other matches. I use Spring REST controller for returning my Entity.

What do I do?


texts: Installing scipy.spatial with ArcPy gives DLL load failed error on import?

texts:
I am running python 2.7.6 from an ArcGIS installation. I have tried to install scipy from the Python Exchange Packages here, using the scipy-0.14.0.win32-py2.7.exe executable. I get no error importing scipy, but when I try to import scipy.spatial I get the following error:

In [1]: import scipy

In [2]: import scipy.spatial
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-2-6f8fbf892cc9> in <module>()
----> 1 import scipy.spatial

C:\Python27\ArcGIS10.2\lib\site-packages\scipy\spatial\__init__.py in <module>()
     89 from .kdtree import *
     90 from .ckdtree import *
---> 91 from .qhull import *
     92 from ._plotutils import *
     93

ImportError: DLL load failed: %1 is not a valid Win32 application.

Anyone had success install scipy with spatial functions on top of Esri's Python site-package (ArcPy)?


texts: Get Feature id in OpenLayers from GeoServer/PostGIS

texts:
With OpenLayers, I'm creating layers like ol.layer.VectorTile

var projection_epsg_name = 'EPSG:900913';
var geoserverLayerName = 'teste:map_shape2';

const tileLoadURL= 'http://127.0.0.1:8081/geoserver/gwc/service/tms/1.0.0/' + geoserverLayerName + '@' + projection_epsg_name + '@pbf/{z}/{x}/{-y}.pbf';

// const tileLoadURL = "http://127.0.0.1:8081/geoserver/gwc/demo/teste:map_shape2?gridSet=EPSG:900913&format=application/vnd.mapbox-vector-tile"
// const tileLoadURL = "http://localhost:8081/geoserver/gwc/service/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&LAYER=" + geoserverLayerName + "&STYLE=&TILEMATRIX=" + projection_epsg_name + ":{z}&TILEMATRIXSET=" + projection_epsg_name + "&FORMAT=application/vnd.mapbox-vector-tile&TILECOL={x}&TILEROW={y}";

var vectorSource = new ol.source.VectorTile({
    tilePixelRatio: 1,
    tileGrid: ol.tilegrid.createXYZ({maxZoom: 19}),
    format: new ol.format.MVT(),
    url: tileLoadURL,
});

var vectorTile = new ol.layer.VectorTile({
    style:styleFunction,
    source:vectorSource,
})

map.addLayer(vectorTile);

I need id to perform a query in the database, but I do not have this attribute. My data comes from PostGIS, and GeoServer is configured in such a way that it displays the id. How can I have this attribute using OpenLayers?


texts: Syntax for converting tuple into string with Merge in ArcPy

texts:
I have a tuple of unknown elements, which I want to convert into a string for using in the arpy command
arcpy.Merge_management([Input1, Input2, ..., Inputx], "D:\\test.shp")

In my example I have only three elements in the tuple and tried to convert this element into one string as the inputs for the command. I do this with following command:
Input = "['" + "', '". join((str(n) for n in args)) + "']"

If I display the variable [Input], I get a plausible result.
['D:\\Buildings\\Gemeinden\8537\8537_Change_Detection_TLM_LAZ_pre.shp', 'D:\\Buildings\\Gemeinden\8537\8537_Change_Detection_TLM_AVS_pre.shp', 'D:\\Buildings\\Gemeinden\8537\8537_Change_Detection_TLM_DEM_pre.shp']

Unfortunately I got always following error message, when I tried to integrate the variable into the Merge command:
arcpy.Merge_management(Input, "D:\\test.shp")

ERROR 000732: Input Datasets: Dataset '['D:\\Buildings\\Gemeinden\8537\8537_Change_Detection_TLM_LAZ_pre.shp', 'D:\\Buildings\\Gemeinden\8537\8537_Change_Detection_TLM_AVS_pre.shp', 'D:\\Buildings\\Gemeinden\8537\8537_Change_Detection_TLM_DEM_pre.shp']' does not exist or is not supported

Can anyone help to solve the hopefully small problem?


texts: I am trying to make mosaic of different landsat images. First I merge the required bands of each raster than try to merge the compiled bands but in the end getting nodata value at the overlapping area. Can anyone help to resolve it?


texts: Automate publishing of WMS or WFS Layer in GeoServer using OpenLayers 3?

texts:
How to automate publishing of a WMS or WFS Layer in Geoserver given any data source (vector/shape-files/GeoJson...) using OpenLayers 3 script

There is a possibility to automate Geoserver operations using REST API, but I couldn't found how to link it with OpenLayers 3.


texts: Convert coordinates to geometry (polygon) using QGIS

texts:
I have an Excel sheet with grid data including geographical coordinates. I know that the coordinates represent areas, i.e., a multi(polygon). However, the data lacks geometry information. I'm looking for an open source way of converting these coordinates into the corresponding geometries, mainly polygons and points. I'm hoping that I can do this using QGIS. I've looked around for information on how to achieve this, but no luck so far.

I'm not even sure this is possible. I'm still researching, but any idea is welcome.

UPDATE: I know that for a polygon, there are supposed to be at least 3 coordinates. But these are just given as (x,y). Nothing that conforms to a linear ring, for instance.I'm thinking that the coordinates in my Excel sheet represent the center of the polygon. Is there away of telling what the polygon is. This sounds far fetched, perhaps because I've never done anything like this before. Your opinions will be highly appreciated.

UPDATE: Sample Data

I plan to separate the 'Geolocation' column into x and y columns. Having thought about it for sometime, I think it's impossible. But any idea that might help will be appreciated.


texts: Linear Referencing using ArcGIS Desktop?

texts:
I'm a fairly experienced ArcGIS 10 user, but am new to Linear Referencing.

I have a road segment dataset, stored as a simple line feature class in File Geodatabase, and I have a table of "events" in the same geodatabase.
The road lines have a Segment_ID attribute, and also a Section_Length.
The events table have the same Segment_ID along with other information, including a "Length" attribute which we want to use to dynamically segment each road segment at that chainage.

My understanding (which may be wrong!) is that we have to use the "Create Routes" tool followed by "Create Route Event Layer" which will segment the lines at the chainages defined in the "event" table.

The problem seems to be that the direction of the lines (which are highways with specific directions of travel) in the original feature class are ignored when creating the Routes, apparently necessitating a manual edit to ensure they are correct.

How can this be avoided to perform the linear referencing?


texts: Superimpose OTB application via Sextante cannot load parameter

texts:
I am trying to perform a pansharpening using OTB via sextante. According to the cookbook I first need to superimpose the Pan and XPS images:

http://www.orfeo-toolbox.org/CookBook/CookBooksu36.html

However I get the following error in my log:

ERROR: Troubles loading parameter, please check your line argument...
This is the Superimpose application, version 4.2.1
Using available image metadata, project one image onto another one

I am using Qgis 2.4.0 and OTB both installed via Osgeo4W 32 bit installer.
My images are WorldView2.

Does anyone know what is causing this error