texts: Cannot read postgis RASTER into ArcGIS using the interoperability extension, either directly by Query Layer

texts:
This post Accessing PostGIS raster data using Query Layer in ArcGIS Desktop? was created in Jul-12. Since then, I have not found any other related-post published.

Using QGIS with DBManager ON, my postGIS raster tables can be displayed nicely so far. That's cool and it proved that my raster tables are OK and not corrupted. However, I need to use ArcMap 10.1 and read those rasters, without having to create views or multy-polygons in the database as the above post suggested.

I have tried two ways:

Setting a Database connection in ArcCatalog directly to the Postgress/PostGIS database. I can display vector features, but only see the attribute table of the raster, so not displaying the raster as a whole. I guess this has not been implemented so far by ESRI, has it?.
By using the Interoperability extension and its postGIS file type, I can see my geometries within my database, anything new, so I am able to add those raster tables into my ArcCatalog tree. Interesting, at least this time, the icon of the element instead of being a table icon seems to be a point feature icon, so something different is going on. But when I try to display these point geometries (a kind of cells?), it says "Error reading Interoperability Data. The operation failed because the database table does not have a numeric number".

If I were able to add this numeric number into my raster table at the database - still I don't know how to do it- would I be able to finally visualize the raster?, or would I see just only the table again?. Would version 10.2 of ArcGIS improve this a little?. Is there any other way round to solve this within ArcGIS?

In the other hand, am I having these issues because I have not done something at the database level which I was supposed to do with my raster files? Someone has succeeded at adding postgis rasters into ArcGIS?


texts: Set new layer to SVG marker in plugin QGIS 2.18

texts:
Im trying to create a new layer and change it to SVG marker and set blue-marker.svg as the marker but its not working.

Here is what i have:

        iface.addVectorLayer( "Point", "marker_layer", "memory")
        mlayer = self.iface.activeLayer()
        symbol = QgsSymbolV2.defaultSymbol(mlayer.geometryType())
        svgStyle = {}
        svgStyle['name'] = '/usr/share/qgis/svg/symbol/blue-marker.svg'
        svgStyle['outline'] = '#000000'
        svgStyle['outline-width'] = '0'
        svgStyle['size'] = '6'
        symbol_layer = QgsSvgMarkerSymbolLayerV2.create(svgStyle)
        symbol.changeSymbolLayer(0, symbol_layer)

Im not getting any errors but the new layer is not set as SVG neither is the default marker to blue-marker.svg

EDIT:
here is the solution:

        layer = self.iface.activeLayer()
        svgStyle = {}
        svgStyle['name'] = '/usr/share/qgis/svg/symbol/blue-marker.svg'
        svgStyle['outline'] = '#000000'
        svgStyle['outline-width'] = '0'
        svgStyle['size'] = '12'
        symbol_layer = QgsSvgMarkerSymbolLayerV2.create(svgStyle)
        layer.rendererV2().symbols()[0].changeSymbolLayer(0, symbol_layer)


texts: uDig window to change between historical dates?

texts:
I'm new to uDig, but I'm experienced Java developer. I have points with temperature measurements that change overtime in SQL database. I would like to give user ability to change between historical data.

Is there window in uDig which allows this?
If not: Can you point me where to start to write a plugin allowing this?


texts: Improving appearance of vector tiles in Mapbox at various zoom levels

texts:
I am wondering how I can better render large datasets in Mapbox for viewing purposes?

Currently I use an ogr2ogr script to convert large ESRI shapefiles to MBTiles.

ogr2ogr -f MBTILES -dsco MINZOOM=10 -dsco MAXZOOM=22
C:\users\myself\path_to_output\target.mbtiles
C:\users\myself\path_to_shp\source.shp

From this I then upload these MBtiles to Mapbox. However one issue I am noticing is the improperly rendered features at higher zoom levels.

This is an example of the Mapbox preview I upload the MBTile at zoom levels <16. This is an example of a shapefile that contains 200k point features within it. At zoom levels <16 there is a strange single point pattern than occurs rather than the more natural clustering that you see at the top.

However at zoom level 17, this is no longer an issue and is appears to be as normal.

Ideally there would be a way to display the datasets with a feature simplification at higher zoom levels that look more natural as opposed to the "dot grid" pattern in the first image.

One option I looked at was tippecanoe however is doesn't seem to be built natively for windows and the recommended install is via MacOS.

Would it be worth the hassle to install tippecanoe or is there an alternative available for Windows users for better rendering of large datasets in Mapbox for viewing?


texts:
I'm new to CityEngine and struggling to put together code for a zoning rule that requires buildings to have a graduated height based on what's called the context transition edge. Max height is a ratio in relation to distance from the property line.

I understand basic CGA but need help thinking this through.


texts: Showing PostgreSQL table comments in ArcGIS Pro?

texts:
Comments are a useful way to annotate database objects such as fields and tables in PostgreSQL. They show up as comments when a PostgreSQL layer is loaded in QGIS, as evidenced in the information pane:

QGIS also shows field-level comments:

My question is, does ArcGIS Pro have the same capability? When I load the layer there, the 'table metadata' shows up as blank:

I might not be looking in the correct place, but I'm not seeing where else it might be hiding.
Alternatively, is there another method to include PostgreSQL comments when loading layers into ArcGIS Pro?


texts:
When I try to run NVIZ from QGIS I get an error which states 'Can't find usable init.tcl in the following directories" and then goes on to list a number of directories. (I did have a screenshot here but couldn't post it as a newbie.)

However when I try to run NVIZ directly from GRASS it runs properly.  I have searched for the init file which the QGIS error report says cannot be found and it is actually in "C:\Program Files (x86)\Quantum GIS Lisboa\lib\tcl8.5".  According to the QGIS error report this is not one of the directories actually searched.  I do have a workaround running NVIZ from GRASS but that is clearly not as satisfactory as running from within QGIS.  I have noted the reply to an earlier question on this topic which recommends installing tcl/tk via OSGEO4W. However I do not believe that these components are actually missing since I've found the init file.

Can anyone suggest a solution, please? I am a new user of QGIS and have downloaded 1.8.0 Lisboa, running under Windows 7.


texts: How to do simple self join and touches function in postgis?

texts:
I am having trouble with this query and neighbor calculation in PostGIS. I want to find out which are the neighbors of Sweden in the table COUNTRIES which has column 'the_geom' for geometry and 'cntry_name' for countries. This is what I tried to do but with little success. Can someone explain to me what am I missing?

SELECT s1.cntry_name
FROM countries AS s1, countries AS s2
WHERE s1.cntry_name = ‘Sweden’ AND
touches(s1.the_geom, s2.the_geom);

After editing the query this got me two results, Finland and Norway!

SELECT s1.cntry_name as sweden_neighbors
FROM countries AS s1, countries AS s2
WHERE touches(s1.the_geom, s2.the_geom) = 'True'
AND s2.cntry_name = 'Sweden';


texts: Working with python add-in while editing is on

texts:
I created a python add-in which makes analyses like spatial join, update cursor and make feature layer. But when editing is open for a layer. I get this error.

Line 289, in onClick
    for rec in recs:

RuntimeError: The operation is not supported within an edit session.

targfields3 = ['SHAPE', 'yapi_id', 'mahalle_id', 'YOL_ID', 'Yol_id_butun', 'YOL_ADI']
        with arcpy.da.UpdateCursor('KAPI selection', targfields3) as recs:
            for rec in recs:
                keyval = rec[0]
                if joindict3.has_key(keyval):
                    rec[1] = joindict3[keyval][0]
                    rec[2] = joindict3[keyval][1]
                    rec[3] = joindict3[keyval][2]
                    rec[4] = joindict3[keyval][3]
                    rec[5] = joindict3[keyval][4]
                else:
                    rec[1] = 0
                    rec[2] = 0
                    rec[3] = 0
                    rec[4] = 0
                    rec[5] = 0
                recs.updateRow(rec)
        del rec, recs

This is the operation mentioned in error.

While making this analyses manually, when edit is open on arcmap, I get no error.

I have to work with this python add-in while edit is on. So, solving this error is crucial for me.


texts: Using QGIS Queries With ArcGIS Server FeatureServer

texts:
I am trying to load USGS Transportation Vector Data from the US National Map data set.

REST Link: https://carto.nationalmap.gov/arcgis/rest/services/transportation/MapServer

If I use the feature server I am able to load layers

However, clicking "Build Query" throws the following error

Protocol "" is unknown

Presuming that this feature is functional for other users, I am wondering what I am doing wrong to be able to build queries for this Feature Server?

There are two things that tell me this should be possible and I am missing some URL construction.

I can do explicit queries as JSON vector layers
Checking the box "Only request features overlapping the current view" works as intended. The mechanism by which it works is through a query...

I have tried adding ?query&f=pjson and a few other combinations of urls that worked as explicit layers.

This is a very large dataset so must be queried, there are many of them so I would prefer the convenience of the FeatureServer layer dialog. Without query there is no way for me to load some of the more detailed layers.

The feature server advertises as queryable

"capabilities": "Data,Map,Query"


texts:
We are just looking into moving from MapServer to GeoServer for some of out WMS needs, one of which is to serve OS MasterMap to our users.
Before I have a go at creating the SLD styles for MasterMap does anyone have any they have created that they are willing to share?


texts: How to view/edit M-values of PointZ shapes ArcMap 9.3.1

texts:
I realize this is a very basic question, but however I can't find any basic Q&A on something so simple. I added three new points to a shapebatch in ArcMap. By default when I added them, it made the M-Value the equivalent of a minimum double value (-1.79769313486232E+308). I know this because I've written software to read ShapeFiles. For testing purposes I need to be able to change the M-values of these points to see how it handles certain values.

I've tried opening the 'Edit Sketch Properties' windows. It's grayed out and doesn't even contain an M column. All tutorials/guides I've come across pertain to PolyLineZ.

Is the M-value computed on something within the shapefile?


texts: Display a .KAP file and associate coordinates to it

texts:
I am not an image format type expert and I am working on a project that needs to display a map and associate it to GPS coordinates. The map files are downloaded from a website that gives me two zipped folder for each map. The first one contains a .KAP and a .BSB file. The second one contains a .tif and a .db file.
What should I do in order to display the map and associate a GPS position to it?


texts: Is the QGIS 3.10 "fix geometries" batch-processing interface missing a parameter?

texts:
My goal is to write five features classes from one FGDB to five tables in one GeoPackage using QGIS 3.10 (Dutch version). I'm trying to do this using the QGIS function fix geometries (Dutch "geometrieën repareren).
If I 'fix' one feature class at a time, everything is fine. But when I try to run the fix geometries function as a batch process through the QGIS user interface (see attached image), QGIS tries to write five GeoPackages with the same name. The name of the output Geopackages can be set as a parameter (GeoPackageInputlayer), but my goal is to write all the output of this batch in to the same GeoPackage. This option isn't supported in the user interface.
Is there a way to write the five tables into one GeoPackage using the GQIS batch process interface of the FixGeometries function?
This is the input for the batch processing of the function 'fix geometries':

And this is the error (file already exists):


texts:
I'm composing a map of Big Island, Hawaii in QGIS 3.14.0 and would like to display major roads and cities + towns. Therefore I downloaded necessary information via QuickOSM. However, as many towns lie close together, the label density is quite high:

Comparing the map with the representation from openstreetmap.org, I highly prefer the representation/ density from osm:

How can I achieve the same density in QGIS? The downloaded data (which should be the same for both representations) does not contain any population information or other fields which could be used for filtering out unwanted towns. I've found following workaround which reduces label density but does not lead to the same representation.
Edit:
I'm aware of the place tag which discriminates between different place types and that they can be styled differently. Asking about displaying towns I was actually referring to place = town, sorry for the ambiguity. I marked the same group of towns in both maps. With provided fields from QuickOSM I'm not able to discriminate further between those towns to achieve the same display density. As the underlying data is the same I think it should be somehow reproducible.
Edit 2:
I'm accepting JoshC's answer that maps shown on the OSM website indeed don't have any additional data which would further discriminate between towns. The different representation stems from the used OSM styling parameters. Reproducing with GIS could work but seems not straightforward, as

tile styling is a little different from the way other GIS programs and
services handle labelling.


texts: Add an image to one feature within ArcGIS Online

texts: Using the web map pop up config, is there a way to add an image for one asset (polygon)? It appears the only option is to add an image for all the assets in that layer.


texts: Does Geoserver update layer changes automatically from PostgreSQL/PostGIS?

texts:
If someone makes changes to tables in PostGIS (that have already been published in Geoserver), will Geoserver notice the changes automatically in its stores or do the layers need to be re-published?

Would there be a way to schedule a bulk layer update every 24 hours?


texts: QGIS Map Composer, right hand pane has gone!

texts:
in the QGIS map composer, I accidentally closed the right hand pane, the pane that includes the command history, item properties, etc. Now I can't get them back. There are no icons to click to get them back, nothing. If I close QGIS and reopen, the pane is still gone.

What do I do to get the pane back?


texts: Simplifying Polygons and MultiPolygons for Visualization

texts:
I completed a computation with Canadian Census data and discovered, when trying to visualize the results with Bokeh 2.0.2, that the Polygon and Multipolygon paths were very high resolution with a point every 30 meters or so.

My first strategy was to use GeoPandas' built-in Shapely.Simplify to simplify the geometry but the test results were less than satisfactory, and the full rendering of the dataset still wouldn't complete in my notebook regardless.

gdf_simplified["geometry"] = gdf.simplify(tolerance=50000, preserve_topology=True)

I was wondering what some next steps could be to try to visualize this geometry? Should I try matlibplot? Are there are other topology preserving algorithms to try?


texts: Is there an open standard for raster information?

texts:
For geographic vector data (Features) there is an OGC standard that defines geometries and their different annotations such as WKT and WKB.

For raster image data however, I was not able to find comparable standards. PostGIS uses something like WKB represenations in their PostGIS Raster, but even for the PostGIS raster I did not find a document that would describe the Raster Type other than the PostGIS manuals which are more user oriented and do not explain the standard in the same way as a standardization document would.

Does anyone know if there is such a standard for Raster?

Links to detailed definitions of the PostGIS raster would also be useful here if there is no standard.


texts: Is it possible to log in to Esri Geoportal Server programatically?

texts:
Esri Geoportal Server has the option to secure datasets, such that it's necessary to log in to view them:

My questions:

Is it possible to interact with this login functionality programatically?

I'd like to provide the login form on another web page, then pass the parameters to Esri Geoportal in the background, without displaying the above user interface.

Can I detect whether the user is logged in (or not) via the Geoportal's REST API? The documentation doesn't mention user management.

