texts: Create layer using selected features from WFS layer using PyQGIS

texts:
I am making a plugin that loads a WFS layer and lets the user filter out the information they need and create a layer based on that filtering. The new layer is created from the original layer based on expressions that are created from the users filtering. This prosess can be very slow as the WFS layer is rather large. I would therfore like to know if there is any faster or betther ways of doing this.

Here is an example of my code:

expr_string = "\"kommune\"=301 AND \"funksjon\"='kultur'" #example of an expr_string
expr = QgsExpression(expr_string)
it = vlayer.getFeatures( QgsFeatureRequest( expr ) )
ids = [i.id() for i in it]
vlayer.setSelectedFeatures( ids )
selectedFeatures = vlayer.selectedFeatures()
tempLayer = QgsVectorLayer("Point?crs=epsg:4326", "newLayer", "memory")
QgsMapLayerRegistry.instance().addMapLayer(tempLayer)
temp_data = tempLayer.dataProvider()
attr = vlayer.dataProvider().fields().toList()
temp_data.addAttributes(attr)
tempLayer.updateFields()
temp_data.addFeatures(selectedFeatures)

This code is based on an anwer from this post: Quickly make server side spatial filters for WFS layers based on other layer's selected geometries in QGIS


texts: Variable Buffer on Street Centerline using ArcGIS

texts:
I have a street centerline file where all road segments have a functional class of Arterial (A), Collector (C), Local (L). I also have a point file of each intersection.

I would like to do a crash analysis using 150' buffers at each intersection. The trick is that the buffer I want is not just a uniform 150' buffer off the intersection point.  I want to buffer the centerlines for a specified length from the intersection point based on the functional class: A = 400', C = 250', L = 150'.

For example where an Arterial (A) intersects a Collector (C), I want a 150' buffer of a 400' line segment from the intersection point for the Arterial roadway and a 150' buffer of a 150' segment for the Local Roadway.  The image shows what I am trying to get at.

I'm thinking I need to first find a way to create line features of lines with of the specified length coded with the functional class. That's where I'm stuck. Then it would be just a simple matter of buffering that result.

Anyone have any thoughts? I have ArcGIS 10.6.1, Network Analyst and Spatial Analyst.

Here are examples of the results using FelixIP steps below.

In areas where the buffers do not overlap it works great:

For intersections that are complex or composed of dual lines it does not work too well. Currently looking for a way to simplify these areas.

Here is a good example of the final buffer I'm trying to achieve:

Results I have using Emil's script.


texts:
I am using QGIS and C++.

When I create a layer of polygons with labels in the center I am able to use a RuleBasedRenderer to control the background colour of the polygon. How do I similarly control the colour of the label? The background colours are limited so the combinations will not be too long. So when for instance, the background is black, I want the label to be white.

Is it possible to have white labels on black polygons and black text on white polygons exist on the same layer?


texts: Can't convert a WKT with a Z value into a SQL Geometry value

texts:
I'm using nettopologysuite in c# to build a nettopology geometry object, and then save it as a geometry object in SQL.  Generally this works fine, except i recent came up with (from a shapefile) the WKT:

  POINT Z(3.59140026404647 50.7568317130073 0)

which broke everything.  Both if i do it in code:

var tmp = SqlGeometry.Parse(new SqlString(geometryObject.AsText()));

or if i do it in the SQL procedure itself:

declare @tmpGeom as Geometry = null;
    if (@GeometryObj is not null)
        Begin
            set @tmpGeom=geometry::STGeomFromText(@GeometryObj ,4326)
            if (@tmpGeom.STIsValid() = 0)
            Begin
                set @tmpGeom = @tmpGeom.MakeValid()
            end
        end

neither work.  I get an error similar to:

  24142: Expected "(" at position 6. The input has "Z".
  System.FormatException Source = Microsoft.SqlServer.Types at
  Microsoft.SqlServer.Types.WellKnownTextReader.RecognizeToken(Char
  token) at
  Microsoft.SqlServer.Types.WellKnownTextReader.ParsePointText(Boolean
  parseParentheses) at
  Microsoft.SqlServer.Types.WellKnownTextReader.ParseTaggedText(OpenGisType
  type) at
  Microsoft.SqlServer.Types.WellKnownTextReader.Read(OpenGisType type,
  Int32 srid) at
  Microsoft.SqlServer.Types.SqlGeometry.GeometryFromText(OpenGisType
  type, SqlChars text, Int32 srid) at
  Microsoft.SqlServer.Types.SqlGeometry.Parse(SqlString s)

Any ideas how to convert the POINT WKT (which does read as 'IsValid' in nettopology) into something SQL can use?


texts: Lines overlaping in same layer - PostGIS

texts: I have lines overlapping in one layer as in image below. How do I delete one of these lines and leave the other one in the dataset?


texts: Union, Dissolve, Buffer -- multiple layers, all in one step?

texts:
I'm a new user of Quantum GIS (1.7.3) and I need to repeat a task multiple times, so I'm wondering if there is a quicker way to do the following (with three or more polygon layers - each layer has a single polygon):

Union
Dissolve
Buffer

The effect will be to have three separate single-polygons shapefiles merged into one (larger) single polygon shapefile, with a buffer around it.

Answers in PostGIS would be useful too as my data is stored there.  Thanks!

Edit: Adding PostGIS information

The multiple polygons I want to Union-Dissolve-Buffer are the result of a catchment basin point table that I've created in PostGIS, converted (in R) to an alpha-hull (i.e. concave hull) polygon shapefile using the CRAN alphahull package.  To create the original catchment table in PostGIS, for example, I do this:

CREATE TABLE catchment_10304 as
select
    id,
    geom_4269,
    (select sum(cost) from (
       SELECT * FROM shortest_path('
       SELECT gid AS id,
          start_id::int4 AS source,
          end_id::int4 AS target,
          newcost::float8 AS cost
       FROM network',
       10304,
       id,
       false,
       false)) as foo ) as cost
from node;

CREATE TABLE catchment1_10304 AS
SELECT * FROM catchment_10304 WHERE cost < 1364;

Then I do it for three different nodes (representing the boundaries of a region).  To create the total catchment area for the region, I want to Union-dissolve-buffer with the three or more individual catchment area tables.


texts:
I want to clip a line layer that I have, only within a polygon. Used Vector -> Geoprocessing Tools -> Clip, and I do get most of line I want that is EXACTLY within the polygon layer. Thing is the line sometimes goes a little bit outside the polygon, but I want to include those parts of the line as well, to get a continuous line. In the picture, the green line is the clipped line, while the pink line is an example of the line that isn't clipped but I want it to be.

Any way to set the clip function to get line with certain radius off the polygon? Or clip manually those parts and then add them to the clipped layer?


texts: ArcMap Symbol levels get mixed up in multiscale representations

texts:
I have created several representations for a road network on different scales. As many of the roads are displayed as cased lines, I use symbol levels to control the drawing order. I then have saved everything in layer files to import or to process with further tools (overpass etc.)

I created a template document with multiple layers for the different scales. They show/don't show depending on the scale set, ie

Scale 1 visible at 1:3,000,000 --> Representation 1
Scale 2 visible at 1:1,000,000 --> Representation 2
...

A representation has different rules and I sort them in the symbol level window corresponding to the hierarchy of the road and attributes like bridges and tunnels, for instance:

Representation 1:

Motorways (Bridges)
Primary Roads (Bridges)
...
Tertiary Roads
...
Pedestrian (Tunnel)

and so on.

These rules and their order are properly saved in the layer files and show just fine when I load a single file into a document or the template. However loading more of my multiscale layers into it, they start to mix up. Have a look here on the first load of the scale of 1:500,000:

Now if I add another layer file, for instance 1:50,000, these entries totally mix up and consecutively mess up the drawing and the export:

That is certainly complete rubbish as I do not want to re-sort them everytime (they do get updated once a week) and furthermore will not save properly after a third of forth layer is added.

I tried a clean document and all other tricks that I knew, but I cannot resolve this issue. Some scales will maintain their original order, some won't. But I cannot try every combination of 7 scales and different representations for roads, railways, etc.

I would be grateful for any hint or help how to workaround (as I am pretty sure it is caused by buggy ArcMap itself).


texts: Creating Raster Colorizer with ArcPy for ArcGIS Pro

texts:
I'm having trouble trying to write the Colorizer code for an NDVI Raster in ArcGIS Pro.
The raster is single band, and if I set the discrete colorizer through the Symbology window, I then can access the colorizer through ArcPy after creating a lyr.symbology object, and manipulate it from there.
However I'm at a loss as to how to create the colorizer from scratch.
What workflow will work?
I presume it it has to do with the _getcolorizer option in symbology, but that's as far as I've got.
    n_sim = lyr.symbology
    n_sim._getColorizer('RasterClassifyColorizer')
    TypeError                                 Traceback (most recent call last)
    In  [8]:
    Line 1:     n_sim._getColorizer('RasterClassifyColorizer')

    TypeError: _getColorizer() takes 1 positional argument but 2 were given

This returns the above error
    n_sim._getColorizer = 'RasterClassifyColorizer'

This doesn't return an error, but the workflow escapes me.


texts:
I added the new vector layer in QGIS and I want to draw a circle of 500M or 1000M from center.

I also added the pluginn Rectangles ovals digitizing in QGIS. But the options of this pluginn is not active.

How do I use this plugin?


texts: Show all layers stocked on a GeoServer workspace in a LayerTree (Openlayers3/ExtJS6/GeoExt3)

texts:
I'm building a layer tree and I would like it to be updated automatically when the page is loaded with layer that are in a specific workspace on GeoServer.

In my layer tree I have a group containing background layers such as Bing satellite and OSM, and another group that will contain the GeoServer layers.

How can I request the layers from the workspace and add them automatically?


texts: How to automatically input a raster value into Raster Calculator?

texts:
At Arc 10.5.1 I have a raster containing a single cell, with a value of 9300.  From it I need to create an output raster whose cell values decay away from the input raster's cell location, starting with the initial value of 9300.

Using the raster calculator, this equation works correctly:

"out_raster" = 9300/EucDistance("in_raster")

Notice that I have manually typed in the in_raster value, 9300.

Instead of manually entering the value, I now want to use the input raster itself, with its value of 9300.  This would allow me to place the equation in ModelBuilder and automate the process with additional single-cell input rasters, whose value may vary from layer to layer (i.e. they're not all 9300).  So I replace the text value 9300 with the input raster:

"out_raster" = "in_raster"/EucDistance("in_raster")

I expected the output to be identical to manually entering the input raster's value, but this output raster is useless; it does not display any cells, and its Layer Property values range from -3.40282e-038 to 3.40282e+038, which I think are the storage limits of the raster.

So how can I "automate" this raster equation without having to manually type in each input raster's value?


texts: Converting KML to layer in ArcGIS Desktop?

texts:
I would like to convert KML to layer, but I did not find the tool (from KML to layer) in the conversion tools in my ArcMap 9.3 version.

What is the issue and how can I solve it?


texts: Seeking guidance for ArcGIS customization with VBA

texts:
I am looking for guidance, suggestions, tips, or insights in performing a sequence of operations using VB macros for ArcGIS.

A) Creating a new table

How to write a macro in VB that creates a new table with its associated fields, loads the table to the current data frame, and adds the required records?  For instance, a parcel table with fields

P_ID (Long)
L_Use_Code (Long)
L_Use_Name (string)
Owner_ID (Long)

B) Creating a new feature class

How to write a VB macro that creates a new shapefile feature class called “Parcel” and adds it to the map window?

C) Relating and joining data to layers

How to write a macro in VB that relates a table (as in (A)) to a layer (as in (B))?

D) Identification

How to write a macro in VB that shows a feature cursor operation to display the values of a given field for each row in a feature layer table?


texts:
I have this attribute table, I would like plot on QGIS Atlas, a line graphic with last 7 days the fire events, from the day of the Atlas Page Name.

For example;

Atlas Page name is '2018-08-25', the graphic draw line of last 7 days automatically like the image below the table

Line graphic exemple


texts: Ogr2ogr does not support file name with non-ASCII character?

texts:
I tried to use the following script to convert shapefile to mif:

ogr2ogr -f "MapInfo File" D:\docs\arcgis_data\shp_mif\BOU.mif D:\docs\arcgis_data\shp\BOU.shp

it works.

But I found that if the name of the shapefile has non-ASCII characters like this:

D:\docs\arcgis_data\shp\臺北.shp

Then the covert will fail with message:

Unable to open datasource `D:\docs\arcgis_data\shp\臺北.shp' with the following drivers
......

Any idea to fix it?


texts: Exporting ASCII file from 'stars' package

texts:
I'm playing around with the stars package (0.4-1) in R, but I can't get how to export a grid object in ASCII format, as the traditional writeRaster(r, filename=file.path(dir.list[j], 'test.asc'), format="ascii") would do.

On one side I tried

write_stars(r, dsn=file.path(dir.list[j], 'test.asc'))

but I get

Error in CPL_write_gdal(mat, file, driver, options, type, dims, from,  :
  driver not recognized.

On the other hand, if I try to convert it to a RasterLayer using the provided function I get different errors:

r = as(r)
Error in .identC(thisClass, Class) :
  "Class"  is not specified and doesn't have a default value

r = stars::as(r)
Error: 'as' is not an object from 'namespace:stars'

The help page is not that helpful and I don't find any other info in the web.


texts: Three different distance values, which one should I use?

texts:
I'm implementing a script that (ideally) has a threshold of <1 meter. However, I have three different distances being pulled from these points that differ by multiple meters. I was wondering which I should go by.
From my understanding, geopy uses a model of some sort for a more accurate number whereas mpu assumes a sphere and uses a "haversine" formula. The other distance derived referred to as "delta" in the script is what was there previously and I'm not sure how to interpret that one in terms of comparing it to the other two. Is there one you would recommend I use? It is reading a random set of survey points, usually around the coast of US.
deg_to_m = 111139
lon_scale = np.cos(np.deg2rad(new_coords[1]))
delta_x = lon_scale * deg_to_m * (old_coords[0] - new_coords[0])
delta_y = deg_to_m * (old_coords[1] - new_coords[1])
delta_dist = np.hypot(delta_x, delta_y)

print(f"delta distance from point {new_coords} and {old_coords}= {delta_dist}")
print(f"geopy distance from point {new_coords} and {old_coords} = {geopy.distance.geodesic((new_coords[1],new_coords[0]), (old_coords[1],old_coords[0]) ).m}")
print(f"mpu   distance from point {new_coords} and {old_coords} = {(mpu.haversine_distance((new_coords[1],new_coords[0]), (old_coords[1],old_coords[0])) )*1000}")

Sample Results
delta distance from point (-73.54103769993867, 40.967273802785186, 17.68) and (-73.55888019993945, 40.96432750278418, 18.35)= 1532.7166288425676
geopy distance from point (-73.54103769993867, 40.967273802785186, 17.68) and (-73.55888019993945, 40.96432750278418, 18.35) = 1537.1838851634725
mpu   distance from point (-73.54103769993867, 40.967273802785186, 17.68) and (-73.55888019993945, 40.96432750278418, 18.35) = 1533.5205811777234

delta distance from point (-73.54103769993867, 40.967273802785186, 17.68) and (-73.55625289993931, 40.96710510278517, 18.32)= 1276.9866000176955
geopy distance from point (-73.54103769993867, 40.967273802785186, 17.68) and (-73.55625289993931, 40.96710510278517, 18.32) = 1280.9052835621858
mpu   distance from point (-73.54103769993867, 40.967273802785186, 17.68) and (-73.55625289993931, 40.96710510278517, 18.32) = 1277.6308281008874


texts:
Possible Duplicate:
  DBF creation and manipulation without excel 2003

Hi -

I'm trying to convert an excel file to a shapefile in arview 3.3. using suggestions by Bill Huber (on the ESRI forums site).  I've completed the recommened 'recipe' of requirements on this link Bill  referred:

http://www.quantdec.com/SYSEN597/GTKAV/section4/chapter_15a.htm#Recipe

but I am now stuck trying to convert .xls to .dbf because there is not an option in MS Excel 2007 to convert  to .dbf (dbase III or IV) that I can see - its all variations of .xls then .csv, .txt and other types.

Perhaps .csv is the file type I  need, but the 'recipe' says dbase III or IV, neither of which are listed  as options in excel.  Does anyone know  where to go from here?  If I can just get it to dbf I should be good (assuming my file formatting worked).


texts: Create time series and export it to csv GEE

texts:
I'm new to GEE and I'm trying to create time series of MODIS LST for a selected point. I first chose a point on the map and my script is:

var POI = geometry;
// Get  temperature data for 1 year.
var modis = ee.ImageCollection('MODIS/006/MOD11A1');
var modisLST = modis.filterBounds(POI)
    .filterDate('2003-12-25', '2004-02-25')
    .select('LST_Day_1km');

// Convert temperature to Celsius.
modisLST = modisLST.map(function(image) {
  return image.multiply(0.02).subtract(273.15);
});

var tempTimeSeries = ui.Chart.image.seriesByRegion({
  imageCollection: modisLST,
  regions: POI,
  reducer: ee.Reducer.mean(),
  xProperty: 'system:time_start',
  seriesProperty: 'label'
});

When I run this script, nothing happens and after 2 minutes I get:

Too many concurrent aggregations

How do I export the time series from here to csv or xlsx files? and also (but not super important) how do I display it on the GEE platform.


texts: Geoserver not recognizing jsonb attribute in SQL view of PostGIS db

texts:
I have an instance of Geoserver 2.8.2 running and a PostGIS database. I am attempting to create a FeatureType based on an SQL view of the following table (viewed through psql \d):

                 Table "public.rapid_feature"
       Column       |          Type           | Modifiers
--------------------+-------------------------+-----------
 uid                | text                    | not null
 geom               | geometry(Geometry,4326) |
 bbox               | geometry(Polygon,4326)  |
 properties         | jsonb                   | not null
 create_timestamp   | time without time zone  |
 hash               | text                    |
 modified_timestamp | time without time zone  |
 layer_id           | text                    |

However, when I create an FeatureType based on an SQL View of this table, the properties attribute disappears:

Not sure what I need to do to make the properties column appear as an attribute in Geoserver.


texts: Convert Lat/Long to UTM specifying which zone to use (Java)

texts:
I have a set of coordinates in decimal degrees spread across two UTM zones in my Java application that I want to export to use in Petrel, so they need to be in the same UTM zone.

Is there a java library out there where I can specify what zone I want to use for the outputdata? tried quite a few libraries now and they only have utils for converting to the correct UTM zone.


texts: Please please please help me figure out what the deal is with my area units!!When calculate geometry for census tracts in km^2, I am given back some absurd, meaningless number such as 0.000696186


texts: Creating "Inner Buffer" using ArcGIS for Desktop?

texts:
I am currently searching for how to create an "Inner Buffer" in ArcGIS for Desktop?

I have a lot of polygons (houses) in one feature class, and one big polygon (in another feature class).

I want so select all houses contained in that one big polygon within a distance to the outer border (e.g. 1000 meters).

Select Layer By Location could be one possibility to solve this but I can't figure out how.


texts: Adding New Value for Field In Attribute Table prevented?

texts:
I have a polyline feature class in ArcGIS (version 10.3.1). There are only four options for the field "Type" (five if you include Null):

The symbology for this feature class was imported from a template feature class. I now want to add a new value for "Type" - "Footpath Redirected". I've added this to the "Unique Values" in the "Symbology" of the "Properties" of the feature class:

However, when I try to change a feature's type in the attribute table, I still have only the original four options (five if you include Null), as per the first screenshot.

What is the error in my approach?


texts:
I'm looking for a plugin that can make a route from given GPS points and then calculate the distance of this route.
I'm going to create GPS points and import it to Qgis as a GPX file.
I know I can calculate distances/polylines alone but I wanted to know if there was a plugin that did both.


texts: Suggestions for how to secure a custom ArcGIS extension using a license file or install key

texts:
I'm looking for a simple methodology for licensing custom ArcGIS extensions with an installation key or license file.  When the user attempts to enable the extension I want to check that they have a valid license before enabling the tools.   I’m not sure if this is the best forum for this question but I bet some of my fellow GIS developers have tackled this before.

Using an installation key would be the preferred approach because it seems less hassle, both for me and for the user.  I’m sure there are many of ways to accomplish this, but I’m trying to work out the basic process and I have a number of questions about it.

I’m assuming that once the user enters the key during the installation a value gets written to the registry.   The registry value is then checked to enable the extension.   Is there a better way?
What are the options for generating the installation key?  What rules can you create to determine a valid installation key vs. an invalid installation key?
What value should you write to the registry?  Perhaps encrypt the installation key and write the encrypted value to the registry.  Then when checking that key to enable the extension, decrypt the registry key, and then make sure it meets the same requirement of a valid installation key?  I’m not looking for Fort Knox here, but it seems you want the registry key value to be something not completely apparent to someone snooping around.

I've written extensions for ESRI software for years, (back to ArcView 3), without concern for validating licensing.  Usually the extension goes to a handful of clients with an associated license agreement text file, and that was it.  And, honestly, I prefer it that way.  If someone wants to break a licensing agreement they will find a way, so why waste time and energy trying to stop it.  However, for a variety of reasons I now need to apply a basic licensing scheme on a few extensions.  So, any suggestions or guidance would be greatly appreciated.


texts: Trying to make many polygons from points?

texts:
I am trying to turns points into polygons, but I encounter this challenge :

My data set has 1000 entries. Each entry has 360 associated columns, each with a coordinate pair. When connected, these coordinate pairs should all draw a specific circle or circle-like shape. So, the end result I am looking for is 1000 polygons, one for each entry.

It's a fairly large set to work with, so I am trying to strategize the best approach. What would you do?

It doesn't seem feasible to use Display XY Data --> Point to Line --> Line to Point for each entry. I'm looking for a more elegant solution.


texts: Converting File Geodatabase feature classes to shapefiles in ArcGIS Desktop?

texts:
How can I convert a File Geodatabase feature class into a shapefile?

I found the help for ArcGIS Desktop 10.1 and it tells me:

  On the main menu, click Customize > Toolbars > Production Mapping

but when I follow this menu, Production Mapping does not exist.


texts: How to correct link a set of colors to a rasterLayer

texts:
I have downloaded a raster with land cover data from http://www.diva-gis.org/gdata.

library(tmap)
library(raster)

grid_cov <- raster('VEN_cov/VEN_cov.grd')
grid_cov

class       : RasterLayer
dimensions  : 1416, 1656, 2344896  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -73.5, -59.7, 0.5, 12.3  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84
data source : /Users/altons/Dropbox/datacamp/geospatial course/VEN_cov/VEN_cov.grd
names       : VEN_cov
values      : 1, 22  (min, max)

Values in raster are:

unique(values(grid_cov))
20 15 14  2 13 16 12  1 19 22  7 17 18  8

which according to this document(pages 19 & 20)represent the type of land cover.

Now I have tried to map specific colors to land cover codes but my attempts have failed so far:

Attempt 1: pass color in ascending order as per land cover codes

cover_palette <- c("#008B00", "#00CD00", "#7CCD7C", "#98FB98", "#FFEFDB", "#FAEBD7", "#A4D3EE", "#CDC673", "#FF7256", "#FF7F50", "#EEE8CD", "#1874CD", "#F0F8FF", "#CD2626")

tm_shape(grid_alt) + tm_raster(legend.show = F,palette=cover_palette,style='cat' )

Bits in red/blue are mountains... which is wrong based on my definition of colors.

Attempt 2: Pass colors in the same order as per values(grid_cov)

cover_palette <- c('#F0F8FF',
'#CDC673',
'#A4D3EE',
'#00CD00',
'#FAEBD7',
'#FF7256',
'#FFEFDB',
'#008B00',
'#1874CD',
'#CD2626',
'#7CCD7C',
'#FF7F50',
'#EEE8CD',
'#98FB98')

tm_shape(grid_alt) + tm_raster(legend.show = F,palette=cover_palette,style='order' )
save_tmap(filename="~/Desktop/attempt2.png",dpi=75)

The last plot - the red bits are mountains which I thought I have colored in white....

Anyway - I think my question is: If I created a data.frame with all the colors then how can I merge it with the RasterLayer object? and therefore use colors$color instead of guessing!

my data frame looks like this:

lc_id = c(20,15,14,2,13,16,12,1,19,22,7,17,18,8)
cover_palette <- c('#F0F8FF','#CDC673',
'#A4D3EE','#00CD00',
'#FAEBD7','#FF7256',
'#FFEFDB','#008B00',
'#1874CD','#CD2626',
'#7CCD7C','#FF7F50',
'#EEE8CD','#98FB98')

colors <- data.frame(id = lc_id,color=cover_palette)


texts: Difference between processing.runalg() and general.runalg()

texts:
In the discussion from a question I asked several days ago (Using the processing framework in a PyQGIS standalone application), someone suggested that I use general.runalg() instead of processing.runalg() for running QGIS algorithms.

I've now used both, and haven't noticed a difference in either performance or syntax between them.  I'm just wondering what sort of differences there are between these two methods, and which one I should be using.


texts: Writing a custom GeoServer WPS process or Data Store

texts:
I have a third party custom written database that I want to access where I can only return JSON from queries (with a location lat/lon) .  I want to be able to visualize the results on a map. It will have minimal styling associated with it.

Could someone please guide me whether I need to write a custom data store or can I do this with a WPS process?

If possible I'd really like to write this in Python and utilize Geoserver Scripting?  Any tips/tutorials on this?


texts: How to download Planet Labs images for a large number of scenes efficiently?

texts:
I would like to download all of the reasonably high-quality Planet Labs images for each of about 200 areas of interest. What is the best way to accomplish this task?

I would like to use the Python API rather than the online interface for the sake of efficiency. However, it does not seem practical to do so using the approach provided in the API v1 Cookbook (https://www.planet.com/docs/api_cookbook/intro_part_3_download_aoi/) because of the need to activate the images and periodically check their availability.


texts: How can I delete an entry of a Geopackage

texts:
I created a geopackage in QGIS 3.18 that contains different entries (layers).
Now I have a very simple question (I hope): How can I delete just on entry of this Geopackage?
Cheers,
Markus


texts: Log or skip row causing error in PostGIS

texts:
In my simple PostGIS query to find the lengths of lines that fall within buffers (polygons)

SELECT a.name, b.route_id,  ST_LENGTH(ST_Intersection(a.buffer, b.wkb_geometry))
FROM buffered a, lines b
WHERE ST_Intersects(a.buffer, b.wkb_geometry)

The query runs endlessly and when I send a cancel request, I am returned an error:

  ERROR:  Error performing intersection: InterruptedException: Interrupted!`