(There is a sample Geoportal Server instance at http://gptogc.esri.com/geoportal/catalog/main/home.page which may be useful for testing.)


texts: Difference between QGIS developer version and QGIS stable version for Python programming?

texts:
I am using python and PYQGIS and I have installed QGIS via OSGEOW. I chose the complete advanced install with all packages, libraries and GIS programs.

There is a 2.99 developer's version and the latest stable version (2.18.7).

What additions does the developer's version feature, and what differences does it make for python programming?


texts: Adding coordinate system to layout in QGIS composer?

texts:
I would like to know if there is a way to add a dynamic coordinate system field into a QGIS composer layout.

I can use an expression to add many items such as document name, file path, etc but I can;t seem to figure out how to add the coordinate system?


texts:
I have 45 tif tiles each of 20-50MB and a single shp file with about 4000 polygons. I want to run zonal histogram in Arcmap 10 to get the number of pixels with each value for each polygon.
If I convert a single polygon from the shp to a raster and use it as the zone input, I can get an output that is correct. However, for the whole shp file either as polygons or as rasters, the tool runs and reports that it completes (very fast, a few seconds), but produces no table.

I think it's because there are too many zones, maybe?

My plan had been to run a batch of zonal histograms, one for each of the 45 tifs.

Am I right that I have too many zones?

And what options for workarounds do I have?


texts: QGIS: Looking for an 'Update' function similar to ArcGIS

texts:
I am doing an analysis of a road network.  The attribute table holds information on the road name, type (e.g. A road, motorway etc).  I would like to group all the minor roads as a single field.  I have tried the dissolve function in QGIS, selecting all the minor roads and dissolving the selection.  This process has been running for over 3 hours and is only on 7%!  Forums suggest other people have found this tool incredibly slow.

I ran the Line dissolve from the processing toolbox which took about a minute top and saved the local roads as a new layer.  I'm hoping I can overlay this somehow onto the original layer as a work around to the dissolve function not working properly.  I believe the function I would need in ArcGIS would be 'Update' but I can't find a similar tool in QGIS.

Any suggestions?


texts: How to make a WFS Layer read only in Geoserver?

texts:
So I imported some postcode data into a Postgis table and published the layer using Geoserver.

As it is postcode data, there is absolutely no need for users to edit, so I would like to make the layer read only.

How to do this?


texts: How to update Geodatabase point file table from Excel data using FME

texts: I have excel spreadsheet with data that I want to load into a GeoDatabase point  table. The table already exists and I want to update the data not append. There is a unique ID to link the Excel data to the point file. When I use Joiner in FME or load data tool in ArcGIS it appends the data rather than updating.


texts: Using Single house subaddress locator in ArcGIS Desktop?

texts:
I am running into a problem with the single house subaddress locator in 10.4.

I have two different 'groups' of addresses with units/suites/etc.

The first 'group' has addresses with the same house number and street name/type, and different units:

1234 Smith Ave Unit A

1234 Smith Ave Unit B

The second 'group' is like the first, but there is a base address record without a unit:

5678 Jones Rd

5678 Jones Rd Unit A

5678 Jones Rd Unit B

Each address has a distinct point feature.

My fields:

HouseNum, Name, Type, SubType, SubNum

I configure my locator with SubType = SUBADDTYPE and SubNum = SUBADDRUNIT. All of the other address components are mapped like normal.

Here's what happens:
Addresses from the first group locate perfectly to the point with the correct unit.

Addresses from the second group ALL locate to the base address point that doesn't have a unit.

So it seems like the locator stops looking for unit information (even if it's mapped) once it successfully locates the standard street address.

Any ideas on how to get it to look further, all the way into the unit fields?

I also have the 'additional field' mapped to a field that contains the full, unparsed address.


texts: Popup from InfoTemplate only showing attributes in Firefox, but not in Chrome or Internet Explorer

texts:
I have a map service I am testing on internally with some popups on a dynamic map service layer. They work perfectly in Firefox (version 45.0.1), but not in Chrome (version 50.0.2661.102) or Internet Explorer (version 11.0.96 on my computer). I am using the ArcGIS JS API 3.16.

These are the only meta tags I have in my HTML file:

<meta http-equiv="Content-Type" content="text/html;charset=utf-8"/>
<meta http-equiv="X-UA-Compatible" />
<meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>


texts: spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelWithin; not working

texts:
I have a task of getting the total area of all the polygon features within a boundary.
So I got the envelope of the boundary and did a spatial query to get the features inside that element

but my problem is I am able to get the elements that overlaps the envelope but not exactly inside it..
my code goes like this

IMxDocument mxDoc = (IMxDocument)GetMxDocumentFromArcMap(m_application);
        IEnvelope gridEnvelope = GetGridExtents("Z29");
        IFeatureClass featureClass = GetFeatureClassFromShapefileOnDisk("D:\\Projects\\Working\\BMRDA\\data", "Hoskote_Satellite_Polygon");

        // Create the spatial filter and set its spatial constraints.
        ISpatialFilter spatialFilter = new SpatialFilterClass();
        spatialFilter.Geometry = gridEnvelope.Envelope;
        spatialFilter.GeometryField = featureClass.ShapeFieldName;
        spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelWithin;
        spatialFilter.set_OutputSpatialReference(featureClass.ShapeFieldName, mxDoc.FocusMap.SpatialReference);

        // Set the attribute constraints and subfields.
        // We want to exclude ramps, highways and interstates.
        spatialFilter.WhereClause = "LEVEL3 = 'AGRICULTURE LAND'";

        // Execute the query.
        IFeatureCursor featureCursor = featureClass.Search(spatialFilter, true);
        ESRI.ArcGIS.Geodatabase.IFeature feature;  // Automatically initialized to null. Used to test existence of a feature in the featureCursor
        double dbArea = 0.0;
        int intCount = 0;
        while ((feature = featureCursor.NextFeature()) != null)
        {
            IArea area = (IArea)feature.Shape;
            dbArea = dbArea + area.Area;
            intCount++;
        }

the spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelWithin; does not work it returns me null I want the polygon to be clipped and get the area of it
Am i missing something here?
Please guide..


texts: Using max function to group max value by unique ID in QGIS

texts:
I have been using the maximum function to group max value by unique "id";
maximum("Ni_ppm_1", group_by:="Hole_id")

The function pulls the maximum "Ni_ppm_1" value per unique "Hole-id".
Below is a snip of the dataset

Is it possible to get the function to pull out the "from" and "To" associated with the maximum Ni value?
The "From" / "To" is the depth of the sample.


texts: Get all dates of image collection when it was created/loaded on sentinel 2 dataset

texts:
I am using "COPERNICUS/S2" image collection for satellite imagery on Google Earth Engine using nodejs, I wanted to highlight the dates in calendar for particular month on which the image was created/loaded in dataset. Like the one working on EO Browser(see attached screenshot)

How can I get those dates with imageCollection?


texts: Digitized vectors don't match while zooming

texts:
I'm trying to digitize in QGIS with Google Maps as the background. Digitizing with constant zoom works fine. But if the zoom is changed, the digitized vectors don't match the Google Maps background image. Thanks for your help!

Regards
Hans


texts: How to know what value use to convert meter in degree using Google Maps info

texts:
I am trying to create a mathematical formula to convert meters to decimal degrees.

Reading this article Decimal Degrees, I thought this generic formula:

x = (Value_in_meter * 0.00001)/1.1132

But I know it is not 100% correct, I should use the other values according where is my point.

I am using Google Maps Api, so, how to discover if my point is at 23N/S, 45N/S or 67N/S?

Complementing

I did this function:

public static double convertMeterToDegrees(double meter, double latitude){

    double quotient;
    double degree = Math.floor(latitude);
    double modDegree = Math.abs(degree);

    if (modDegree == 0){
        quotient = 1.1132;
    } else if (modDegree <= 23){
        quotient = 1.0247;
    } else if (modDegree <= 45){
        quotient = 0.7871;
    } else {
        quotient = 0.43496;
    }

    return (meter * 0.00001)/quotient;
}

Is it correct?


texts: Is there any open source implementation of GeoXACML based access control mechanism for OGC Web Services?

texts: I am looking for an open source implementation of GeoXACML PEP, PDP to provide access control to my OGC web services. I have already set-up SAML based authentication. I found a non-opensource implementation from Secure Dimensions but I want to know if there is any opensource implementation.


texts: Write variable content into specific cell of datasheet (Mapbasic)

texts: I have a value stored in a variable. How can I write this value to a table (for example into Row2 Column4) using Mapbasic?


texts: Export (to any format) leaves out pieces of DEM

texts: I have a map comprised of three different sources of mosaiced DEM (two Lidar and one 1/3 arcsec float nationalmap DEM) which show up entirely intact in ArcMap (10.2.2), when I export to PDF at a DPI over 300 or even as a lossless TIF file, parts of the DEM drop out and are white. The vector graphics (line or area) are intact over the white space, just the DEM is not there. It does not affect one source of data exclusively, and the missing sections are not consistent from export to export. Also, in my elevation legend, one of the ticks for the elevation numbers comes through looking fatter than the other ticks... All of these are circled on the attached screen shot: ArcMap on left, PDF and zoomed in screenshot of PDF on right. Does anyone have any idea how to combat this? 200 dpi is not acceptable!


texts: Deleting row value from Multiple point feature in different sub-folder using ArcPy

texts:
I have different month folders, each folder contain a shapefile called, For January folder Merged001.shp, For FEBRUARY Merged002.shp so on,

That shapefile having one field name 'Solar', and i am trying to delete specific  value from Solar field or we can say whichever value < 15 its will be deleted from Solar field (which is < 15.)

I have made a script to loop through all the month folder and read it the shapefile's Solar field and delete it < 15 value.

Below is my script and error massage

    import arcpy, os, calendar
    from arcpy import env
    from arcpy.sa import *
    arcpy.CheckOutExtension("Spatial")
    topWorkspace = r'D:\Delete_Test'
    env.workspace = topWorkspace
    env.overwriteOutput = True

    months = {calendar.month_name[i].upper(): str(i).zfill(3) for i in range(1, 13)}
    for folderPath in arcpy.ListWorkspaces():
        baseName = os.path.basename(folderPath).upper()
        if baseName in months:
            monthNumber = months[baseName]
        arcpy.env.workspace = folderPath

        shapeFiles = (arcpy.ListFeatureClasses('Merged_*.shp')[0])
        print shapeFiles

    with arcpy.UpdateCursor(shapeFiles, ["Solar"]) as cursor:
 for row in cursor:
    if row [0] < 17:
        cursor.deleteRow(row)
del row

error massage :

Traceback (most recent call last):
  File "<module2>", line 19, in <module>
AttributeError: __exit__


texts: Overpass query Paris metro stations with line numbers?

texts:
I understand from the wiki article of Paris Metro, that there are 16 lines and 302 stations. Could you help me run an overpass query that can give me the data? i.e I can see which station belongs to which line?

I have seen this. This API is built by Overpass DE but do not know whether this is pulling data direct from OSM.

Edit: unless I am not able to spot it, neither of these two does the job in R using osmdata package.

paris = opq(bbox = "paris") %>%
  add_osm_feature(key = "station", value = "subway") %>%
  osmdata_sf()

paris2 = opq(bbox = "paris") %>%
  add_osm_feature(key = "type:RATP", value = "metro") %>%
  osmdata_sf()


texts: How to make polygons clickable in Leaflet.js / folium

texts:
I'm trying to create a clickable leaflet.js map using folium and Python.

After a previous query - answered here Geopandas/folium map not displaying - I've added the polygons.

But I'm getting a little bit lost in the documentation on what to add to the output javascript to create a pop up for each polygon using the dataframe.

Can anyone give me a steer in the right direction, please?


texts: Change a raster from longitude display (-180, 180) to (0, 360)

texts:
Using R software, I want to modify the attached raster in order to visualize it centered on the Pacific Ocean. Indeed, I need to show the entire Oceania: in my opinion I have to shift from (-180,180) to (0,360)...

Coud you provide me the appropriate code?


texts: How can I show points on the same coordinate as e.g. a bigger symbol in QGIS?

texts:
I use QGIS 1.8 for mapping archaeological findings. Some of the findings within one layer are located on the same site and therefore have the same coordinates. Therefore you just see one point instead of the actually overlapping points. There is the tool - in German its called "Punktverdrängung" (maybe point displacement in english?), with which you can circle overlapping around a central point, but that does not really fit my needs.

Is there  the possibility of showing more than one point on the same coordinates as a bigger symbol or anything like that?

Thanks for your help!


texts: How can I convert my scale bar from milimetres/ degree into metres or kilometres using QGis?

texts: I am making a map in print composer. I am trying to add the scale bar and it's appearing as a box with a zero digit inside it. What should I do to make it display properly?


texts: Purpose of filterEl.addEventListener('keyup', function(e)

texts:
I'm learning the mapbox gl js examples and have trouble understanding a function.
The function can be found here https://www.mapbox.com/mapbox-gl-js/example/filter-features-within-map-view/ it's in fact the "filterEl.addEventListener('keyup', function(e)" in the example.

I already removed the complete function and noticed no difference in behaviour.

This is the code:

 filterEl.addEventListener('keyup', function(e) {
        var value = normalize(e.target.value);

        // Filter visible features that don't match the input value.
        var filtered = airports.filter(function(feature) {
            var name = normalize(feature.properties.name);
            var code = normalize(feature.properties.abbrev);
            return name.indexOf(value) > -1 || code.indexOf(value) > -1;
        });

        // Populate the sidebar with filtered results
        renderListings(filtered);

        // Set the filter to populate features into the layer.
        map.setFilter('airport', ['in', 'abbrev'].concat(filtered.map(function(feature) {
            return feature.properties.abbrev;
        })));
    });


texts: Select features intersecting a linestring in QGIS

texts: I wrote a python plugin that grabs all (vector, linestring) features between two points on the map.  At the moment, I select all features intersecting the rectangle around these two points, and then use shapely to filter out features not actually intersecting the line between the two points.  This is inefficient, and I'm wondering if there is a way to perform the correct operation directly in QGIS (regardless of the data-source backing the layer)?


texts:
GOAL: Take a raster and reproject it from one co-ordinate set to another and also change the pixel size while I am at it by using GDAL.

WHY: I am trying to use gdal to speed up multiple rasters having to be reprojected which is taking a long time in ArcPy.

Current progress: I can get this to happen nicely by using ArcPy and the following segment of code:

arcpy.ProjectRaster_management(in.tif, out.tif, "PROJCS['WGS_1984_UTM_zone_48N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['false_easting',500000.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',105.0],PARAMETER['scale_factor',0.9996],PARAMETER['latitude_of_origin',0.0],UNIT['Meter',1.0]]", "NEAREST", "100", "", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")

The GDAL code: I am using the segment of code below to hopefully try to obtain the same result.

call(['C:\Program Files (x86)\GDAL\gdalwarp.exe',"-s_srs", "EPSG:4326","-t_srs","EPSG:32648","-tr","100","100","-r","near","-srcnodata","0.0","-dstnodata","0.0","-of","GTiFF",in.tif, out.tif])

Actual result from gdal: A dataset that is reprojected to the correct co-ord system and pixel resized to the desired size.

The problem: The raster dataset when opened and worked with keeps showing that the min and max values it contains is the nodata value of +-3.40282e+038

The key point: Does anyone have any idea why this is not working correctly in gdal?

While trying different different options I found that by running

call(['C:\Program Files (x86)\GDAL\gdalwarp.exe',"-s_srs", "EPSG:4326","-t_srs","EPSG:32648",,"-srcnodata","0.0","-dstnodata","0.0","-of","GTiFF",in.tif, out.tif])

I found that the resulting cell size from the projection is ~1000x1000 so I don't know why it won't let me scale to a smaller cell size?

Extra checks:

I have checked the spatial reference website to ensure that the EPSG co-ord I am using as the correct ones.
And as far as I can tell from the gdalwarp info page the code should work nice.


texts: QGIS 3x can't run the custom crs for Ordnance Survey maps

texts:
Recently I solved the problem regarding the mismatch between British National Grid CRS and WGS 84 in the query below:
QGS 3x Ordnance Survey map doesn't overlay the Google Maps properly
I downloaded the .gsb file and copied it into the relevant folder like shown in the instructions.
However now I see, that this .gsb file (with custom CRS) is valid for one layer only. I can't see it for other Ordnance Survey Maps...

You can see the layer above with "Unknown CRS" (located in unknown Coordinates System section) which comes from the copied .gsb file unlike another layer below, which doesn't have this layer in the search box at all.


texts: Can OpenStreetMap be used as a WMS service in MapServer?

texts:
I want to produce a png image that uses OpenStreetMap data and has a vector overlaid on the base data. I am able to produce the vector overlay from within Mapserver, but I have difficulty understanding how to include an external data source like OpenStreetMap.

I would prefer to use an external service over loading the OpenStreetMap data directly as I wish to use it only as a basemap.

Is it possible to do this from within MapServer?


texts: Why doesn't ArcGIS rasterize every polygon

texts:
I am trying to rasterize elements from a vector map but somehow ArcGIS only lets me rasterize 15 objects. I added a new field so that every object would get the same value. The original layer had 87 objects in it. Why do I now only get 15 in the rasterized layer?

The operation

The target area

The original table

The table after rasterizing

I've done this many times before and I can't figure out what is going wrong..


texts: Creating a temporay layer out of a selection in QGIS

texts:
I want to create a tempory layer out of an existing layer and only get those features from the existing layer, which are selected. I know that ArcGIS has this functionality that you can create a layer from selected features.

However, I could not find any working solution for this in QGIS. I am aware of this question which describes, how to save a selection into a new file which is not my goal here.

Checking the context menu for "save as..." did not yield any results. In the plugin section with the search word "selection" did not give any meaningful results either.

While trying around I came up with a solution with pyqgis but I would prefer one which does not involve any coding because I need to show this to a colleague.

*pyqgis way:

layer = iface.activeLayer() # layer which has a selection
feat = layer.selectedFeatures()
fields = layer.fields()

selection = QgsVectorLayer('Point', 'temp', 'memory')

dp = selection.dataProvider()
dp.addAttributes(fields)
dp.addFeatures(feat)

selection.commitChanges()
selection.updateExtents()

QgsMapLayerRegistry.instance().addMapLayer(selection)


texts: How to forbid zoom map under scroll page?

texts:
I wrote code at ArcGis JavaScript.
There is a map on web page at full display with a fixed height  (400px)
When i scroll-down page and cursor misses in zone map then browser works for zooming map, not scroll page further. How to forbid?


texts: How to batch Check geometry in ArcGIS for a list of datasets from a CSV?

texts:
We have around 50+ gdb's(in different folders) which has 450 datasets in total. Among the 450, we have a list of 100 feature dataset's for which we have to do QAQC. So is it possible to batch Check geometry process ? From ArcGIS's help page I've come across this http://pro.arcgis.com/en/pro-app/tool-reference/data-management/check-geometry.htm This code works great for a single gdb to check geometry of all the datasets within that gdb. I would like to run Check geometry process only for particular datasets within particular gdb's.

Eg: D:\QAQC\f_c.csv This will have the list of datasets for which Check geometry process has to be run. CSV will have file path like

E:\Chk\A\ddd.gdb\Feature_dataset_1

E:\Chk\B\eee.gdb\Feature_dat_1

I know the file_path should not be represented like this. Just to throw some idea on what I want, I've mentioned the path like this

Edit:

This code worked like a charm

import csv
import arcpy

csvpath = r"path\to\csv"
with open(csvpath, "r") as csvfile:
    fcs = [r[0] for r in csv.reader(csvfile)]
arcpy.CheckGeometry_management(fcs, "outtable")


texts:
I have imported SRTM DEM into Oracle 11g using FME.The block size is set to 128X128.I am able to view the data in QGIS,Tatuk and play with it using GDAl command line.Unable to view the same in ArcMap. Never.Never faced any issue with vector though.Why Raster loading is not as simple as vector ? Apparently , all documentation that I have come across has evasive responses providing no clue at all on what to do in simple steps.

ArcCatalog shows me The raster and RDT tables.I get the error "The specified raster column does not exist", when try to drag and drop.I have verified that my DBTUNE parameters are set SDO_GEORASTER.

Any clues?


texts:
i have some vector layers in geoserver with style (SLD)

When i add these layers though OL3 i want to use same SLD.

right now i am giving style to layer as below

style: new ol.style.Style({
                    stroke: new ol.style.Stroke({
                      color: '#f00',
                      width: 1
                    }),
                    fill: new ol.style.Fill({
                      color: 'rgba(255,0,0,0.1)'
                    })
                  })

How can i get the same SLD which is present for the layer in geoserver and use in OL3?


texts: How to programatically replace ArcGIS Online image with new imagery in mxd's?

texts:
I have many mxds that previously used ArcGIS Online Bing imagery. I want to replace those with newly purchased imagery. I know I need to use the arcpy.mapping.MapDocument module to do this though not sure how. However when I try to identify the source for the online imagery it is listed as 'source unknown'. Script snippet below shows how I am trying to identify sources in each mxd which I am writing to a .csv

for map in mxd_list:
mxd = arcpy.mapping.MapDocument(map)
layers = arcpy.mapping.ListLayers(mxd, '*')
for i in layers:
    try:
        ds = i.dataSource
    except:
        ds = 'source unknown'
    writer.writerow([map, i.name,  ds ])


texts: Different icon of the nearest marker to current location in Leaflet

texts:
I'm working with Leaflet and JavaScript for the first time, so I feel a bit lost here. I have an app that displays user's location and a couple of points of interest from a GeoJSON on a map. The points have a popup bound to them that displays information when clicked.

Now I would like the marker for the closest POI to the user be different than the other markers. I have also managed to find the closest point using L.GeometryUtil.

//GEOJSON LAYER
//----------
function onEachFeature(feature, layer){
    if (feature.properties && feature.properties.name) {
        layer.bindPopup(feature.properties.name);
    }
}

var geojsonLayer = L.geoJSON(pivovary, {
    pointToLayer: function (feature, latlng) {
        return L.marker(latlng, {icon: beerMarker})},
    onEachFeature: onEachFeature
}).addTo(map);

//USER LOCATION
//----------
var closestpt = null
map.locate({setView: true, maxZoom: 16});
function onLocationFound(e){
    L.marker(e.latlng, {icon: youAreHereMarker}).addTo(map),
    closestpt = L.GeometryUtil.closestLayer(map, geojsonLayer.getLayers(), e.latlng).latlng
}

map.on('locationfound', onLocationFound)

I tried to add a new marker at the closestpt coordinates in onLocationFound, but then I am just placing a new marker on top of one that is already there and I cannot use the popup that is bound to my geojsonLayer anyway.

Is there a way to render the markers after I find the closest point?


texts: How to return absolute values using ArcMap field calculator?

texts: I have a field of number ranging from -50 to 50. Is there a way using the ArcMap field calculator to populate a field with all of these standardised (positive numbers) so I have a range of 0 to 50 (where something that was previously -1 is now 1).


texts: How can a leaflet SVG overlay be placed outside of the leaflet-overlay-pane?

texts:
I have SVG on top of GeoJSON in my map. How can a leaflet SVG overlay be placed outside of the leaflet-overlay-pane while keeping GeoJSON in place in the leaflet-overlay-pane? (For example, how could I put just the SVG into the leaflet-shadow-pane and then turn pointer-events off on the leaflet-shadow-pane... thus making it possible to interact with the GeoJSON?)

My apprehension is this:
Leaflet map layers are stored in panes:
http://leafletjs.com/reference.html#map-panes

However, it is not clear to me just exactly how these panes are ordered and what each of them contains.

according the Leaflet's css (http://cdn.leafletjs.com/leaflet-0.7.2/leaflet.css), the order of these panes makes more sense:

.leaflet-map-pane    { z-index: <"?">; } //Pane that contains all other map panes. (presumably all the below map panes)

.leaflet-tile-pane    { z-index: 2; } //Pane for tile layers.
.leaflet-objects-pane { z-index: 3; } //Pane that contains all the panes except tile pane. (presumably all the below panes)

.leaflet-overlay-pane { z-index: 4; } //Pane for overlays like polylines and polygons.

.leaflet-overlay-pane svg { z-index: <"?">; }
.leaflet-overlay-pane <"?"> { z-index: <"?">; }

.leaflet-shadow-pane  { z-index: 5; } //Pane for overlay shadows (e.g. marker shadows).
.leaflet-marker-pane  { z-index: 6; } //Pane for marker icons.
.leaflet-popup-pane   { z-index: 7; } //Pane for popups.

I am specifically interested in .leaflet-overlay-pane svg  and any subsequent panes that reside inside .leaflet-overlay-pane (i.e. GeoJSON pane).

Thanks


texts:
Is it possible to draw polylines with one point (x,y) and given specific angles e.g. for 45 degrees, and given length?
In other words: without specifying the end point.


texts:
I have a number of images, (these could be vectorised) that are label classifications. There are only 2 classes, 1 or 255. Each image is from a different date. For each date I also have data that can be in csv or point form.

I wanted to try and use logistic regression to get a future, binary prediction using the images, as the dependent variable and the point data associated with each pixel as the independent variables.

I've done this with one pixel, I need to try and extend the model to a wider area and I'm having trouble visualising how to do this. I had thought of having it all as numerical data and looping through the process for each pixel but this would take a very long time.

I am not too familiar with GIS so I could be missing what's possible, but I had also thought of converting the raster to points and populating the attribute tables for each point by joining tables by geolocation. I could convert this values in to rasters as well using spatial interpolation.

I believe there's also some methods in R or python. Or maybe there's some data fusion method I'm missing. I'm a bit stumped on the matter so I'm looking for ideas on how this might be done.

Maybe someone has experience with a similar problem?

Basically I want to input images and text data/point data in to the model and have my output as a predicted image. Maybe this is impossible in GIS.


texts:
I want to perform a simple calculation with the Field Calculator using this piece of script (note that all the elif statements are attempts for another syntax and non of the worked):

def change(OldFor_Com):
    if (OldFor_Com) == "OHWH":
        return shape.area@Hectares / 40
    elif (OldFor_Com) == "OHWH":
        return !shape.area@Hectares / 30!
    elif (OldFor_Com) == "OHWH":
        return (!shape.area@Hectares / 30!)
    elif (OldFor_Com) == "OHWH":
        return !shape.area@Hectares! / 30
    elif (OldFor_Com) == "OHWH":
        return (!shape.area@Hectares! / 30)
    if (OldFor_Com == "OTHH"):
        return !shape.area@hectares! / 40

If the field is populated with OTHH I want another field to be populated with the area of the field, divided by 40, or 30, or... and so on.
I modified my syntax a couple of times (as you can see) but it still says Invalid Syntax.

Is a calculation possible with these if statement or do I need a workaround?

Note that the shape.area@Hectares is an expression to calculate the area of the field, and not another field in my attributetable!

EDIT: 7/15/2014 around 11pm

I modified the script. I added a field to the table and calculate the area of the shape first (with !shape.area@Hectares!) and then I try to run this

def change(OldFor_Com):
    if OldFor_Com == "OTHH":
        return Area_Ha / 40
    if OldFor_Com == "OHWH":
        return Area_Ha / 30
    if OldFor_Com == "OSFH":
        return Area_Ha / 375
    if OldFor_Com == "PINE":
        return Area_Ha / 10
    if OldFor_Com == "OOFH":
        return Area_Ha / 375
    if OldFor_Com == "None":
        return Area_Ha
    else:
        return 0

It bounces back with the following:

ERROR 000539: Error running expression: change(u"OTHH", 7)
Traceback (most recent call last):
File "<expression>", line 1, in <module>
TypeError: change() takes exactly 1 argument (2 given)

Failed to execute (Calculate Field (10)).

Is it not possible to do any calculations in the if statement?

EDIT: 7/15/14 12:30pm

This how the attributestable looks like before I add the field SizeScore

After I add the field SizeScore (FLOAT).

And the following screen comes up if I use this code:

def change(OldFor_Com, Area_Ha):
    if OldFor_Com == "OTHH":
        x = Area_Ha / 40
        return x
    if OldFor_Com == "OHWH":
        x =  Area_Ha / 30
        return x
    if OldFor_Com == "OOFH":
        x = Area_Ha / 375
        return x
    if OldFor_Com == "OSFH":
        x = Area_Ha / 375
        return x
    if OldFor_Com == "PINE":
        x =  Area_Ha / 10
        return x
    if OldFor_Com == "None":
        return Area_Ha
    else:
        return 0

It only populates the rows where it does not need to calculate.
The only statement that is working is that one.

if OldFor_Com == "None":
     return Area_Ha

But for the others are no error messages coming back!

SOLVED, SEE ACCEPTED ANSWERE
Thanks a lot to everybody who went through this with me!


texts: Normalizing raster using raster calculator in ArcGIS 10.3.1?

texts:
I have been trying to 'range standardise' or normalise my raster values from 0-100 and I know how to do this mathematically:

  ("raster"-min("raster"))/(max("raster")-min("raster"))*100

and I have input this equation into the raster calculator using both the min() and max() and the .minimum and .maximum syntax.

Why are these not working for me? ArcGIS 10.3.1


texts:
I am trying to make a selection of my "layer" and create a "new_layer" of the selection.
Next I want to append data to the memory layer. But nothing get appended with this script, nor get I any messages.
After the data is appended I will write it to a permanent layer.
Can I append data to a memory layer like I am trying to do or what seem to be wrong?

try:
# Make memory layer from subset of feature class
    new_layer = arcpy.MakeFeatureLayer_management("layer","sub_layer","Contract = AB23" )

except:
   print(arcpy.GetMessages())

DataToAppend = "c:/data.gdb/MyLayer"
appendData = []
appendData.append(DataToAppend)

# Append data to layer from input
  arcpy.Append_management(
    inputs=appendData,
    target=new_layer,
    schema_type="NO_TEST",
)

# save layer permanently
arcpy.CopyFeatures_management("sub_layer", out_location)

-----------------------EDIT---------------------------

This script work. But its feels a bit bloated and I am need to delete the copied feature class if appendData fails somehow. I will make this dynamic and add it as a script to a toolbox

# Execute CopyFeatures to make a new copy of the feature class
arcpy.CopyFeatures_management(inFeatures, outFeatures)

# Create temporary layer based on the copied FC
arcpy.MakeFeatureLayer_management(outFeatures, "tempLayer")
print("Created feature layer")

# Select Subset of temp layer
arcpy.SelectLayerByAttribute_management("tempLayer", "NEW_SELECTION", "Contract IN ('AB123')")

if int(arcpy.GetCount_management("tempLayer")[0]) > 0:
    print("Selected items")
    arcpy.DeleteFeatures_management("tempLayer")

new_data = c:/data.gdb/MyLayer
appendData = []
appendData.append(new_data)

# Append data to layer from input
arcpy.Append_management(
    inputs=appendData,
    target=outFeatures,
    schema_type="NO_TEST",
)


texts: ArcHydroTools FillSinks Error HRESULT E_FAIL in ArcMap

texts:
I'm trying to use ArcHydroTools in ArcMap 10.2.2 and this error always appears in the finish of the process:

File "<stdin>", line 1, in <module>
  File "C:\Program Files (x86)\ESRI\WaterUtils\ArcHydro\bin\ArcHydroTools.py", line 6805, in FillSinks
    raise e
arcgisscripting.ExecuteError: System.Runtime.InteropServices.COMException (0x80004005): Error HRESULT E_FAIL has been returned from a call to a COM component.
   at ESRI.ArcGIS.Geoprocessing.GeoProcessorClass.Execute(String Name, IVariantArray ipValues, ITrackCancel pTrackCancel)
   at ESRI.APWR.ApHydro.ArcHydroOp.FillSinks(ApLayers& aplayers, Boolean useIsSink, Boolean fillAll, Double fillLimit, ITrackCancel trackCancel, IGPMessages messages, String& exMessage) in C:\Builds\HydroSolutions\10.2_ArcHydro\Sources\ArcHydroTools\src\ApHydro\ArcHydroOp.vb:line 64688
Failed to execute (FillSinks).

My code is:

import arcpy
import ArcHydroTools

from arcpy import env
ws = env.workspace = r"C:/python_results"
input_name = "C:/python_results/project.tif"
output_name = "C:/python_results/Fil"
ArcHydroTools.FillSinks(input_name, output_name)

I tried this options (which I found in another answer here) but nothing changed:

moved all of my data under C:\data\
created new file geodatabase (named testing.gdb) under that folder
increased the max number of unique values to render
set out the targets (C:\data for raster and C:\data\testing.gdb for vector files) using ApUtilites
set out the processing extent to the extent of my DEM (Geoprocessings>Environments>Processing extent ) )
saved my project as testing.mxd under C:\data
set testing.gdb to be default workspace. (this finally solved the last errors I was getting)
also named my data frame simply to "data" since ArcHydro -tools tend to add data frame name to output-file paths. I'm not sure if step this is needed.

I know that this error normally is about the communication between machines, but maybe is about computer memory, I don't know...


texts: Table output from QGIS graphical modeler

texts:
I have developed my own model, which processes input vector layer. The final output is vector layer that contains many attribute fields, and of course, also the geometry, but I need a simpler output - geometry and most of those fields are unnecessary.

With this model, how can I get a simple table output with selected attribute fields of the output layer?


texts: Difficulties with placement properties in maplex labeling

texts:
I'm working in ArcGIS 10.2.2 and am having two difficulties with labelling using the Maplex Engine.

The first problem is that I have two labels that I can't seem to put any additional space in between. (The circled area on the right side of the attached image.) Both labels are labelling features in the same class in the same layer. The feature weight is 0, the label buffer is maxed out at 50, with a hard constraint, yet they are placed right next to each other. What I would prefer is to get the Starr Street label moved to the left so that it is the fairly open area to the left of Newhall Street, but the label won't seem to move there.

My other issue is that I have one north-south street where the label runs top to bottom instead of bottom to top, as all the other streets are. I assume this is because this is the only street on a curve, and the maplex engine determined that this is the best way to place it using a centered curve. Is there any way to force the label to run the other way instead?


texts: Reorder Shapefile of lines for consecutive ordering (Assemble Lines into LineString)

texts:
I have a Shapefile containing segments of Lines/LineStrings:

Unfortunately the order of the DBF attributed table is random.

How to reorder the segments so entry #1 in the Shapefile will be "Segment 1", and the consecutive entry will be a segment which:

Share same coordinate point geometry (Blue dot example: "Your end is my beginning").

OR

First coordinate point geometry is the closest to the last coordinate of previous segment (Orange & green dots example).

The logic here is that sometimes a physical road LineString explicitly separated into Lines (segments), and sometimes mistakenly a small gap can exist.
(actually I think that most of the times condition 2 will met, because if 2 segments share same Geometry point, They also probably will be connected as single LineString)

Note- The question is just how to find the ordering and reorder (The next step will be to form single LineString either by connecting via the share coordinate or by connecting Last->First coordinates of consecutive lines).

EDIT:

With Join Multiple Lines plugin It's possible to merge segments "As Expected" in the above example, including removing duplicate node (blue).
The problem that it merge all other attributes into single one: For example column "ID" so will be only id=1.
Then applying Explode Lines tool, splits the merged line back to segments at nodes with the correct order, however because of what mentioned at point 2 - the data inside the attributes is lost.


texts: How to unproject a UTM ESRI Shapefile over WGS84 extent boundaries

texts:
Background:

I am working on a collaborative project in Python that uses GDAL/OGR to work with rasters and vectors. There is one function that attempts to project one ESRI Shapefile's coordinate system into another. This is done by creating a coordinate transformation from one spatial reference to another, then reading through the geometries and attempting a Transform, which is then written to an output vector of the desired projection. Read: reprojecting a shapefile.

Question:

I've ran across one application where it is desired to "unproject" from UTM Zone 3N (Alaskan Region) to WGS84 (lat/long). Doing this causes distortion and bad extents for the output, where the returned shapefile's polygon stretches from -180 to 180 degrees. I believe the issue is due to the projected shapefile crossing the extents of the -180 degree boundary. To confirm, I broke the projected shapefile into 3 parts, west of -180 degrees, middle at -180 degrees, and east of -180 degrees. When unprojecting west / east to WGS84, the returned shapefiles were as expected, however, the west portion lived at positive 180 and the east portion lived at -180 degrees. The middle portion was returned distorted. To visualize this better I provide 2 photos below, the first being the projected data as explained, and the second being the broken results as just mentioned:

I would like to know the proper way of handling this unprojection in Python using GDAL/OGR. I'm a developer by trade and don't have a thorough understanding of GIS projections. May be similar in ideas to this post: How to stop gdalwarp creating world-spanning outputs near the dateline?


texts: Opacity control for KmlLayer in Google Maps V3 API?

texts: Using a KmlLayer object within the Google Maps API, is it possible to control the opacity of a KML's image overlay features ?


texts: Moving/Reprojecting raster to a different central meridian?

texts:
I have a raster set to 0-360E projected using WGS:1984.

When I modify the Data Frame Properties GCSWGS1984 to a custom meridian of 180, the raster does not move. The polygon of the continent moved but not the raster.

Also, I want the raster to display at 160. Because at 0-360, UK (of course) is cut lengthwise. So I want it at 160.

I am using ArcMap 10.3.


texts: How can I get the UK Land Cover 2007 25m raster data set from Digimap to work in QGIS?

texts:
I have downloaded the UK Land Cover 2007 25m raster data set from Digimap (which my university has a subscription to). I have uploaded the data into ArcMap where it works fine, with all 23 classes and colours shown. However, when I upload the data into QGIS (which all my other data works on) it says the raster is not a supported file type. It appears in black, grey and white and I cannot find a way to get it to appear with the 23 colours and classes it should have. When I downloaded the data set it came with; AUX file, XML document, OVR file, DBF file and LYR file(which has raster in the name).

Is there a way to get it to work in QGIS? I am open to any suggestions. Could I convert it to a supported file type in Arc and bring it into QGIS?


texts: Optimize ArcServer Service Drawing Speed

texts:
I am managing an ArcGIS service that is used by Geocortex Essentials as a web map app for some users. The feedback I have got so far is that it doesn't refresh fast enough. I'm looking for some answers to these questions:

Do labels with a mask on them draw slower?
What point symbol draws faster: Character or Simple?
Does having a label only draw in its polygon result in slower labelling?


texts: Error compiling PostGIS 3.1 from source : "undefined reference to `GEOSIntersectionPrec'"

texts:
I am just trying to build the new PostGIS 3.1.0 (with GEOS 3.9.0). But I consistently get an error about "undefined reference to `GEOSIntersectionPrec'".
The commands I used was (which was successful with PostGIS 3.0.x, but not with 3.1.0{,alpha3,beta1,rc1}):
wget https://download.osgeo.org/postgis/source/postgis-3.1.0.tar.gz
tar xf postgis-3.1.0.tar.gz
pushd postgis-3.1.0
./configure --with-geosconfig=/opt/geos/bin/geos-config --with-sfcgal=/opt/sfcgal/bin/sfcgal-config
make

The error says:
...
---- Making all in loader
make[1]: Entering directory '/home/user/src/postgis-3.1.0/loader'
/bin/bash ../libtool --mode=link \
  gcc -I ../liblwgeom  -I/opt/geos/include   -I/usr/include/libxml2 -I/opt/sfcgal/include -I/usr/include/json-c  -DNDEBUG  -std=gnu99 -g -O2 -fno-math-errno -fno-signed-zeros  -fPIC -DPIC    shpopen.o dbfopen.o getopt.o shpcommon.o safileio.o shp2pgsql-core.o shp2pgsql-cli.o ../liblwgeom/liblwgeom.la -o shp2pgsql -L/usr/lib/x86_64-linux-gnu -Wl,-Bsymbolic-functions -Wl,-z,relro -Wl,-z,now -L/usr/lib/llvm-9/lib  -L/usr/lib/x86_64-linux-gnu/mit-krb5 -Wl,--as-needed -static ../liblwgeom/liblwgeom.la  -lc
libtool: link: gcc -I ../liblwgeom -I/opt/geos/include -I/usr/include/libxml2 -I/opt/sfcgal/include -I/usr/include/json-c -DNDEBUG -std=gnu99 -g -O2 -fno-math-errno -fno-signed-zeros -fPIC -DPIC shpopen.o dbfopen.o getopt.o shpcommon.o safileio.o shp2pgsql-core.o shp2pgsql-cli.o -o shp2pgsql -Wl,-Bsymbolic-functions -Wl,-z -Wl,relro -Wl,-z -Wl,now -Wl,--as-needed  -L/usr/lib/x86_64-linux-gnu -L/usr/lib/llvm-9/lib -L/usr/lib/x86_64-linux-gnu/mit-krb5 ../liblwgeom/.libs/liblwgeom.a -lm -L/opt/geos/lib -lgeos_c /usr/lib/x86_64-linux-gnu/libproj.so -ljson-c -L/opt/sfcgal/lib -lSFCGAL -lc
/usr/bin/ld: ../liblwgeom/.libs/liblwgeom.a(lwgeom_geos.o): in function `lwgeom_intersection_prec':
/home/user/src/postgis-3.1.0/liblwgeom/lwgeom_geos.c:710: undefined reference to `GEOSIntersectionPrec'
/usr/bin/ld: ../liblwgeom/.libs/liblwgeom.a(lwgeom_geos.o): in function `lwgeom_unaryunion_prec':
/home/user/src/postgis-3.1.0/liblwgeom/lwgeom_geos.c:788: undefined reference to `GEOSUnaryUnionPrec'
/usr/bin/ld: ../liblwgeom/.libs/liblwgeom.a(lwgeom_geos.o): in function `lwgeom_difference_prec':
/home/user/src/postgis-3.1.0/liblwgeom/lwgeom_geos.c:840: undefined reference to `GEOSDifferencePrec'
/usr/bin/ld: ../liblwgeom/.libs/liblwgeom.a(lwgeom_geos.o): in function `lwgeom_symdifference_prec':
/home/user/src/postgis-3.1.0/liblwgeom/lwgeom_geos.c:891: undefined reference to `GEOSSymDifferencePrec'
/usr/bin/ld: ../liblwgeom/.libs/liblwgeom.a(lwgeom_geos.o): in function `lwgeom_union_prec':
/home/user/src/postgis-3.1.0/liblwgeom/lwgeom_geos.c:1044: undefined reference to `GEOSUnionPrec'
collect2: error: ld returned 1 exit status
make[1]: *** [Makefile:108: shp2pgsql] Error 1
make[1]: Leaving directory '/home/user/src/postgis-3.1.0/loader'
make: *** [GNUmakefile:22: all] Error 1

All the components had worked with PostGIS 3.0.x versions. I had the same error with 3.1.0alpha3/beta2/rc1 and thought this might an issue with a development version. But now the released version seems to give the same error.
Does anyone know how to fix the error?
(This is under Ubuntu 20.04 LTS).


texts: How to delete a shapefile from a directory using QGIS?

texts:
I am using the following code to clip a shapefile using processing tools via Python and would like to delete the clipped file from the directory.

Here is my code.

from PyQt4.QtCore import *
from PyQt4.QtGui import  *
from qgis.core import *
from qgis.gui import *
from osgeo import *
import processing
import os

boundLayer = "C:\\temp\\bound.shp"
inLayer = "C:\\temp\\layer.shp"
outLayer = "C:\\temp\\outlayer.shp"
processing.runalg("qgis:clip",inLayer,boundLayer,outLayer)

some process

if os.path.exists(outLayer):
    QgsVectorFileWriter.deleteShapeFile(outLayer)

WindowsError: [Error 32] The process cannot access the file because it is being used by another process: 'C:\\temp\\outlayer.dbf'

How can I delete the file?


texts: Obtaining and counting overlapping line segments after adding ST_ClosestPoints with PostGIS

texts:
I have some LineStrings that overlap in some segments, but not necessarily with coinciding vertices, like shown in the example below. With PostGIS I would like to extract all "segments" and count how many LineStrings are overlapping here (1 would mean that nothing overlaps).

As an example, consider the following situation: Two LineStrings (LINESTRING(0 0, 100 40) and LINESTRING(50 20, 150 60)) and a point P (POINT(75 25)).

Easy to see, that the LineStrings overlap between B(50,20) and C(100,40) as shown here:

My approach is:

Extract all (distinct, in the case of duplicate vertices) vertices (A,B,C,D in a first step, see below) from the LineStrings and make a MULTIPOINT of it
ST_Split the LineStrings and count the geometries

I achieve this with the following query (sample data included):

with
    sample(geom, id) as (values
        (st_makeline(array[st_makepoint(0, 0), st_makepoint(100, 40)]), 1),
        (st_makeline(array[st_makepoint(50, 20), st_makepoint(150, 60)]), 2)),

    splitter(geom) as (select st_union(distinct x.geom) as splitter from (select (st_dumppoints(geom)).*  from sample) as x)

select count(y.geom), st_asText(y.geom) as segment from (select (st_dump(st_split(s.geom, splitter.geom))).* from splitter, sample as s) as y group by y.geom

returning the expected result (2 occurences of B-C):

count | segment
----------------------------
1     | LINESTRING(0 0,50 20)
2     | LINESTRING(50 20,100 40)
1     | LINESTRING(100 40,150 60)

To make things a little more interesting, I add a point P' to LineString B-D (cp. above, LINESTRING(50 20, 150 60)), which is the closest point on B-D to P.

Running the same query again (with the modifications to the sample data mentioned above)...

with
    sample(geom, id) as (values
        (st_makeline(array[st_makepoint(0, 0), st_makepoint(100, 40)]), 1),
        (st_addpoint(st_makeline(array[st_makepoint(50, 20), st_makepoint(150, 60)]),
                     st_closestpoint(st_makeline(array[st_makepoint(50, 20), st_makepoint(150, 60)]), st_makepoint(75, 15)),1), 2)),

    splitter(geom) as (select st_union(distinct x.geom) as geom from (select (st_dumppoints(geom)).*  from sample) as x)

select count(y.geom), st_astext(y.geom) as segment from (select (st_dump(st_split(s.geom, splitter.geom))).* from splitter, sample as s) as y group by y.geom

... I expect (not knowing the exact coordinates of P) six segments, with two of them twice:

count | segment
----------------------------
1     | LINESTRING(0 0,50 20)
2     | LINESTRING(50 20, P'x P'y)
2     | LINESTRING(P'x P'y, 100 40)
1     | LINESTRING(100 40,150 60)

... but I get

count | segment
-------------------------------
1     | LINESTRING(0 0,50 20)
1     | LINESTRING(50 20,69.8275862068966 27.9310344827586)
1     | LINESTRING(50 20,100 40)
1     | LINESTRING(69.8275862068966 27.9310344827586,150 60)

So apparently, C does not split B-P'-D (but P' does) and P' does not split A-C (but B does).

Selection of the vertices works fine:

-- create sample data
select distinct st_astext(x.geom) as vertices from (select (st_dumppoints(geom)).*  from sample) as x

vertices
-------------------------
POINT(0 0)
POINT(100 40)
POINT(150 60)
POINT(50 20)
POINT(69.8275862068966 27.9310344827586)

So ST_Split works only when points lie exactly on the LineStrings, but I would expect that P' matches this condition for both LineStrings (or none).

Does adding P' modify the geometry in a way that ST_Split does not work as expected or what might I am doing wrong here? Or is the whole approach likely to be unsuitable? Or is ST_ClosestPoint unsuitable to create P'?


texts: How to filter data by category, then map only maximum-value row for each?

texts:
I'm building a consolidated map of all transit systems in my area in QGIS. I have a GTFS dataset of bus lines, and it includes shapes for every trip on every route in the system (as it should), but that's far more detail than I need, and very slow to render.

I'd like to filter the data to map only the longest possible trip on each route in each direction (inbound/outbound). What's the best way to do this? Line number, route direction, and shape length are all separate columns. I'd prefer to keep the underlying dataset intact if possible.


texts:
I'd like to add a 'Decimal Number (double)' field to an attribute table, however, this option seems to be missing when I go to add a field. I have looked at some videos for adding fields in QGIS and noticed that others have many more field type options than I do, including the 'Decimal number (double)'. Below is a list of the field types I can currently choose from:

Whole number (integer)
Whole number (integer 64bit)
Decimal number (real)
Text (string)
Date
Date and time

Why are my options more limited than others and how can I get the options I need? I'm using QGIS v3.4 Madeira.


texts: Store Graphics in Local when browser refresh

texts:
Using ArcGIS Javascript API, I am adding different graphics to the map. Later I will save them in a server. If I forgot to save, or due to some reason, the page is refreshed before saving, I lose all the graphics. My requirement is to replace the graphics to that particular session.
Is there any possible way to save the graphics?


texts: QGIS 3.6 - Modifying data defined override-expression?

texts:
Is there a way to alter the data defined override expression? I have a points layer which I label with one of its attributes. The attribute table also comes with a field for an offset from point-value for correct displaying, but the format of this field doesn't work with the override on these two points:

Datatype in the table is int, Override is looking for string
Table values range from 1-9, Override is looking for 0-8 (and starts from top left quadrant, instead of bottom left)

I would like to avoid making a copy of the layer and alter the values, as this is a layer that I sync from an external source.
Here is the current expression:
Data defined override
Active: no   (ctrl|right-click toggles)
Expected input:
int
[0=Above Left|1=Above|2=Above Right|
3=Left|4=Over|5=Right|
6=Below Left|7=Below|8=Below Right]
Valid input types:
string
Current definition ():
undefined

And this is what I would like it to be:
Data defined override
Active: no   (ctrl|right-click toggles)
Expected input:
int
[7=Above Left|8=Above|9=Above Right|
4=Left|5=Over|6=Right|
1=Below Left|2=Below|3=Below Right]**
Valid input types:
int
Current definition ():
undefined

Best regards
/Jeston


texts: Is there a way to provide use access to my scripts without showing them to the user? I want to provide free services without giving away my code.


texts: Why do lines appear after running Slope and/or Aspect tools? ArcMap 10.3

texts:
I am attempting to create a predictive model for my study area. I have 10m DEMs which I have merged together using the Mosaic to New Raster tool (Picture 1). I am now trying to use the Slope and Aspect tools on the Merged DEMs. While the tools run successfully, noticeable lines appear after each time I run it (Pictures 2). I can only add two photos, but he Aspect has the same lines in it.

I am at a complete loss. Did I incorrectly merge the DEMs? I went through numerous renditions of merging before finding something that a) didn't completely wash out the terrain, and b) didn't have noticeable merge lines.


texts: I have two datasets- matching years, different sources. They are matchable by FIPS code, but my problem is 1) that the first data set, from which I am drawing my independent variable data, only has the five-digit FIPS code data, and my other data set has information disaggregated down to the eleven digit level. Any thoughts on combining these more disaggregated data?


texts: Using for loop to print last messages from multiple processes?

texts:
Is there a way to use a "for" loop to print the success messages for all of the following geoprocesses? In this example, I am clipping the entire river system of Kansas to a pre-selected ecoregion, and then calculating the total miles of rivers within that ecoregion.

You can see that I have used the "print" function to display a custom completed message after each geoprocess succeeds, but that can be a lot of work; I think a "for" loop would work better.

You'll notice that I've used the arcpy.GetMessageCount() to print the last message, but it only prints the last message of the last geoprocessing item (I believe)

Here's the code for what I've tried. Once again, I would like to get rid of the print statements underneath each geoprocessing step:

import arcpy, sys
from arcpy import env
env.workspace = "C://whatever_workspace"
env.overwriteOutput = True

inFc = "ks_ecoregions.shp"
clipFc = "ks_major_rivers.shp"
outTable = "out_table"

#Processes
try:
    #Check out extension
    arcpy.CheckOutExtension("Spatial")
    print "Extension checked out\n"

    #Make feature layer
    arcpy.MakeFeatureLayer_management(inFc, "ks_ecoregions_layer")
    print "Feature Layer Made\n"

    #Select ecoregion
    arcpy.Select_analysis("ks_ecoregions_layer", "ks_ecoregions_select", where_clause = "US_L3NAME = 'Flint Hills'")
    print "Region selected\n"

    #Perform buffer around the selected ecoregion
    arcpy.Buffer_analysis("ks_ecoregions_select", "buffer_layer", "10 Kilometers")
    print "Buffer created\n"

    #Clip rivers to region
    arcpy.Clip_analysis(clipFc, "ks_ecoregions_select", "riverClip")
    print "Clip created\n"

    #Add table to attributes, converting "kilometers" to "miles
    arcpy.AddGeometryAttributes_management("riverClip", "LENGTH", "MILES_US")
    print "Table added\n"

    #Get summary of total miles of rivers in selected region
    arcpy.Statistics_analysis("riverclip", outTable, [["LENGTH", "SUM"]])
    print "Summation computed\n"

    #Delete intermediate data
    arcpy.Delete_management("ks_ecoregions_layer")
    print "Intermediate feature class deleted\n"
    arcpy.Delete_management("ks_ecoregions_select")
    print "Intermediate feature class deleted\n"
    arcpy.Delete_management("buffer_layer")
    print "Intermediate feature class deleted\n"

    #Get projection of original dataset
    desc = arcpy.Describe("riverClip")
    print "Data information for original projection: \nName: {0} \nData Type: {1} \nProjection: {2}\n".format(desc.name, desc.datasetType, desc.spatialReference.name)

   #Check extension back in
    arcpy.CheckInExtension("Spatial")
    print "Extension checked back in"