This query has worked with one table of lines, but does not work with another. How can one log the index of the problematic row(s) in order that they be later skipped or removed?


texts: Is it possible to clip a shapefile to an image() in R?

texts:
I am working on mapping some netcdf data in R and have been having a tough time trying to perform a clip with an imported country shapefile that I'm using. Here is my question from yesterday. If you know anything about that, feel free to pop over there and answer that question. :)

BUT I'll start fresh here. So for my work, we are planning on mapping tons and tons of weather data all across the globe. For practice, the guy providing our data provided us with some "example" data to play with. My boss would like for me to create my maps in R. These maps are intended to be return period maps. The data is in netcdf format. I was able to map the netcdf data with the levelplot function from the rasterVis package. I was also able to map the netcdf data with a simple image function. Here's my code (sorry for the funky naming conventions; I've just been playing around with this):

library(maptools)
library(ncdf4)
library(RColorBrewer)
library(raster)

setwd("D:\\stuff")

#Madagascar shapefile
madagascar <- readShapeSpatial("MDG_adm0.shp")
plot(madagascar)

#Call on netcdf file
ncdf.data <- nc_open("swio_rpmaps_200_83.nc")

#Call on the 10 year wind return periods
returns <- ncvar_get(ncdf.data, "y010")

#Define latitude and longitude
lon <- ncdf.data$dim$longitude$vals
lat <- ncdf.data$dim$latitude$vals

#plot image of windspeeds
image(lon, lat, returns, col = cm.colors(9,alpha=.6), add = TRUE)
#plot Madagascar on top of windspeeds
plot(madagascar, add = T)

Here is the result image:

Is there any way to clip the Madagascar shapefile to this image? Or another way to manipulate the data so that I can perform a clip? Ideally, I would only have the data associated within the bounds of whatever country I'm mapping.


texts:
I have an annotation layer where each annotation has it's own font size, position etc.

Accidentally in ArcCatalog I deleted Symbol with an ID 0. That caused annotations to display as red rectangles. Now if I update SymbolID column with existing Symbol ID then annotation is displayed, but it automatically updates most columns in that row (for example FontSize). I absolutely need to keep existing values.

Is there any method of defining Symbol with ID=0 or to update SymbolID column without auto updating other columns?

I tried SymbolID -1, but ArcMap will not let me.


texts: Use svg file in openlayer to represent a layer

texts:
I want to use an SVG file as a source for one of the layers for my map in OpenLayers. The svg file is in a coordinate system such that it overlays my coordinate system of my map. I have tried the code below but it gives me

  AssertionError "Assertion failed. See
  https://openlayers.org/en/v5.3.3/doc/errors/#7 for details.

let Tunnel = new vectorLayer({
          source: new VectorSource({
            projection: projection2326,
            url: 'assets/Tunnel.svg',
            wrapX: false
          }),
          extent: [813967.0186951972, 825986.8061076261,814705.0826232458, 817842.6498878849]
        });

 this.map = new Map({
      target: 'map-container',
      layers: [
        new TileLayer({
          source: new OSM({
          wrapX: false
        })
      }), Tunnel
      ],
      view: new View({
        projection: projection2326,
        center: getCenter(extent),
        zoom:13,
      })
    });


texts: QGIS export vector PDF/SVG and import to Scribus: suspect map frame size

texts:
I would like to prepress my maps in Scribus (I need to convert them in to CMYK-Colorspace for my printing agency). So I need to export my maps as vector file from QGIS and import it into Scribus. So far so good. The vectors and labels are imported very well but the map size and the map frame has changed. You can see it in the second screen shot:

The red outside frame is my DIN A4 page from Scribus, the other red frame (selecting frame) is that one of my imported qgis-map. The black frame is my final map. It's scaled. But why? I think it could have something to do with the masking...

(Importing into Inkscape works well, but the "inkscape-reexport" behaves like the qgis-pdf.)

EDIT: I found out, that it makes a difference, how you import the PDF into Scribus. If you use File -> Import -> Get Vector File... or use drag and drop, the result will be like I described, but if you use  right click -> Paste File the imported PDF will placed in it's original size and borders. For me its enough to go further...


texts: What is the meaning of TL in LiDAR point clouds data sets?

texts:
I am exploring some USGS LiDAR data sets and found several occurrences of data sets with the same name except that one had the acronym "TL" in its name and the other not. Example of this is:

USGS Lidar Point Cloud NM RioSanJose 2016
USGS Lidar Point Cloud NM RioSanJose TL 2016

What is the meaning of "TL"? And how will it make the point cloud different from the the other data set?

My online searches gave nothing and the only possibility I can think of would be Topographic Lines but I could not confirm so far.

Edit: The links do not point to the same tile, but more examples are found in this FTP


texts: Line width and labels based on blended values of features sharing same geometry in QGIS

texts:
I want to appropriate line segments that share the same geometry with a particular line width and labels based on the aggregated values of one field "GIS@BUS_5_Min".

My attribute table and polyline objects look as following, see image below

"BUSLINIE" means a bus line, "Overlap" means if the line has overlaps, i.e. it's geometry repeats.

Values of the field "GIS@BUS_5_Min" for these three objects has to be aggregated as one value and represented as a label sum("GIS@BUS_5_Min").

And the width of this line (selected in yellow) will be sum("GIS@BUS_5_Min")/100 = 19.10/100 = 0.1910

How can I achieve such output?

I assume it is possible with Layer Properties > Width > Data defined override > Edit..., but I do not know how to check the same geometry and how to sum the values of certain features with CASE WHEN condition THEN result END conditional function?

Dissolve function does not work in my case because I do not have any key-fields.

References:

QGIS Line width transition
Blend differing line thicknesses in QGIS?
Changing line width by classification in QGIS
How to sum up values of overlapping polygons in QGIS?
Merging attribute and geometric features in QGIS?


texts: Import self-written ArcView 3.2 extension to ArcMap 9.3 or higher

texts:
I've found an older self-written extension (.avx) for ArcView 3.2 from a retired co-worker. Unfortunatly I'm not a programmer, but I'd like to use this extension for ArcMap 9.3. or 10.x.
Is there any way to import/convert this extension?


texts: Querying the currently visible layer in Openlayers

texts:
My map has several layers, and I have a check box list that can turn individual layers on and off (ie. hidden/visible).

I'd like to have querying capability, but with so many layers, I think the easiest way might be to only allow querying on whatever layer is visible, and if there are a few visible, only the top layer is queryable.

The tutorials I've seen are only for querying one layer. How would I go about doing this for multiple layers and for only those that are visible? This is a poor attempt, and only works for one layer:

  var displayFeatureInfo = function(pixel) {
    var features = [];
    map.forEachFeatureAtPixel(pixel, function(feature, layer) {
      features.push(feature);
    });
    var container = document.getElementById('information');
    if (features.length > 0) {
      var info = [];
      for (var i = 0, ii = features.length; i < ii; ++i) {
        info.push(features[i].get('County_Name'))+info.push(features[i].get('Name'))
      }
      container.innerHTML = info.join('- ') || '(unknown)';
    } else {
      container.innerHTML = '&nbsp;';
    }
  };

        map.on('click', function(evt) {
    var pixel = evt.pixel;
    displayFeatureInfo(pixel);
  });


texts: How to send multiple features in the SOE response for the multi point input locations?

texts:
For multi point input locations, I need to send multiple features in the SOE response. How do we combine the multiple byte[]

result = FindNearFeatures(layerID, location, distance.Value);


texts: Removing unused references in ArcObjects vba?

texts: Can anyone suggest me how will I find and remove unused references used in a project done in ArcGIS 9.3.1 vba.


texts: QgsProcessingParameterString(), Store user input?

texts:
I'm adding some stuff to a QGIS plugin and trying to figure out how to work with user inputs. In this case, the user input is a string, which seems simple to work with, but I can't find an example of it anywhere.

I am not using Qt Designer to design my interface for the dialog box, I'm instead using various methods from the QgsProcessingParameters class. Specifically I am using QgsProcessingParameterString.

In the plugin, I'm using the initAlgorithm and adding a line for the user to input a string. Here is the code:

def initAlgorithm(self, config):

#Add String as input
    self.addParameter(
        QgsProcessingParameterString(
            self.INPUT_TITLE,
            self.tr('Input title string (Optional)')
        )
    )

Now in the processAlgorithm I use a method .parameterAsString

def processAlgorithm(self, parameters, context, feedback):
    """This actually does the processing for creating the print layout and exporting as .pdf"""

    log = feedback.setProgressText

    input_title = self.parameterAsString(parameters, self.INPUT_TITLE, context)

    log(print the user input)

I'm just trying to print the user input in this example.
How do I get the value the user input (a string)? I want to create a variable which stores the user input as a string. I can use the string later in my program.

And are there any examples of using QgsProcessingAlgorithm.parameterAsString?


texts: Can gdal_calc operate conditionally with values from two tifs

texts:
I am trying to correct some water table .tifs (WT) to the ground level elevation (DEM) so that water levels are not incorrectly presented as above ground. I have looked at many similar questions and answers but it seems to me that the conditional capability of gdal_calc is limited to subbing in a static floating point number. I need it to use values from a DEM in the event that the water table is above ground. Pythonically, this would be

if (DEM - WT) > 0:
   output = WT
else:
   output = DEM

However, I can only seem to get gdal_calc to sub in a static floating point number, not a reference to a more dynamic .tif file where the value would depend on the grid location.

Is this possible with gdal?  I have tried the following command line calls:

python "C:\Program Files\GDAL\gdal_calc.py" -A DEM.tif -B WT.tif --calc="B*((A-B)>0)" --outfile WT_belowGround.tif

This produces a .tif with values only where WT is below ground with No Data gaps where it is above ground. This is close to what I want but I want data coverage in all locations for producing contours.

python "C:\Program Files\GDAL\gdal_calc.py" -A DEM.tif -B WT.tif --calc="B*((A-B)>0)" --calc="A*((A-B)<0)" --outfile WT_aboveGround.tif

This produces a .tif with values only where WT is above ground with No Data gaps where it is below ground. Again, this is close to what I want but I want data coverage in all locations and it seems the second --calc call overwrites the first with No Data values.

Trying to combine both of these I have tried several things without the desired result. The most obvious thing to me is to assign the NoDataValue to the desired DEM.tif as is shown below.

python "C:\Program Files\GDAL\gdal_calc.py" -A DEM.tif -B WT.tif --calc="B*((A-B)>0)" NoDataValue=DEM.tif --outfile diff.tif

While this runs and produces a .tif file, it gives an error and again, it produces a .tif with values only where WT is below ground with No Data gaps where it is above ground. The error is below

RuntimeWarning: invalid value encountered in greater

Is this something that gdal can do or do I need to find another way to approach this? One thought I have is to produce the two .tif files of WT above ground and below ground and combine them but I am not sure how to combine them either.


texts: Maximum likelihood classification, Landsat: Why to exclude thermal band?

texts:
(transfered from https://stackoverflow.com/questions/39348977/maximum-likelihood-classification-landsat-why-to-exclude-thermal-band)

I search for an argument (which I could cite, ideally) to support my decision to exclude Thermal band 6 from Maximum likelihood classification (MLC) of Landsat (5-7) imagery.
To my knowledge, the thermal band 6 is suggested to exclude from MLC because of its coarser spatial resolution (~ 120 m), comparing to another bands (30 m).

Can you suggest me some papers resulting this?

The only one I've found is:
http://www.isprs.org/proceedings/XXXIV/part1/paper/00077.pdf
... a classification utilizing all of the TM bands
(excluding the thermal band due to its coarse spatial resolution) ...

Do you have suggestions for another paper/handbook?


texts: Select Layer By Location gives ERROR 000840 ({select_features} as arcpy.Polygon)

texts:
I tried to select SC27060054 point object using polygon.
r"Project\project.DBO.Qo" is a point layer.

import arcpy

p1 = arcpy.Point(254426.380,  4664270.258)
p2 = arcpy.Point(254432.730,  4664270.258)
p3 = arcpy.Point(254432.730,  4664264.967)
p4 = arcpy.Point(254426.380,  4664264.967)

poly = arcpy.Polygon(arcpy.Array([p1, p2, p3, p4]))

arcpy.SelectLayerByLocation_management(r"Project\project.DBO.Qo", "CONTAINS", poly)
Runtime error  Traceback (most recent call last):   File "<string>", line 1, in <module>
File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 7279, in SelectLayerByLocation
raise e ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Boolean. Failed to execute (SelectLayerByLocation).

desc = arcpy.Describe(r"Project\project.DBO.Qo")
print desc.shapeType
Point


texts: Modifying Font Style simple marker according to character value

texts:
Currently I am writing a standalone Python script and I am trying to modify a simple point marker using font style. There are some interesting fonts with different characters that aren't exactly on a normal keyboard. Does anyone know how to use these characters?

For example, I currently have:

symbol = QgsMarkerSymbolV2.createSimple({'color': 'transparent', 'outline_color': 'transparent'})
layer.rendererV2().setSymbol(symbol)
# create the fontStyle symbol
fontStyle = {}
fontStyle['color'] = 'gray'
fontStyle['font'] = 'DejaVu Serif'
fontStyle['chr'] = '0x65'
fontStyle['size'] = '8'
# create a layer for the fontStyle
symLyr1 = QgsFontMarkerSymbolLayerV2.create(fontStyle)
symbol.appendSymbolLayer(symLyr1)

I can't type in the actual character in fontStyle['chr'] but I do know the actual value pertaining to the desired character. Is there a 'value' property such that I can say something like:

fontStyle['value'] = '0x65'


texts: arcpy.QuickExport_interop crashing python when exporting to gml

texts:
arcpy.QuickExport_interop is crashing python when exporting to a gml.  We are running ArcGIS 10.4.1 on a Windows Server 2012 R2.  The same script runs successfully on two different desktops.

Problem signature:
  Problem Event Name:   APPCRASH
  Application Name: python.exe
  Application Version:  0.0.0.0
  Application Timestamp:    5560ae54
  Fault Module Name:    arcgisscripting.pyd
  Fault Module Version: 10.4.0.5524
  Fault Module Timestamp:   567b1378
  Exception Code:   c0000005
  Exception Offset: 000000000008c1b0
  OS Version:   6.3.9600.2.0.0.272.7
  Locale ID:    1033
  Additional Information 1: 6fd2
  Additional Information 2: 6fd24e2e58701489213ba4c83af23142
  Additional Information 3: 652c
  Additional Information 4: 652c12f5220101ef60f75474ab1ca3b6


texts: Spatial Join of 2 line layers offset by minuscule amount

texts:
Let's say I have 2 line layers (green & red, in the image below) and I need to transfer attributes from one layer (red) to another (green) by doing a Spatial Join. The layers are identical, but offset by a minuscule distance, b/c of some reprojection & geographic transformation -- something like 0.005 meters (5 mm).

Of course this is not even a meaningful offset for any real-world analysis issue, but for a Spatial Join, it seems it would prevent the lines from actually overlapping perfectly & thus the Intersect operation would transfer the wrong attributes.

Using "Within a Distance" w/ some similarly tiny distance would invariably capture not only the parallel green line that the red line was meant to be on top of, but the perpendicular green line as well.

It's hard to believe I'm really out of luck on such a simple operation & have to resort to some much more complex procedure, simply due to a less than 100% absolutely perfect overlap between 2 otherwise completely identical line layers. And I can't really imagine any workaround. Any ideas?? Thanks very much!

Image below!:


texts:
I have a data set of roads which is line and its projection system is WGS UTM (EPSG:32639) I have added it into my database and created a topology using this query

SELECT pgr_createTopology('ways', 0.00001, 'geom_2d', 'gid',rows_where := 'gid>0', clean := TRUE);

The above query works fine. But When I try to route using any alogorithm such as this query

    SELECT id1 as path,b.geom as the_geom
FROM pgr_dijkstra('SELECT gid as id, source, target, st_length(geom_2d) as cost FROM ways',
                           483,507,FALSE ,FALSE) a,
ways b
WHERE a.id1=b.gid

It returns some lines which are not a route at all. Here is an example

In above image the yellow lines is my edges, points are vertexes which are created by pgr_createTopology and blue lines is the result of routing
I have tested different values of tolerance such as 1,10,100,0.00001,0.001 and for some of them there is no result and for some f them the result is exactly the same, I think there is something wrong with layers topology but I can not figure out how to fix it.

Can you help me to fix my problem?

edit 1:
here is the result of this query

SELECT * FROM "ways" WHERE gid=483 or gid=507


texts: Using UpdateCursor for joined field in ArcPy?

texts:
I have a huge point vector data (around 3 million point; in a gdb) with attributes and I have table (created with Near analysis.). I want to copy one of the field values to data. So I joined the table and tried the field calculator to do the trick, but after 16 hours I realized it just froze (it worked with similar data).
So the next that could work I think is the UpdateCursor. Unfortunately I got an Error: "TypeError: cannot update join table"

tried this:

field1 = '!temp_Features.leftd!'
field2 = '!leftt:NEAR_DIST!'

with this:

cursor = arcpy.da.UpdateCursor("layer2",(field1,field2))

this:

fields = ['temp_Features:leftd', 'leftt:NEAR_DIST']

with this:

cursor = arcpy.da.UpdateCursor("layer2",fields)

Same error every time.

(I have ArcGIS with advanced license)


texts:
I have been working in Python with the Acaconda Python distribution to create code that accomplishes the following algorithm:

I have multiple variable netCDF4 files [NCEP Reanalysis tmax(K), tmin(K), shum(%); prate(kg/m^2/s) daily values for an entire year(s)]. Each file spans between (Jan 1 2006- Dec 31 2010). They are 16(lat) by 16(long) grids. Data was downloaded from ESRL: PSD: NCEP Reanalysis. The dimension is time in days. What I'm trying to accomplish is if any given day (n) that satisfies conditions for all variables at each corresponding lat/lon:

tmax and tmin: 18°C ≤ T(n) ≤ 32°C [conversion to K: 291.15K ≤ T(n) ≤ 305.15K]

shum(Q): 20 ≤ Q(n) ≤ 80

prate(R): 1.5mm ≤ R(n) ≤ 20mm [conversion to (kg/m^2/s): R(n) / 84600]

I want to store the number of days PER YEAR in a 4 new netCDF4 - or ArcGIS 10.1 raster/compatible - files (one for each year).

From my understanding, I have not been able to find the correct function to loop through various time steps at specific lat/long/time variables. Perhaps, I have not used the correct phrasing in my search, but I am a novice in programming beyond simple routines.


texts: QGIS Raster calculation fails in Graphical Modeler

texts:
I'm trying to use the graphical modeller to automate a certain process (extracting North-South facing slopes from a DEM).

The process is as follows:

Calculate Slope
Calculate Aspect
Extract North-South Slopes using the raster calculator expression:

Slope*cos((3.141592/180)*Aspect)^2

This works flawlessly when I do it step by step, using the default Processing toolbox raster calculator, Slope and Aspect.

When I try to input this in the Graphical modeller however I get an error: "pointer cannot be cast to type 'float'"

The aspect and slope are calculated as expected, but the raster calculator fails.

I've tried this in QGIS 3.8 and 3.10.

I've tried assigning the CRS to the slope and aspect layers before the calculation, and tried using the GDAL calculator (which returns a "TypeError").

My question is: What am I doing wrong here? The same steps work perfectly outside the graphical modeller, so I don't see why they should fail here.


texts: How to maintain connection to ArcGIS Online through ArcPy without having to have the ArcMap application open?

texts:
I am working with ArcMap 10.5 on a corporate network though a VDI. We have been trying to implement a batch script which will pull data from an API to update features in a gdb, which then updates the souces in an MXD. This all works just fine. However, the script then tries to automate publishing the MXD to overwrite a featureservice in an ArcGIS Online developer account, and even though we are passing our login credentials to the arcpy.SignInToPortal_server function, we get this error message:

Traceback (most recent call last):
  File "(path to python script)/UpdateDashboardService.py", line 95, in <module>
    "PUBLIC", "SHARE_ORGANIZATION", "")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\server.py", line 1509, in UploadServiceDefinition
    raise e
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Server: Dataset My Hosted Services does not exist or is not supported
WARNING 001404: You are not signed in to ArcGIS Online.
Failed to execute (UploadServiceDefinition).

After doing some searching online we found a workaround where if you have an ArcMap document open through the GUI and signed in to the account and then run the batch, the stand alone script will work and the features will be updated to the developer account. This does work, however we really are trying to make this completely standalone and not require more than just running the batch the script is part of so that this process is as seamless as possible. Requiring the user to manually open ArcMap and sign in to the service is not ideal.

I'm wondering if there is an ArcPy function or tool that will force the ArcGIS Connection Utility to display we are logged in without having to have the MXD open.


texts: pyqgis: Get the Click on canvas standalone app

texts:
I added a subwindow in a MDI area containing a QgsMapCanvas.
Defined the clic tool QgsMapToolEmitPoint.

But when I click on the canvas there's no response, al least it should print 'clicked'.
¿¿May be something is missing??
(debug with pycharm on windows)

from PyQt5.QtCore import *
from PyQt5.QtGui import *
from PyQt5.QtWidgets import *
from qgis.core import *
from qgis.gui import *
from qgis.utils import *

from resources.main2 import *

class myMainClass(QtWidgets.QMainWindow, ):
    def __init__(self, parent=None):
        super(myMainClass, self).__init__()
        self.ui = Ui_MainWindow()
        self.ui.setupUi(self)

        self.swMap = QgsMapCanvas()
        self.swMap.setCanvasColor(Qt.white)
        self.swMap.enableAntiAliasing(True)
        self.swMap.setSelectionColor(QColor(255, 255, 0, 50))
        self.ui.mdiArea.addSubWindow(self.swMap)

        pointTool = QgsMapToolEmitPoint(self.swMap)
        pointTool.canvasClicked.connect(self.changeColor)
        self.swMap.setMapTool(pointTool)

    def changeColor(self,point, mouse_button):
        print('clicked')
        if mouse_button == Qt.RightButton:
            self.swMap.setCanvasColor(Qt.Red)
        if mouse_button == Qt.LeftButton:
            self.swMap.setCanvasColor(Qt.Blue)

if __name__ == "__main__":
    print('start')

    app = QgsApplication([], True)
    app.initQgis()
    myapp = myMainClass()
    myapp.show()
    sys.exit(app.exec_())

The class main2.py

from PyQt5 import QtCore, QtGui, QtWidgets

class Ui_MainWindow(object):
    def setupUi(self, MainWindow):
        MainWindow.setObjectName("MainWindow")
        MainWindow.resize(258, 137)
        self.centralwidget = QtWidgets.QWidget(MainWindow)
        self.centralwidget.setMinimumSize(QtCore.QSize(0, 0))
        self.centralwidget.setObjectName("centralwidget")
        self.verticalLayout_2 = QtWidgets.QVBoxLayout(self.centralwidget)
        self.verticalLayout_2.setObjectName("verticalLayout_2")
        self.mdiArea = QtWidgets.QMdiArea(self.centralwidget)
        self.mdiArea.setObjectName("mdiArea")
        self.verticalLayout_2.addWidget(self.mdiArea)
        MainWindow.setCentralWidget(self.centralwidget)
        self.menubar = QtWidgets.QMenuBar(MainWindow)
        self.menubar.setGeometry(QtCore.QRect(0, 0, 258, 21))
        self.menubar.setObjectName("menubar")
        MainWindow.setMenuBar(self.menubar)
        self.statusbar = QtWidgets.QStatusBar(MainWindow)
        self.statusbar.setObjectName("statusbar")
        MainWindow.setStatusBar(self.statusbar)
        self.actionOpenSubWin = QtWidgets.QAction(MainWindow)
        self.actionOpenSubWin.setObjectName("actionOpenSubWin")

        self.retranslateUi(MainWindow)
        QtCore.QMetaObject.connectSlotsByName(MainWindow)

    def retranslateUi(self, MainWindow):
        _translate = QtCore.QCoreApplication.translate
        MainWindow.setWindowTitle(_translate("MainWindow", "MainWindow"))
        self.actionOpenSubWin.setText(_translate("MainWindow", "OpenSubWin"))


texts: Creating a table of multiple rows and columns from attribute table

texts:
I am trying to create a street index a City. I have a table that has all the street names and their grid references. I am looking to create a street index using that. I think turning off the undesired fields off leaving only the main three columns I need and then outputting that would work. That instead  creates a really long table. I would like to break that table into columns so I can add it to street maps. In the end I need an index that looks like this.

That is the end goal. If there is a way through python I am willing to attempt that as well. That right now is created by polylines and a white polygon in the background and an annotation layer for the text. All done by hand. I need to have it automated so once a new input is made the user does not have to create and enter in a new annotation and create the correct size of a table.

I am using ArcGIS 10.3 on a Windows 7 computer. Let me know if any additional information is needed.


texts: adding x and y into single column in QGIS

texts:
I have always use the field calculator to retrieve my x and y coords ($x, $y), but this ends up with 2 separate columns, which I then have to manually concatenate.
Is there a way to do this easily within QGIS? - perhaps using the expression box? I feel like there must be some easy expression similar to concatenate in excel (e.g. $x, ", ", $y).


texts: QGIS3 How to export shapefile symbology to scale?

texts:
I'm looking how to export in QGIS3 a vector layer with a determined scale.
It was straightforward in QGIS2 when I was exporting to a shapefile I had the option symbology scale but it disapear in QGIS3 or maybe I'm missing something.
The intent is to use the shapefile in another program.


texts:
I notice that a GPX file contains not only lat + lon but also elevations:

  <trk>
    <name>Some route</name>
    <trkseg>
      <trkpt lat="50.64536" lon="3.05657">
        <ele>23.6</ele>
      </trkpt>
…
      <trkpt lat="50.63836" lon="3.06491">
        <ele>32.7</ele>
      </trkpt>
    </trkseg>
  </trk>

I was wondering if someone knew of a web site where I could simply upload a GPX file and have it drawn on a map with the route colored differently based on the elevation, ie. where the flat parts are drawn in green while the steep parts are shown in red?


texts: QGIS Print Layout Grid Reference Numeric Format

texts:
Using QGIS3.14
The default style of grid labels is single size digits, e.g. 251000E and 5897000N.
I wish to show the grid labels in a format similar to Ordnance Survey with the 51 of the easting and the 97 of the northing in a larger size.
Can this be done within the Print Layout Grid properties or is there a plugin to achieve this?


texts: Setting correct coordinate system when rxporting from Draftsight

texts:
I'm trying to import a DXF into GIS (I have QGIS, Arcmap 9.2, GVsig )
The file has been given a coordinate system & appears to be the right size, just in the wrong place.
In Draftsight (or any other software that can read DXF) how can I set a correct coordinate system?  I've been looking at the Draftsight resources but can't work it out as I have limited CAD knowledge.

I can't really affine this in QGIS as I have no base mapping, but I am sure I know the existing CS is out by 400,000m on x & 100,000 on y (OSGB Tile code SU) so if I could do this in Draftsight & re-save it it would come out correct.


texts: Unexpected slope returns from terrain analysis

texts:
I'm writing with a question regarding unexpected returns on a slope analysis using DTM LiDAR data at the scale of 2m in QGIS.

This is for my MSc dissertation (UK) which is due in about 6 months, so I want to be thorough.

I carry out the analysis I transformed the DTM raster data into slope raster data using QGIS slope Raster Terrain Analysis function - output format = GeoTIFF with Z factor of 1. I then converted the raster slope data into vector data (polygon) using the raster to polygon function. When I review some points on the map, they don't seem to align with what I know of the land (a few degrees off).

I should note that I have my map (including vector point data) in the WGS 84 Pseudo Mercator CRS whereas my DTM in data requires OSGB 1936 to overlap.

Perhaps I'm missing something, or are the slope outputs are correct?

LiDAR data was sourced from link provided below.
https://data.gov.uk/harvest/gemini-object/33e78ae5-27d9-44db-9e2c-2e75c049e09c


texts: Casting igxobject to ifeatureclass in ArcObjects C#?

texts:
In above Image text box show the file path. I selected file from the igxdialog and file stored in igxobject. Now to how to casting from igxobject to ifeatureclass of selected file... but igxobject is in button1 and validate button can't access igxobject..

code:

       private void btnbrowse_Click(object sender, EventArgs e)
    {
        ESRI.ArcGIS.CatalogUI.IGxDialog dlg = new ESRI.ArcGIS.CatalogUI.GxDialogClass();
        ESRI.ArcGIS.Catalog.IGxObjectFilter objfltr = new ESRI.ArcGIS.Catalog.GxFilterFeatureClassesClass();

        objfltr = dlg.ObjectFilter;
        ESRI.ArcGIS.Catalog.IEnumGxObject Sel = null;

        dlg.DoModalOpen(0, out Sel);
        ESRI.ArcGIS.Catalog.IGxObject obj = Sel.Next();
        txtdatapath.Text = obj.FullName;
    }
    private void button2_Click(object sender, EventArgs e)
    {
        System.IO.StreamWriter sw = new System.IO.StreamWriter(@"D:\ARC OBJECTS\nullerrors.txt");
        ESRI.ArcGIS.Catalog.IGxDataset dset = txtdatapath as ESRI.ArcGIS.Catalog.IGxDataset;
        ESRI.ArcGIS.esriSystem.IName nam = dset.Dataset.FullName;
        ESRI.ArcGIS.Geodatabase.IFeatureClass fcls = nam.Open() as ESRI.ArcGIS.Geodatabase.IFeatureClass;


texts: How to import label classes to make symbology classes?

texts:
In Answer to an earlier Question I was able to write SQL Expressions to create label classes.

Now, I wish to use those label classes to create symbology classes.

It is possible to make label classes from symbology classes, but I want to do it the other way around.


texts: Creating feature template in blank map of ArcMap

texts:
I understand how to create features with existing templates, but sometimes a map has no templates to begin with and no layers to use to create them; perhaps just a basemap.

Is there a way within ArcMap to do this, or must I create features in Catalog first?


texts: Garmin Communicator Plugin not working on any browser

texts:
When I access the garmin plugin compatibility site.

Which one? This one: http://software.garmin.com/es-ES/gcp.html

I get the following notice: [Browser name] No Longer Supports this Plugin.

On any browser.

This happens on different laptops.

I followed the instructions to enable the plugin on Internet Explorer, and it won't work.

Anybody else having the same issue?


texts: Best GIS system for high performance web application - PostGIS vs MongoDB

texts:
I am working on a web/mobile application based on location data. Since i am already familiar with MongoDB, i found the geospatial indexing of mongo is quite suitable for my needs. As i am mainly dealing with simple/short location points, Mongo 2d indexing is good for me.

Along the way i picked PostGIS, because of its stable/mature way. And its awesome feature set. But my main concern is performance since my data is heavily dependent on location(mostly 70 - 80% of the db calls deal with the location).

I like mongo because its used by high performance web apps like foursquare already. But i have seen PostGIS is mainly used in government/enterprise projects (mostly non web/mobile apps). So i am little confused right now to choose the right GIS db for my web/mobile app? Got any suggestions?


texts:
I am trying to create random locations nearby my location. What i want is to create random latitude/longitude pairs inside a 200 meters circle surrounding my location.

This is the formula i came up with (with the help of people at StackOverFlow):
(Random number between -1 and 1)*radius + (old longitude) = new longitude within radius of old longitude

(Random number between -1 and 1)*radius + (old latitude) = new latitude within radius of old latitude

The thing is that something weird is happening with  my implementation because all the random locations are too near of my location center, it seems that the formula does not cover the whole radius.

Any idea of what could be wrong with my formula?

Edited to show the current java implementation:

public static Location getLocation(Location location, int radius) {
    Random random = new Random();

    // Convert radius from meters to degrees
    double radiusInDegrees = radius / METERS_IN_DEGREES;

    double x0 = location.getLongitude() * 1E6;
    double y0 = location.getLatitude() * 1E6;
    double u = random.nextInt(1001) / 1000;
    double v = random.nextInt(1001) / 1000;
    double w = radiusInDegrees * Math.sqrt(u);
    double t = 2 * Math.PI * v;
    double x = w * Math.cos(t);
    double y = w * Math.sin(t);

    // Adjust the x-coordinate for the shrinking of the east-west distances
    double new_x = x / Math.cos(y0);

    // Set the adjusted location
    Location newLocation = new Location("Loc in radius");
    newLocation.setLongitude(new_x + x0);
    newLocation.setLatitude(y + y0);

    return newLocation;
}

I am not sure what i am doing wrong, because the new locations are created in the middle of the sea.

Any idea?


texts: Change Icon of Single Vector Feature in Layer

texts:
Currently I am trying to change the icon of a particular feature of a vector layer that the user is focusing on.  I add each feature to the map like so:

var point = new OpenLayers.Geometry.Point(pt.lon, pt.lat);
var markerStyle = OpenLayers.Util.extend(OpenLayers.Feature.Vector.style['default'],  {
    externalGraphic: iconURL
});
var marker = new OpenLayers.Feature.Vector(point, attributes, markerStyle);

Later I do the following to update the feature's icon:

var marker = this.findSelectedMarker();
if (marker) {
    marker.style.externalGraphic = newIconUrl;
    this.layer.redraw();
}

But when the layer redraws, all features in my layer use the newIconUrl, not simply the selected marker that I am trying to update.

How can I change the icon of the one selected feature of my layer?  Thanks.


texts: Installing LecoS plug-in for QGIS v3.16 on Mac Big Sur

texts:
I am new to QGIS and was learning some spatial analysis. I have tried to install LecoS plug-in on QGIS v3.16, unsuccessfully.
I have tried the All-in-one version, but when I install the plug-in I got an error saying that PIL is not in the library. I have tried to install it using KyngChaos Python modules (https://www.kyngchaos.com/software/archive/pymod/) but without success.
then
I have tried the alternative build of KyngChaos (https://www.kyngchaos.com/software/qgis/), here a part form many crashes especially when opening Python console, I got an error stating that Scipy was missing in QGIS pathway. Then I have tried to install Scipy but it goes under the directory Python 2.7 keeping the former error as the alternative build work only on Python 3.


texts:
For every cross-section profile, Manning's value is assigned.

There are numerous values in a single cross-section profile but HEC-RAS allows only 20 values per cross-section.

How can I deal with this problem?


texts: Retrieving Aerosol Optical Depth (AOT) from Landsat multispectral data

texts:
I need to retrieve Aerosol Optical Depth from Landsat multispectral datasets.

Has anyone successfully tried something like that before?


texts: Using open source software to merge large shapefiles?

texts:
Some background: I have 2 polygon shape files with approximately 500,000 records each. The two files are identical in terms of field structure.

Using ArcGIS 10.0 I can'nt merge these two files using Geoprocessing / Merge. ArcMap allways crashes somewhere in the middle of the process.

Is there any open source software that can do this job? QGIS has the function: Vector/ Merging shape file into one , but there was a reported bug with this function in the current version 1.7. I know GDAL can do a lot of things and wonder if GDAL or some other opensource software can do this properly.


texts: List table names within PostGIS bounding box

texts:
So I have a bounding box drawn by the user on Leaflet from which i can get features of specific table as follows
 SELECT * FROM $table WHERE $table.geom && ST_MakeEnvelope($test, 4326) where $table is the variable holding table names selected by the user and $test is the BBOX coordinates.

This works fine, I get the list of features within the BBOX of that selected table.

Now if I have to List all tables in my database within the BBOX how would I proceed?


texts: Generation and Visualization of Stacked Maps in Qgis software

texts: We want to generate 3D stacked model of multiple 2D maps in  Qgis software. We installed the PcraterMapstakVisualisation pluggin on the software but havent been able to achieve this. Kindly assist us in achieving this objective in Qgis software


texts: How many bands can a single NetCDF file contain?

texts:
I am collecting some precipitation data (centimeters of precipitation per day) from the Gridmet dataset (http://www.climatologylab.org/gridmet.html). The data comes in a NetCDF format and there is a different band for each day of the year. So the entire raster has 365 bands or dimensions, and at each pixel there are 365 individual measurements for cm of precipitation. Notethat Gridmet has a bunch of data on other climate indicators in separate files.

My question was, how many bands or dimensions can a single NetCDF file contain? I was trying to google this, but did not find anything. Is there a limit on the number of bands (I call them bands because I am used to thinking in tif rasters), or is there a bound on the total amount of data in the file, and so I can have unlimited bands until the data max is reached?

The reason I was asking is that I want to consolidate a bunch of datasets into a single or handful of NetCDF files. So if I want precipitation for 20 years in a single NetCDF file, which would be 365*20 = 7300 bands or dimensions, would that work? That way I can pull the full time-series for any point. If someone has a better idea, let me know.

Any info is appreciated. Thanks.


texts:
I have around 40K raster layers, each having an average size of 300 MB. I am trying to build the raster visualization system in which the user can select some parameters from different dropdown buttons and get the required layers. My question is how can I upload/manage the store/workspace so that the performance of the system will remain the same? What happens, if I upload all the rasters into only one workspace?

Is the performance of the system will slow down due to huge rasters?


texts: How I can disable some languages in Geonetwork

texts:
Please, how I can display just some languages in Geonetwork?
I edit thes files, but I don't see any change in the dropdown menu to select the language, In fact , I want just those languages: French, English and Arabic, but actually I see all languages that are by default in Geonetwork.

config-gui.xml

    <languages>
        <ara/>
        <eng/>
        <fre/>
    </languages>

js/CatControllers.js:

    $scope.langs = {'eng': 'en', 'fre': 'fr', 'ara': 'ar'};
    $scope.langLabels = {'eng': 'English', 'fre': 'Français', 'ara': 'Arabe'};

web-inf/config-spring-geonetwork.xml

    <util:set id="languages" value-type="java.lang.String">
       <value>ara</value>
       <value>eng</value>
       <value>fre</value>
    </util:set>


texts:
Is it ok to respond with a single item inside a FeatureCollection or should I respond with a Feature object instead? I'm using FeatureCollection in a REST-like server, as index view (/things.geojson) as well as show view (things/3.geojson), like this:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates": [ ... ]
      },
      "properties": { ... }
    }
  ]
}

I couldn't find anything in the spec, so if it's undefined I'm asking for best practices and references.


texts: Storing and updating large amount of raster data

texts:
Large amount of data = ~10TB of rasters (multiple projections) with the possibility of this data growing.

The basic architecture I'm planning at the moment is having 2 clusters of GeoServer--one cluster (2 machines) as a WMS and another as a WFS. I'd like for these clusters to point to one repository for their data (raster and vector) and I am trying to figure out whether it's better to store my raster data in a file system or a database.

I've been reading "PostGIS in Action" and I'm leaning towards a DB structure for the rasters, but I'm still on the fence about it.

I have two main goals that the structure must be able to support:

Automate the process of building pyramids to load up onto GeoServer
Target specific tiles to update when new rasters come in

Right now I know I can use gdal_retile to manually build the pyramids--that coupled with GS REST API I can automate the uploading and building (although I am admittedly having trouble getting GeoServer-Manager running). But I feel a bit weary just leaving stuff in a file system given the large amount of data and the updating of files (probably going to be monthly, but could be more frequently). Because of this, the DB approach seems more appealing--but I am just a software engineer; GIS is completely new to me and I'd like some wisdom from people who do this before I start building things that can't scale properly.


texts: How to scale features declining using map units in QGIS

texts:
I have layer in QGIS that I want to style using SVG symbols.

I want the SVGs to increase in size when zooming in the map. I can solve this by using map units instead of millimeters.

Though, I don't want the symbols to increase 1:1.

Say, I want the symbol to be 75 map units at zoom level 1:50.000 but only 37.5 map units at zoom level 1:3.000.

I guess that using QGIS Expressions here is key, but I can't work out a proper mathematic formula that fixes my problem.

Can I use the scale_exp function in the Expression string builder in QGIS?

Can anybody help out?


texts: Clipping with grid and naming layers in QGIS?

texts:
I have created a regular (nominally square) grid using Create grid and can clip points and lines within each grid square using the Clip tool. So now I have 168 layers called Clipped. I'd like to give each layer a unique name (e.g. Cliiped1, Clipped2, ..., Clipped168 will do!).

Is this possible? How?


texts: Different scale of basemap after exporting to PDF/Image

texts:
After creating and exporting my map to pdf/tif obviously there are different scales of the basemap in the exported file compared to the map I created in the printcomposer. See pictures...

I am using QGIS 2.8.1 on a windows 7 (64bit) machine. For the basemap I use basemap.at (WMS). The CRS is the same for all layers (WGS 84 / Pseudo Mercator).

Created Map in Print composer at 1:950000:

And what I get after exporting it to pdf is a more detailed basemap...

Does anybody know why this happens?


texts: How send Geoserver WMS GetFeatureInfo query to an area

texts:
I am using GetFeatureInfo query to WMS service. When I click a WMS object on the map, the response is coming.

http://localhost:8080/geoserver/wms?
request=GetFeatureInfo
&service=WMS
&version=1.1.1
&layers=topp%3Astates
&styles=
&srs=EPSG%3A4326
&format=image%2Fpng
&bbox=-145.151041%2C21.73192%2C-57.154894%2C58.961059
&width=780
&height=330
&query_layers=topp%3Astates
&info_format=text%2Fhtml
&feature_count=50
&x=353
&y=145
&exceptions=application%2Fvnd.ogc.se_xml

But I want to send a request a polygonal area. Not a single x,y location. How can I do this? Is possible?


texts: Using less parameters for Python script tool?

texts:
I am looking to create a user-input script for finding surrounding features of a given point. However, the parameters for it are very sensitive and I would like for it not to have to be such a process to create results.

This is the script that I have so far:

import arcpy
import sys

inFc = sys.argv[1] #input feature class (for some reason must be shapefile because the feature classes in a geodatabase don't have FIDs
nearFc = sys.argv[2] #near feature class for near table analysis (feature  class)
outTable = sys.argv[3] #output table for near table analysis, also used as input for field name change so that the "IN_FID" field changes to "FID" to match the input and so the selection can be performed (table parameter)
distance = sys.argv[4] #search radius for near table analysis
outLyr = sys.argv[5] #X/Y event layer so that the near table analysis can be presented on the map
projection = sys.argv[6] #projection for the x/y coordinates
selId = sys.argv[7] #Query expression obtained from input feature class so the user can choose one feature from the input to highlight all the surrounding near features from

arcpy.GenerateNearTable_analysis(inFc, nearFc, outTable, distance, "LOCATION", "NO_ANGLE", "ALL", "", "GEODESIC")
arcpy.AlterField_management(outTable, "IN_FID", "FID")
arcpy.MakeXYEventLayer_management(outTable, "NEAR_X", "NEAR_Y", outLyr, projection)
arcpy.SelectLayerByAttribute_management(outLyr, "NEW_SELECTION", selID)

As you can see, there are too many parameters and it would not be easy for a user to perform this quickly for a demonstration. Taking into consideration my comments, are there any possible ways that this could be streamlined and adjusted?


texts:
I've written following code to make a summary list of no of features having same "ID" values so that I can later remove duplicate "ID" values from them if any.

import arcpy
from arcpy import env
mxd = arcpy.mapping.MapDocument("CURRENT")
layers = arcpy.mapping.ListLayers(mxd)
for layer in layers:
    lyr = layer
    desc = arcpy.Describe(layer)
    path1 = str(desc.path) + r"\outstat.dbf"

    env.workspace = desc.path
    arcpy.Statistics_analysis(lyr, path1, ["SUM"], "ID")

However this script gives error 000800 whenever I run it.
Exact error message is:

Runtime error  Traceback (most recent call last):   File "<string>", line 11, in <module>   File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\analysis.py", line 1160, in Statistics     raise e ExecuteError: ERROR 000800: The value is not a member of SUM | MEAN | MIN | MAX | RANGE | STD | COUNT | FIRST | LAST.

How can I resolve the error I am making?

There is only one layer in the mxd.


texts: FCC Database in PostGIS, display in ArcMap

texts:
I used a modified version of the scripts at: https://bitbucket.org/dirkcgrunwald/fcc-uls to pull the FCC databases down and put them into a PostGIS environment.

I created a view so I didn't have to build a query every time I wanted something. No doubt it'll evolve over time as I'm pretty green to this.

SELECT row_number() OVER () AS oid,
fr.unique_system_identifier,
en.call_sign,
en.entity_name,
fr.frequency_assigned,
fr.frequency_upper_band,
(lo.lat_degrees + lo.lat_minutes / 60)::numeric + lo.lat_seconds / 3600::numeric AS latitude,
((lo.long_degrees + lo.long_minutes / 60)::numeric + lo.long_seconds / 3600::numeric) * (-1)::numeric AS longitude,
fr.power_output,
fr.power_erp
   FROM pubacc_fr fr
     JOIN pubacc_lo lo ON fr.unique_system_identifier = lo.unique_system_identifier
     JOIN pubacc_en en ON fr.unique_system_identifier = en.unique_system_identifier
  WHERE fr.class_station_code = 'WDX'::bpchar;

I added the row_number bit in an attempt to have a column with a unique ID so that I could open the table in ArcMap. So far my attempts have been futile. While a select * from weather_radar; does indeed show the OID column, I don't see it when I browse the view in phpPgAdmin and ArcMap doesn't see it either when I connect to it.

Thoughts?

My drive to open the attribute table is to figure out why lat\long are limited to five digits when there's at least a dozen showing in phpPgAdmin when I browse the view. This has the unfortunate side-effect of lumping all points in a given metro into a small area nowhere near the correct point for any of them.


texts: Selecting features to project to different coordinate system?

texts:
I've got data from all over Europe in one feature class projected in Web Mercator Auxiliary Sphere. Due to the fact, that length and area values are not usable I'm now working on a python script (ArcGIS 10.0) that projects the features to the corresponding UTM zones and writes the calculated values back to the orginal fc.

First the script assigns the corresponding UTM zone to each feature via spatial join with an UTM zone fc. Within the next step I'd like to select all features within one UTM zone and project them (e.g. all features with value 32 will be projected to UTM ETRS89 zone 32 and so on). My problem is, that project_management does not have a "where_clause" so what will be the best way to make this selection? I've got lot's of data in one feature class and many feature classes, so speed does matter...

Does anyone know a better way to get precise length and area values for data in Web Mercator projection?


texts:
I am trying to call a WMS using OpenLayers. I have wrote the code but somehow I think I am missing something here.

My code is:

<!doctype html>

    <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.2.1/css/ol.css" type="text/css">

    <title>OpenLayers WMS</title>

My Map


texts: How to count near intersecting polylines and merge them?

texts:
I have a dataset with several polylines created by GPS points of car trips. I am trying to create a simple heatmap where roads with more traffic (where polylines intersect) get darker.

I already tried overlapping each polyline and then applying a sligh opacity to each one and that works. The problem is that with huge amounts of data, the request to the server for the data is too big.

My idea is, in the server, process each trip and merge each overlapping polyline into a single one with and associated weight. Is there any way of doing this?


texts: Displaying symbols properly in ArcGIS Desktop?

texts:
I opened one point shapefile with ArcGIS 10.2 and I want to change the symbol of it.

When I click on the symbol, the ESRI symbols in the Symbol Selector don't show properly, but like "!" and "&".

What can I do to show them properly?


texts: Numbering Watersheds using ArcGIS Desktop?

texts: I am working with a watershed map, and I was wondering if there is a function or a way that I can get ArcMap to number the watersheds. For instance, if we have watershed X numbered as 0, I want the watersheds that surround/or share a boundary with watershed X to be numbered 1, and the watersheds that surround/ or share a boundary with the number 1 watersheds to be numbered 2 and so on. Is there a way this can be done by building a model in ModelBuilder? I would like to save it so when I change a watershed layer and run the model it will do the same in a new watershed area.


texts: Cannot Set TOC View/Tab to Display with VBA?

texts:
I have some code from 9.3 that would check the CurrentContentsView.ContextItem to ensure that the TOC was in Display (or Drawing Order) view.  If not it would set it to Display via:

If Not pDoc.CurrentContentsView.ContextItem = 1 Then
    Set pDoc.CurrentContentsView.ContextItem = 1
    pDoc.ContentsView(1).Refresh Nothing
End If
pDoc.ContentsView(1).Refresh Nothing

Now that the TOC is a DockableWindow (which I can access) and there are several IContentsViews, this no longer works.  It's surprising that this is impossible to find in the forums and the VBA SDK as as brief and cryptic as could be.

How do I access the ContextItem (which seems stuck at NULL) to set it? Also, is the fact that "currentContentsView" automatically reverts to a lower-case "c" when the object apparently shouldhave an upper-case "C" significant?

Any help would be MUCH appreciated. For some reason, my code runs 3x faster in Display View and I want to automatically set it if the user forgets (or doesn't know).

Oh...and yes, I AM aware of VBA's demise and the existence of Python. Noted.

Thanks in advance.


texts: Calculate the footprint along the track between two distant squares

texts:
I am trying to achieve this kind of operation programatically, using preferably PostGIS (or else GDAL)

Is there any well-known operation doing this?

My use case is the following:

I have two squares representing the footprint of the pictures that a satellite shot at the start and the end of a semi-orbit.
I would like to be able to reconstruct the polygon of the area that has been scanned along the path between these two snapshots.
The width of the swath is of course constant.
The start and end squares can not be aligned, and therefore the resulting polygon shall be calculated on the envelope (and not only a dilation of the path). In that sense the above picture doesnt represent exactly what I need (ie: the resultant shape shall be wider in its center, since resulting of the 'dragging' of the first shape to reach the second).


texts: How do we import multipolygon relations from OpenStreetMap properly with osm2pgsql?

texts:
We are using osm2pgsql to import the OpenStreetMap's planet.osm file to GeoServer. Then we try to get features from the polygon layer using GDAL/OGR's WFS driver in C#.

We've tried the import with and without the --multi-geometry option.

The multi-polygon I'm using for testing is http://www.openstreetmap.org/relation/69663. I try to fetch the relevant features with
layer.SetAttributeFilter("(osm_id=69663)OR(osm_id=28414267)OR(osm_id=30407439)");
The above ID's are the relation ID of the multi-polygon, the way ID of the outer way and the way ID of the inner way, respectively. This query returns only a single feature with a single geometry with 22 nodes and no nested geometries. This feature matches the outer way in ID, number of nodes and its tags. This means that not only is it missing the inner way; it is also missing the tags from the relations, which is what identifies it as a building.

What happened with the inner way? Was it discarded because it doesn't have any tags? What happened with the relation? Most importantly, how do I get the result I seek?

EDIT: iant wanted to know what was in the geometry column. The geometry column is named "way" and for the osm_id 28414267 it contains:
0103000020E61000000100000016000000F988F3CBBB012A40ECEB0AC73DCD4B40B24CBF44BC012A400F046D173FCD4B40015B142CC4012A40399787F13FCD4B40F8CDB348C9012A4050357A3540CD4B407F750B13FC012A400E6EC6C441CD4B405BE4E83B09022A40263ACB2C42CD4B40DB3E9AA03C022A4060F01FC143CD4B406DECC84741022A4001E8E6D143CD4B403BAC70CB47022A409C22765C43CD4B40F107A8154B022A408A5B60EA42CD4B40B46DBD4B4E022A405BF1683E42CD4B4048BD4CB967022A404015376E31CD4B40A5CD829A6B022A409BDC39DE2ECD4B40D0807A336A022A401E79D6242ECD4B409EFB613A63022A4095DD825B2DCD4B406BC1E677F5012A40FA3D569522CD4B4063EFC517ED012A40FA3D569522CD4B40E92D7997E6012A407D1C72D822CD4B4057804AF0E1012A40F9C1548E24CD4B40DD8431D8C3012A40BF93BB7438CD4B40A977A7F1C1012A406B3304B639CD4B40F988F3CBBB012A40ECEB0AC73DCD4B40
We converted it to text and it contains exactly what I got with OGR: 22 nodes, no nested geometries.

EDIT 2: The above hexadecimal data as text:
POLYGON((13.0033859 55.6034478,13.0033895 55.6034879,13.0034498 55.6035139,13.0034888 55.603522,13.0038763 55.6035696,13.0039767 55.603582,13.0043688 55.6036302,13.0044043 55.6036322,13.004454 55.6036182,13.0044791 55.6036046,13.0045036 55.6035841,13.0046976 55.603071,13.0047272 55.6029928,13.0047165 55.6029707,13.0046633 55.6029467,13.0038259 55.6026179,13.003762 55.6026179,13.0037124 55.6026259,13.0036769 55.6026781,13.0034473 55.6032854,13.0034328 55.6033237,13.0033859 55.6034478))
I'm too lazy to check all of them, but the eleventh node (13.0045036 55.6035841) matched the first node on OpenStreetMap.


texts: Moving all shapefiles (SHP) in MXD to new File Geodatabase (*. gdb) using ArcGIS Desktop?

texts:
What I need to do is to move all the shapefiles in an mxd into a newly created gdb, checking prj, and update mxd reference.

I have ArcGIS 10 SP4, license ArcView.

I studied a bit of Python and get familiar with basic, then I copied some written code and pasted all together, attached is the code.

With some help from ESRI forum I ended up with the following script where I get an error in the arpy.SpatialReference line.

# Import system modules
import arcpy
from arcpy import env
import os
import traceback
# Load required toolboxes...
arcpy.AddToolbox("C:\\Program Files (x86)\\ArcGIS\\Desktop10.0\\ArcToolbox\\Toolboxes\\Data Management Tools.tbx")

# Set environment settings
env.workspace = "C:\\Users\\a391665\\Documents\\AAAAPERSONALE\\rugareto\\GIS\\test"
print("env.workspace completed successfully")

#check GDB exist
if os.path.exists("c:\\park\\test4.gdb"):
    arcpy.Delete_management("c:\\park\\test4.gdb")

# Execute CreateFileGDB
arcpy.CreateFileGDB_management("C:\\park", "test4")
print("CreateFileGDB completed successfully")

# Set local variables
outWorkspace = "c:\\park\\test4.gdb"
print("outWorkspace completed successfully")

# Iterate
def fcs_in_workspace(workspace):
    arcpy.env.workspace = workspace
    for fc in arcpy.ListFeatureClasses():
            dsc = arcpy.Describe(fc)
            if dsc.spatialReference.name == "WGS_1984_UTM_Zone_32N":
                print os.path.join(workspace, fc, "prj: WGS_1984_UTM_Zone_32N")
                arcpy.FeatureClassToGeodatabase_conversion(fc, "c:\\park\\test4.gdb")
                print os.path.join(workspace, fc, "importaed in GDB")
            elif dsc.spatialReference.name == "Unknown":
                print os.path.join(workspace, fc, "UNKNOWN")
                #HERE I'D LIKE TO IMPORT THE SHP IN GDB ADDING "_UNKN" TO THE FEATURE NAME
            else:
                print os.path.join(workspace, fc, "DIFFERENT PRJ!")
                # Determine the new output feature class path and name
                outfc = os.path.join(outWorkspace, fc)
                # Set output coordinate system
                print "PRJ CREATED!"
                # Create a spatial reference object using a projection file
                outCS = arcpy.SpatialReference('WGS_1984_UTM_Zone_32N')
                print "outCS CREATED!"
                arcpy.Project_management(fc, outfc, outCS)
                print "PROJECTED!"
    for ws in arcpy.ListWorkspaces():
        fcs_in_workspace(os.path.join(workspace, ws))

fcs_in_workspace("C:\\Users\\a391665\\Documents\\AAAAPERSONALE\\rugareto\\GIS\\test")


texts:
I have a set of points defined by their longitude and latitude, and I would like to check if they lie within a "polygon" on sphere. By polygon on sphere, I mean area in between 4 lines going between 4 points that I have. As I understand it, Polygon will give me straight lines that are cast onto a sphere, and I want arcs that are the shortest possible connections between the points. So how may I do that?
edit: by do that, i mean this:

SELECT * FROM ST_Distance(ST_GeographyFromText('POLYGON((70 -40,70 -39,71 -39,71 -40,70 -40))'), ST_GeographyFromText('POINT(70.48 -39)'), False)

should return 0, because the POINT lies on the boundary arc.
I'm using postgresql 9.4 with postgis 2.1.8


texts: Using the results of geoprocess function in another geoprocess function

texts:
I have a list to hold the iterated results of running arcpy.MakeFeatureLayer and 'arcpy.CopyFeatures` geoprocesses.

selectList = []
fip = '13245'
for fc in featureclasses:
     select = arcpy.MakeFeatureLayer_management(fc, "temp_lyr","FIPS_CODE = '{}'".format(fip))
     arcpy.env.workspace = "C:/Workspace/Sandbox/selects.gdb"
     selectFeatures = arcpy.CopyFeatures_management(select, "select_{}".format(fip))
     selectList.append(selectFeatures)

When I print(selectList) it gives me a weird list of results:

[<Result'C:\\Workspace\\Sandbox\\MapChangeProject\\selects.gdb\\select_13245'>, <Result 'C:\\Workspace\\Sandbox\\MapChangeProject\\selects.gdb\\select_06037'>]

What I want to do is merge all of these results together into one feature like this:

arcpy.Merge_management(selectList, "merge")

What am I doing wrong here and how do I get my results into a format that will work for my merge?


texts:
I am currently working on a computer with windows 7 on it and i need to instal the Pg-Trgm module for postgresql 9-5. I have to admit that i am kind of lost here since i don't have any idea how to do this on a windows workspace.

I need it for matching text based columns between them and for some slight manipulations.


texts: How do I change the default color of the canvas in QGis 3.4?

texts: The default color is white. How may I switch it to another color ? To be clear, I'm talking about the background, not the menus, nor panels or toolbars, etc.


texts: No FeatureTypeList in WFS published by QGIS Server!

texts:
The WFS URL I'm using to get service from QGIS Server has no data about the layers and features type. It should be noted that the WMS URL is ok. Is there any problem with the URL I'm using? Should I add any further tag?
The WFS URL:
http://localhost/cgi-bin/qgis_mapserv.fcgi.exe?SERVICE=WFS&VERSION=1.0.0&REQUEST=GetCapabilities&MAP=E:/test.qgz


texts:
Currently we are working in a geological project..in which all the Contacts, Linear & some other geological features have to be captured.

The thing is most of the features have to be Captured as "Line" and then it has to be converted to "Polygon". So that we don't have to repeat the digitization for features like Contact....

1) Is it not at all possible in Qgis to Snap a line with the same line ? like that of in AutoCAD, Mapinfo..??  (Image 1)

2) I've captured a polygon as two separate lines...and when i use "Line to Polygon" (Enclose in the case of Mapinfo)....it actually creates two polygons instead of one single entity. Is it a bug ? (Image 2)

I am trying to convince some of colleagues (and the management as well) to switch over to QGIS..atleast in few systems...So i am trying to figure out the available functionalities in QGIS as that of MapInfo. So any sort of feedback will be of great help.

Update: Polygonize command does the trick of converting the two separate lines into one single polygon feature...(Image 2).


texts: are these two proj4 utm projection strings equivalent?

texts:
I have a list of proj4 calls. Some of them have one format and some have another, I've been told they are equivalent. Is that so? If so, which one should I use?

for example, for ETRS89 / UTM zone 20N:

+proj=utm +zone=20 +ellps=GRS80 +towgs84=0000000 +units=m +no_defs
+proj=utm +zone=20 +ellps=WGS84 +datum=WGS84 +units=m +no_defs


texts: Speed up intersection, union in PostGIS?

texts:
I have three tables: water, wetlands and parcels. I want to update parcel table with area that is covered by water or wetland (they can overlap).

So I want to

Merge water with wetland
Dissolve and multipart -> singlepart
Intersect this with parcels to calculate area:

Number of parcels is ~400 000, water 200 and wetlands 10.

The Query runs for 35 min and update ~2200 rows with correct results. This seems very slow so im probably doing something inoptimal. What can I do to speed it up?

update parcels

set water_or_wetland_area = sub.area_m2