except arcpy.ExecuteError:
    pass

severity = arcpy.GetMaxSeverity()

if severity == 2:
    print "Error occurred: \n{0}".format(arcpy.GetMessages(2))
elif severity == 1:
    print "Warning raised: \n{1}".format(arcpy.GetMessages(1))
else:
    count = arcpy.GetMessageCount()
    print arcpy.GetMessage(count-1)


texts: How to create a range domain in ArcGIS with an interval of 10

texts:
I need to create a range domain in a file geodatabase. The user will be able to select a value ranging from 0 and to 200 with an interval of 10. For example, 0, 10, 20, ..., 200.

How can I make and configure a domain type like that?


texts:
There is no way to visualize a convexhull in qgis 1.8 lisboa?
I have tried following views in Postgis 1.5

I was told to use an id column to publish it in Qgis.
If you not use an id column, u cannot load it anyway because of an error.

1.

create or replace view con_hull_baeume3(bid) as
select '0', ST_ConvexHull(ST_Collect(geom)) FROM baeume;

2.

create or replace view con_hull_baeume4(bid) as
select bid, ST_ConvexHull(ST_Collect(geom)) FROM baeume group by bid;

but this one make a convexhull on each of the geometrypoints.
I mean just use 'select ST_ConvexHull(ST_Collect(geom)) FROM baeume;'
works as sql perfectly. But i must define it as a view to load it into a postgis layer of QGIS.


texts: Delete Files in a folder and disregard an existing folder

texts:
I did a lot of geoprocessing steps and now I like to delete all files including result, result, _join, _interand reclass in their name. I can do this with the following

import glob
import os
for fl in glob.glob(r"D:\Users\julia\erste_aufg\delete\*result*" ):
    os.remove(fl)

for f2 in glob.glob(r"D:\Users\julia\erste_aufg\delete\*result*" ):
    os.remove(f2)

for f3 in glob.glob(r"D:\Users\julia\erste_aufg\delete\*_join*" ):
    os.remove(f3)

for f4 in glob.glob(r"D:\Users\julia\erste_aufg\delete\*_inter*" ):
    os.remove(f4)

for f5 in glob.glob(r"D:\Users\julia\erste_aufg\delete\*reclass*" ):
    os.remove(f5)

see here.

But before deleting I created a resultfolder in the folder(workspace) with the files and copied the results of the geoprocessing steps which I need furthermore to this resultfolder. Because of this resultfolder the code to delete is not working any more. What can I do that the delete-code disregards the resultfolder?

(The code down is just to make the question more clear)

import arcpy
import os
from arcpy import env
env.workspace = r"env.workspace = r"D:\Users\julia\urban_A_sel""
#new folder
# Set local variables
out_folder_path = env.workspace
out_name = "resultfolder"