from parcels,
     (select parcels.parcel_id, sum(round(st_area(st_intersection(waterwetland.geom,
     parcel.wkb_geometry))::numeric,2)) as area_m2
     from parcel,
         (select (st_dump(st_union(st_union(water.wkb_geometry,
         wetland.wkb_geometry)))).geom
         from water, wetland
         ) as waterwetland
     where st_intersects(waterwetland.geom, parcels.wkb_geometry)
     group by parcel.parcel_id) as sub
where parcels.parcel_id = sub.parcel_id


texts: Toggling data source in Chloropeth map with Leaflet

texts:
I'm starting with the Choropleth map example from Leaflet. (http://leafletjs.com/examples/choropleth/)

I have a very large geoJSON file (all US counties), containing data points over time.  To match the example cited, let's say the data series is population density in 1980, 1990, etc.  I would like to map all of this data.  An easy solution would be to create a separate layer for each decade, but the file would take too long to load.  I'm hoping for a solution with one geoJSON file.

var statesData = {"type":"FeatureCollection","features":[
{"type":"Feature","id":"01","properties": {"name":"Alabama","density2000":94.65, "density1990":90.00, "density1980":80.000},"geometry":    {"type":"Polygon","coordinates":[[[-87.359296,35.00118]

I plan to have radio buttons at the top of the map where the user can select which decade he wishes to see on the map, using an onclick event function.

The code to update the information text box on the upper right hand corner of the map is here:

info.onAdd = function (map) {
    this._div = L.DomUtil.create('div', 'info');
    this.update();
    return this._div;
};

info.update = function (props) {
    this._div.innerHTML = '
US Population Density
' +  (props ?
        '' + props.name + '
' + props.density + ' people / mi2'
        : 'Hover over a state');
};

info.addTo(map);

This function is called when the user hovers over the state here:

function highlightFeature(e) {
    var layer = e.target;

    layer.setStyle({
        weight: 5,
        color: '#666',
        dashArray: '',
        fillOpacity: 0.7
    });

    if (!L.Browser.ie && !L.Browser.opera && !L.Browser.edge) {
        layer.bringToFront();
    }

    info.update(layer.feature.properties);
}

And when the map is initially created...

geojson = L.geoJson(statesData, {
    style: style,
    onEachFeature: onEachFeature
}).addTo(map);

Basically, I want the props.density to change to props.density1980 or props.density1990 based on what the user selects.  This sounds simple conceptually - assign a variable based what the user selects, and use "if" statements accordingly.  I'm having trouble coding this though - I think I'm getting mixed up with the onClick procedures.


texts: Piping custom processing algorithm output into another in QGIS?

texts:
I am using QGIS 3.0.

Following up on Calling a QgsProcessingAlgorithm from a Processing python script in QGIS 3.0?, I am writing a Processing python script based on the template provided here.

It takes a few additional parameters in addition to the 'INPUT' and 'OUTPUT' from the template. I am using these parameters to create some lines and everything works out fine.

Now I need to call a Processing Algorithm, namely qgis:clip to clip the lines I just created according to the Polygon entered as 'INPUT'.

Is there a way to pipe my output into qgis:clip ? Or do I need to create a temporary layer that I will then feed into qgis:clip and let it's 'OUTPUT' become mine ?

Either way, I'm not sure where to start.

I've tried defining a QgsFeatureSink or a QgsProcessingFeatureSink but the first is an abstract thingy and the second just won't be recognized (I can't even import it from qgis.core)

I've seen Using processing using other output (result from another algorithm) as input and Import error for qgis.core when running OSGeo4w shell script but I couldn't figure out a working solution to my problem so far.