# Execute CreateFolder
arcpy.CreateFolder_management(out_folder_path, out_name)
newpath= env.workspace + "\\" + out_name
out_gdb_path = newpath
arcpy.CreateFileGDB_management(out_gdb_path, "output.gdb")
outWorkspace= newpath +"\\"+"output.gdb"
resultList = arcpy.ListFeatureClasses("*_result*")
arcpy.FeatureClassToGeodatabase_conversion(resultList, outWorkspace)


texts:
I have two raster layer, the first one r is my zonal mask as a raster grid, and I created as:

ext <-extent(Yraster)
r <- raster(nrows=32,ncols=32,ext)
r[] <-(1:ncell(r))

and project it as the 'Yraster' one

projection(r)<-"+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

  > r
class       : RasterLayer
dimensions  : 32, 32, 1024  (nrow, ncol, ncell)
resolution  : 428.465, 430.3401  (x, y)
extent      : 346090.2, 359801.1, 8783943, 8797714  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : layer
values      : 1, 1024  (min, max)

> Yraster
class       : RasterLayer
dimensions  : 459, 457, 209763  (nrow, ncol, ncell)
resolution  : 30.00192, 30.00192  (x, y)
extent      : 346090.2, 359801.1, 8783943, 8797714  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : acre_2010_non1
values      : 1, 1  (min, max)
Raster Attribute Table
     fields : ID Rowid    COUNT
        min :  1     0 13821927
        max :  1     0 13821927

but when I do the zonal statistics, an error pop out:

> list <- zonal(Yraster,r, na.rm=TRUE)
Error in compare(c(x, z)) : ncols different

How can I make the zonal statistics work?


texts:
When we run Generate Near Table tool in ArcMap 10.4.1, we get IN_FID and NEAR_FID.
I ran it through my point features and have my data.
But, I also want coordinates of IN_FID points and NEAR_FID points in the same table.
Suggest a clue. I use ArcGIS 10.4.1


texts: How do I separate/visualize polylines that overlap? (Show them next to each other)

texts:
I have a network of polylines between multiple nodes. All nodes are connected with lines in a two-way system (A->B and B->A). I want to show all the lines at the same time, but the problem is that they overlap and only one correlation are shown (A->B but not B->A). This means that the correlation just seems to be one-way instead of two-way.

How do I show them next to each other using QGIS?


texts: What is the most suitable way to get the boundary of a state or a country or the whole world dynamically?

texts:
I am working on a project where I have to create filled map charts based on the data client has provided. The data has latitude and longitude information. Based on the selected level (city or state or country) I have to get the boundary shapes of those locations and then filled them with appropriate color.

I have tried nominatim framework api. They are giving response in geoJson format which I am rendering in client side using leaftlet.js. As for a single country the size of geoJson is around 1-1.5 MB, I cannot use this approach.

I tried searching for any api service which can serve boundary data in geoJSON but still no luck. I have seen people suggesting to download and store shapefile(shp, shx), and convert them to geoJSON using QGIS manually. But as I am going to use it for a web-application and need to get/generate shapefiles on request from client, I don't know how this kind of approach will fit in it.
Suggest me some tools and libraries which can be used to convert this data to a supported format for a web application.

I know some of the things is sounding naive, as I am very new to this field. Let me know if you need any clarity on anything I said about my project.

Please help me to find a better approach to solve this use-case.


texts: Select interaction for each layer is not working in OpenLayers

texts:
I have multiple layers in my OpenLayers application. I need layer based specific info Windows. Each layer will have its own select interaction. This is my problem. I do not want to add layers in one select interaction, because it is not my solution.

var feature1 = new ol.Feature(new ol.geom.Point(ol.proj.fromLonLat([21.54967, 38.70250])));
var vectorSource1 = new ol.source.Vector({ features: [feature1] });
var vectorLayer1 = new ol.layer.Vector({ source: vectorSource1 });
var selectClick1 = new ol.interaction.Select({
    condition: ol.events.condition.click,
  layers: [vectorLayer1]
});
selectClick1.on("select", function(e) { alert(1); })

var feature2 = new ol.Feature(new ol.geom.Point(ol.proj.fromLonLat([21.84967, 39.70250])));
var vectorSource2 = new ol.source.Vector({ features: [feature2] });
var vectorLayer2 = new ol.layer.Vector({ source: vectorSource2 });
var selectClick2 = new ol.interaction.Select({
    condition: ol.events.condition.click,
  layers: [vectorLayer2]
});
selectClick2.on("select", function(e) { alert(2); })

var feature3 = new ol.Feature(new ol.geom.Point(ol.proj.fromLonLat([21.44967, 37.70250])));
var vectorSource3 = new ol.source.Vector({ features: [feature3] });
var vectorLayer3 = new ol.layer.Vector({ source: vectorSource3 });
var selectClick3 = new ol.interaction.Select({
    condition: ol.events.condition.click,
  layers: [vectorLayer3]
});
selectClick3.on("select", function(e) { alert(3); })

var map = new ol.Map({
  target: 'map',
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    }),
    vectorLayer1,
    vectorLayer2,
    vectorLayer3
  ],
  view: new ol.View({
    center: ol.proj.fromLonLat([21.54967, 38.70250]),
    zoom: 8
  })
});

map.addInteraction(selectClick1);
map.addInteraction(selectClick2);
map.addInteraction(selectClick3);

The working code is here.


texts:
I want to do something like this using Python and QGIS:

And I want to save what user input in each field in 4 variables.
Here is what I tried but it is not giving what I want:
import sys

class Dialog(QDialog):
    NumGridRows = 3
    NumButtons = 4

    def __init__(self):
        super(Dialog, self).__init__()
        self.createFormGroupBox()

        buttonBox = QDialogButtonBox(QDialogButtonBox.Ok | QDialogButtonBox.Cancel)
        buttonBox.accepted.connect(self.accept)
        buttonBox.rejected.connect(self.reject)

        mainLayout = QVBoxLayout()
        mainLayout.addWidget(self.formGroupBox)
        mainLayout.addWidget(buttonBox)
        self.setLayout(mainLayout)

        self.setWindowTitle("Form Layout - pythonspot.com")

    def createFormGroupBox(self):
        self.formGroupBox = QGroupBox("Form layout")
        layout = QFormLayout()
        layout.addRow(QLabel("Street:"), QLineEdit())
        layout.addRow(QLabel("Apartment:"), QLineEdit())
        layout.addRow(QLabel("Post code:"), QLineEdit())
        layout.addRow(QLabel("City:"), QLineEdit())
        self.formGroupBox.setLayout(layout)

app = QApplication(sys.argv)
dialog = Dialog()
sys.exit(dialog.exec_())

But gives this:


texts: How to preserve NetCDF dimensions with gdal.warp

texts:
I have a netcdf file that is structured with precipitation grid values with dimensions (time, x, y)

When I perform gdal.warp (in python if that matters) I get a new file that has multiple bands that have dimensions (x, y), with the number of bands = len(time) of the previous file and time is not part of the new file.

Is there a way to preserve the dimensions, precipitation name, and time in the new file?

import gdal

warped = gdal.Warp(
                destNameOrDestDS=destNameOrDestDS,
                srcDSOrSrcDSTab=srcDSOrSrcDSTab,
                srcSRS=srcSRS,
                dstSRS=dstSRS,
                dstNodata=float(srcNodata),
                format="NETCDF",
                xRes=xRes,
                yRes=yRes,
                targetAlignedPixels=targetAlignedPixels,
            )

So my source netcdf file looks like:

crs
QPE
time
x
y
z

where crs is the coordinate system and QPE is precipitation and then my output looks like:

Albers
Band1
Band2
Band3
x
y

I would like to preserve the original structure so combine the Bands into a single 3 dimensional array with time as the 3rd dimension.


texts:
I would like to add a white halo to labels using ArcObjects. I think I need to reference the IFormattedTextSymbol interface, but I'm not sure how to implement this into my current code (see below).

Dim pDoc As IMxDocument
Set pDoc = ThisDocument
Set pMap = pDoc.FocusMap
' Get the selected layer
Dim pLayer As IGeoFeatureLayer
Set pLayer = pMap.Layer(0)

Dim pFc As IFeatureClass
Dim strOIDName As String
Set pFc = pLayer.FeatureClass
strOIDName = pFc.OIDFieldName

' Get the selected features from the layer
Dim pFSel As IFeatureSelection
Set pFSel = pLayer
Dim pSelSet As ISelectionSet
Set pSelSet = pFSel.SelectionSet
Dim pFCur As IFeatureCursor
pSelSet.Search Nothing, False, pFCur

' Loop through the selected features and create a label
' expression for the selected features

Dim pFeat As IFeature
Dim strSql As String
Set pFeat = pFCur.NextFeature

Do While Not pFeat Is Nothing
   If strSql = "" Then
     strSql = strOIDName & " = " & pFeat.OID
   Else
     strSql = strSql & " or " & strOIDName & " = " & pFeat.OID
   End If
   Set pFeat = pFCur.NextFeature
Loop
Debug.Print strSql

' Get AnnotateLayerPropertiesCollection from layer
Dim pAnnoLayerPropsColl As IAnnotateLayerPropertiesCollection
Set pAnnoLayerPropsColl = pLayer.AnnotationProperties
Dim pAnnoLayerProps As IAnnotateLayerProperties
pAnnoLayerPropsColl.QueryItem 0, pAnnoLayerProps, Nothing, Nothing
pAnnoLayerProps.Class = "LabelSel"
pAnnoLayerProps.WhereClause = strSql
Dim aAnnoVBScriptEngine As IAnnotationExpressionEngine
Set aAnnoVBScriptEngine = New AnnotationVBScriptEngine

Dim pTextsymbol As ITextSymbol
Set pTextsymbol = New TextSymbol

Dim pFont As IFontDisp
Set pFont = pTextsymbol.font

Dim pFormattedtxtsym As IFormattedTextSymbol
pFormattedtxtsym.ShadowColor = True

pFont.Bold = True
pFont.size = 10
pFont.Weight = 2
pTextsymbol.font = pFont

Dim pLELayerProps As ILabelEngineLayerProperties
Set pLELayerProps = pAnnoLayerProps
Set pLELayerProps.ExpressionParser = aAnnoVBScriptEngine

pLELayerProps.Expression = "[Longlabel] & vbnewline & [Latlabel] & vbnewline & [Elevlabel] & _ vbnewline &[Driftlabel]"

Set pLELayerProps.Symbol = pTextsymbol

' Display the lables
pLayer.DisplayAnnotation = True
' Refresh the Data Frame
pDoc.ActiveView.Refresh


texts: Where does the Geonames admin_code1 for Russia come from?

texts:
Looking at the data dump of post codes for Russia (country downloads | Russian download link) I can't understand what the admin_code1 corresponds to.
Looking at the ISO3166 sub regions, I see the codes are characters (whereas in geonames the codes are numeric).  Looking at Wikipedia I can't find a numeric reference for regions anywhere.  I was tempted to order alphabetically and take the row number, but some numbers are missing (2 for example).
In the readme.txt geonames lists admin_code1  as: 1. order subdivision (state) varchar(20), but this is just generic to all countries and gives no source.  I've also checked the Russia page, which links to the St. Petersburg's administrative website, but there's seemingly nothing useful here.
What do the numbers here represent?  I need to link the ISO 3166 subdivisions to the countries in geonames, thus I need to create a lookup, is there an existing resource, or way to do this?


texts: Batch Calculate of loaded Vector Layers using PyQGIS?

texts:
I have a series of polygon that was already batch loaded in my QGIS desktop using this code (courtesy of this post):

import os

source_dir = "/data/brent/LUCID/admin_boundaries_clup/sen2agri_all_muni"
for files in os.listdir(source_dir):
    if files.endswith(".shp"):
        # create vector layer object
        vlayer = QgsVectorLayer(source_dir + "/" + files, files, "ogr")
        print(files)
        # add the layer to the registry
        QgsMapLayerRegistry.instance().addMapLayer(vlayer)

The next steps I want to do to each recently added vectors are:

1) Enter into editing mode. (In preparation for number 2)

2) Create new field, field name as "AREA", field type = real, output length = 10, precision = 2.

3) Calculate area using this string (in field calculator): '$area'

4) Save the edits.

Is it possible to insert codes from the example above?


texts: How to insert data into postgis database with phpPgAdmin?

texts: I have 50 line objects in .shp format in folder on my computer. I need to insert it to postgis database on the distant server.. When I log in in phpPgAdmin, it says it runs on localhost:5433. What is the quick way to insert my .shp data into database?


texts: Specify an additional plugin folder that's located in the local network

texts:
I'm creating several plugins for an organization, that uses a local network to share workflow. But isn't a go choice install and upgrade the plugins for every coworker computer.

How could I create a folder into the network and then specify this path to QGIS? Is there any field to add complements folders such as exists for C++ complements?


texts: Getting ListFields from Autocad DWG using ArcPy?

texts:
I am trying to get LisfFields from AutoCad DWG, to arcpy Add-In, i have tried

    mxd = arcpy.mapping.MapDocument("CURRENT")
    Layers = arcpy.mapping.ListLayers(mxd)
    for Layers in Layers:
        if Layers.name==Target_Layer:
            fieldList = arcpy.ListFields(Target_Layer)
            for field in fieldList:
                fieldName=str(field.name)

It works for Feature Class Layers but not for DWG.

Is there other way to do this?

I have also tried MakeFeatureLayer before LisfFields and it did not work.


texts: after reinstall geoserver how to locate old geoserver files

texts: i have published lots of Arch GIS map on Geo server software. but few days ago because of some problem in my system i have uninstalled Geo server 2.1.4 from my system. but i kept a copy of my Geo server folder. Now i have again installed the same version of Geo server in my system. now the question is how to relocate my work folder in Geo server. i don't want to publish again all those files.


texts: Push JSON to an External WebSocket Error ESRI Geoevent

texts:
Has anyone successfully configured a Push To An External Websocket with ESRI geoevent?  I have an input a connector, service connector, and an output, however when I start the output I immediately get the error:

Status: Error - Could not startup the WebSocket Transport. An initialization problem has occurred. Error: Connection refused: no further information.

I increased the buffer and text bytes to 102400 in the settings  and used the default output settings:
ws://localhost:80/ws also tried ``wss://localhost:6180/ws2`

I am outputting a simulation stream to troubleshoot and I can output UDP so I know GeoEvent is working in some capacity.

I have tested this both on linux redhat 7.1 and windows 7 and the results are the same.  I can attempt any suggestions provided on either OS. And I am running ArcGIS Server 10.3.1 on local machines.


texts: Printing maps from OpenLayers with a Google Maps base layer

texts:
I need to implement an automatic pdf report generator on a web site. In the generated pdf file I need to insert images from maps extracted from OpenLayers with a Google Map base layer.

On top of Google Map I have an overlay with some custom OWS (WMS/WFS etc).

Is this possible? And is it allowed? I was reading the ToS of Google, and it is not totally clear to me.


texts: Using Geodata Service in ArcGIS API for JavaScript?

texts:
Can I use Geodata service in ArcGIS JavaScript API 3.1?

Background: I have simple map service and the all attribute data is stored in SQL Server database. I want to perform query operation on map service in JSAPI (and data should come from database.)

    I have only database access and I cannot create the Geo-database.
    I can't attach the data to shapefile (it's too large)

I searched on Google and asks few people so they told me that publish the data as "geodata service" but according to API documentation JSAPI (3.1) do not support Geodata service.

Also in JSAPI examples it's not mentioned clearly.

So what should I do?

Software: ArcGIS Server 10.1, ArcGIS Java script API 3.1, ArcSDE


texts: Converting raster from WGS_1984_plate_carree to WGS_1984 in ArcGIS Desktop?

texts:
I am a beginner using ArcGIS and I  downloaded a TIFF which is in WGS_1984_plate_carree coordinate system. This image covers all the world and in this coordinate system the extent goes from -180 (left) to 180 (right) and from -90 (bottom) to 90 (top). I need to convert it to WGS_1984 without changing its extent but when i try to convert it, using the project raster tool, i get this extent:

Top      0,000808483755708

Left     -0,00161696751142

Right    0,00161696751142

Bottom   -0,000808483755708

I know that plate carree is a projection while WGS_1984 is a geographic coordinate system but i don't know how to solve this problem


texts:
I have attached geodatabase "Riverside" to databases list in MS SQL SERVER 2012 Enterprise edition  :

After that I created two logins using :

   USE [master]
GO
CREATE LOGIN [sde] WITH PASSWORD='sde_class', DEFAULT_DATABASE=[master], DEFAULT_LANGUAGE=[us_english], CHECK_EXPIRATION=OFF, CHECK_POLICY=OFF
CREATE LOGIN [GISADMIN] WITH PASSWORD='gisadmin', DEFAULT_DATABASE=[master], DEFAULT_LANGUAGE=[us_english], CHECK_EXPIRATION=OFF, CHECK_POLICY=OFF

Then  I tried to add a database connection from catalog window within ArcMap ,but the connection failed :

After that  I went to security folder inside the "Riverside"  I found a user with name "GISADMIN" . I thought I could delete this user and create another user with same name . The deletion failed ,

my last option was mapping this user"GISADMIN" to the login "GISADMIN" .

use RIVERSIDE
go
exec sp_change_users_login 'GISADMIN', 'gisadmin'
go

Nothing changed!
I'm using ArcGIS 10.3 running on Windows 7 machine.


texts: Changing raster values in specific zone only

texts:
I would like to replace raster values only in specific zone of the raster. Say this is my example raster:

library(raster)
r <- raster(ncol=1000, nrow=1000)
r[] <- runif(ncell(r),0,1)
extent(r) <- matrix(c(0, 0, 1000, 1000), nrow=2)

What I would like to do now, is to change all values that equal to val <- 0.8and replace with valnew <- 0.2.
The extent of raster where my function should search and replace values is following (only this square of raster should be searched for val values in raster and replaced with newval)

extentnewvals <- matrix(c(300, 300, 500, 500), nrow=2)

Any ideas how to approach it?


texts: Error 2 Creation of data source "A"."B" failed: ERROR: relation "B" already exist

texts:
Does anyone know what I may be doing wrong?

Foolish question from a newbie perhaps but as I reckon it may be the code for table "B" that I created, how do I edit the code in pgAdmin3?

The key aim of this is to migrate data from one table in MapInfo into an existing table with a completely different table structure in QGIS.

code for current table below:

-- Table: A."B"

-- DROP TABLE A."B";

CREATE TABLE A."B"
(
  "UID" integer NOT NULL,
  roadname character varying(100),
  CONSTRAINT tmp_ga2_pk PRIMARY KEY ("UID")
)
WITH (
  OIDS=FALSE
);
ALTER TABLE A."B"
  OWNER TO postgres;


texts: v.net.allpairs - Missing parameter value: WHERE condition SQL statement without where keyword

texts:
v.net.allpairs is not running in QGIS 2.18.27 (but runs in QGIS 2.14).
QGIS 2.18.27 has responded: Missing parameter value: WHERE condition SQL statement without where keyword.
How to fix that?


texts: Is Data Interoperability Extension for ArcGIS Server required to extract data in .dgn format?

texts:
I am trying to create a Geoprocessing Task in my .mxd to in order to extract data to .SHP and .DGN formats which will be utilized in a data extraction tool for a web application.

Right now, there is a version conflict between ArcMap 10 / ArcGIS for Server 9.3 which will be solved next week the server is upgraded to 10.

Is the Data Interoperability (DI) extension required for ArcGIS Desktop to create the GP task? More importantly, is the DI extension/license required for ArcGIS Server 10 for the .DGN extraction tool to work on the application?