texts: Using PostGIS to calculate distance of route over large area

texts:
I am new to PostGIS and looking to measure the distance of the Pacific Crest Trail which spans western North America from Canada to Mexico.

I have a data set containing the lat lon and dates from GPX data taken on the trail from a hiker. I also have a shapefile of the entire trail that's fairly high resolution. What I'd like to do is to extract the total mileage gained for each of the GPS locations from computing the length of the trail data and then intersecting the two.

I think the steps to doing this as are follows:

process the GPS location data to get lat lon positions at weekly intervals
with additional special events like summiting Mt Whitney, trail start, end, etc.(this part is finished)
snap the lat lon GPS points to the nearest PCT shapefile nodes.
split the PCT shapefile into segments based on an intersection with the GPS points.
convert the total mileage from meters to miles.
total the distance gained by each GPS point as the dates increase.

Some questions I have are:

What projection or method of measuring distance would be most accurate for such a large north - south distance?

I'm thinking of transforming the data into a custom equidistant projection but am unsure about setting the parameters for yielding the most accurate distance measurement.

Or does PostGIS have other functionality that would yield accurate distance measurements without having to transform the data?

Does this sound like a reasonable approach?

Would you have any suggestions or tips that might help?


texts:
I write a code pack all table in Map Window, after pack all table, some tables are removed from map window,

How do I known which table is pack, and open table to map window in the order as original layer.


texts: D3 for maps---at what stage to bring in data to the geo?

texts:
I'd like to map a world choropleth for display with D3, a la:

http://bl.ocks.org/mbostock/4206573
http://bl.ocks.org/jasondavies/4188334

I have a dataset I'd like to display that's keyed to ISO-alpha-3 keys.  So...

danger.csv
iso,level
AFG,100
ALB,0
DZA,12

etc.

Following the instructions on topojson, I know I can do...

wget "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_0_countries.zip"
unzip ne_50m_admin_0_countries.zip
ogr2ogr -f "GeoJSON" output_features.json ne_50m_admin_0_countries.shp -select iso_a3
topojson -o topo.json output_features.json --id-property iso_a3

to produce a worldmap json that is ID'd by ISO3.

My question is: at what point in the workflow should I merge the data from danger.csv onto the geo data?  I had previously worked with qGIS as a GUI, but where /should/ the merge happen?  In the .shp?  After the ogr2ogr?  Dynamically in the browser after the topojson shrink (like here http://bl.ocks.org/mbostock/4060606 http://bl.ocks.org/mbostock/3306362)?

I'm pretty good with python, but pretty new to javascript, and find myself copying and pasting Bostock examples more than actually being a generative coder there.

(I also have a related, but more involved follow-up on Stackoverflow that maybe I should migrate here:  https://stackoverflow.com/questions/18604877/how-to-do-time-data-in-d3-maps)


texts: Getting lat/long for centroids using QGIS

texts:
I have a series of polygons (lakes) in Kansas, Oklahoma and a few in Texas.
I created centroids of each lake, which I need for distance measurements.
These centroids have derived lat/long information when they are clicked, but I am unable to get the same information added to the attribute table.
After doing research, I see that a number of people have had similar problems with the lat/long.  I did manage to get some sort of location in the attribute table, similar to this (Adding coordinates to attribute table), but I'm unable to get them converted to actual longitude/latitude.


texts: "ImportError: no module named qgis.core" in standalone app

texts:
I am developing a standalone app (in Windows) using libraries from QGIS. When I run the script in cmd, I get the ImportError: no module named qgis.core. I read in several posts that I need to configure the PYTHONPATH and PATH but I also read that with OSGEO4W installation there is no need to configure.

How do I do that?

I am very confused and I try several ways but without success. Now I am reinstalling QGIS thorugh OSGEO4W.

So, I need to import QGIS libraries in a script outside QGIS. This script runs inside Python27 folder that I have already installed for other app.

Is it possible to use the QGIS libraries through this directory?

If not what is the solution?

I resolved the import qgis.core error, just following the instructions of how to import qgis modules outside qgis. Now the error is different: ImportError: DLL load failed: %1 is not a valid Win32 application. And the error is related with import sip. I search and find that I use python from OSGEO4W (64bits) and external python (32 bits).

How I can solve it?


texts: Automatic geo-referencing options - mosaic dataset footprints?

texts:
I have 50 Jpgs that represent a route and I want to automatically georeference them. I have previously georeferenced 50 jpgs that exactly match the image extents and created a mosaic dataset. I would simply copy and rename the .jgwx/xml files and this almost works, but I am assuming because the resolutions of the jpegs are different this is offsetting them slightly.

Is there a way of automating this possibly using previous mosaic footprints?

I also can't re-export them from source to match resolutions.


texts: Split polygon with selected line from another layer in QGIS

texts:
I am trying to split a selected polygon ("regions") with a selected line from another layer ("streets") in QGIS 3.8.3.
In the example below, I want to split the selected polygon coloured in yellow with the street pointed at by the blue arrow.

After reading the post Can't Find the Split selected features with selected line from another layer tool, I realized I can't find it either. I installed the plugin, but the tool seems to be missing from my toolbar:


texts: Getting GPS coordinates perpendicular and parallel to AB?

texts:
I already posted a question at Stack Overflow

I got the solution from @Fang but the problem is that the coordinates I need to use are GPS coordinates and apparently GPS needs a different formula than cartesian.

I am working with Google Maps APIv3 polygons and having coordinates of AB (by drawing them with the polygon tool) I click on position C which I need to be moved to D which is perpendicular to AB and CD is parallel to AB

Having

A = 50.88269282423443,6.0036662220954895
B = 50.882753744583226,6.003803014755249
C = 50.88252571592428, 6.003832183778286

- D is perpendicular to AB
- CD is parallel to AB

What would be the formula to get D?


texts: C# library for converting geom field of SQLServer into set of google latlng

texts: Is there any C# library which can convert the geometry field into set of google latlngs? I have stored a shapefile into SQLserver table. From Sqlserver i want to read the geometry field data and want to convert this geom data into set of google latlngs in my C# application. Finally i want to save these google latlngs into another table in SQLServer.


texts: MapServer and OpenLayers: Properly using MAXSCALEDENOM in mapfile (google baselayer)

texts:
My base layer: google map

My points on the map: from an airports.db file, set up as a WFS layer in the OpenLayers code.

As of now:

map loads, points in the correct places on the map.

My goal:

when the page loads in the browser, the points should not display on the initial zoom. When you zoom in to the second level, they should.

I've successfully done this before using just an OpenLayers map, but never a google map as a base layer.

For MAXSCALEDENOM I've tried all sorts of values, but the points are always displaying.

Here is the Airports LAYER in the mapfile.

LAYER
            NAME Airports
            DEBUG 5
            METADATA
                    "wfs_srs" "EPSG:3857"
                    "wfs_title" "Airports"
                    "gml_featureid" "Name"
                    "gml_include_items" "all"
                    "wfs_enable_request" "*"
                    "wfs_encoding" "UTF-8"
                    "wfs_extent" "-180 -90 180 90"
            END
            TYPE point
            STATUS ON
            MAXSCALEDENOM 20000
            LABELITEM 'Name'
            CONNECTIONTYPE OGR
            CONNECTION '/tmp/navhome/navdata/airports.db'
            DATA 'select "POINT(" || LONG || " " || LAT || ")" as WKT_GEOMETRY, * FROM airports WHERE name LIKE "S%"'
            #DATA 'select "POINT(" || LONG || " " || LAT || ")" as WKT_GEOMETRY, * FROM airports WHERE ICAO != "ZZZZ"'
            PROJECTION
                    "proj=latlong"
                    "ellps=WGS84"
                    "datum=WGS84"
            END # PROJECTION
    END # LAYER

Here is the OpenLayers code for the WFS layer:

var wfs = new OpenLayers.Layer.Vector("Airports", {
            styleMap: new OpenLayers.StyleMap({
                    externalGraphic: "red_pin.png",
                            graphicWidth: 18,
                            graphicHeight: 18
            }),
            wrapDateLine: true,
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.WFS({
                    srsName: "EPSG:3857",
                    version: "1.0.0",
                    url: "http://my_url/cgi-bin/mapserv.exe?map=/tmp/navhome/navcore/airports2/airports_viewer.map&SERVICE=WFS&VERSION=1.1.0&srsName=EPSG:3857",
                    featureType: "Airports"
            })
    });
    map.addLayer(wfs);

Any input on this matter is appreciated. I assume it has something to do with the extent, but I'm not totally sure.


texts: Update field properties of a layer definition with ArcGis Python API

texts:
We have a fields in a large amount of AGOL hosted feature layers that we want to make un-editable and others that we do not want visible.

Rather than going into the service definition for each and every one to manually change "editable": true to "editable": false for many fields in many feature layers, and then also do the same for "visible", we were hoping to automate the process via the Python API (which we are new to).

So far we have:

from arcgis.gis import gis
gis = GIS(.....login info......)
mpfs = gis.content.get(......id......)
mpfl = mpfs.layers[0]
replacementCopy = mpfl.properties.fields
for field in replacementCopy:
    if field.name in (list of names):
        field["editable"] = False
for field in replacementCopy:
    if field.name in (list of names):
        field["visible"] = False

But the field properties seem to be missing the "visible" property, and we are not sure how to go about enacting the changes once we have the update field properties.


texts: How to customizing the default map scale list in ArcGIS for all users?

texts:
How do I customizing the default map scale list in ArcGIS? I can set a list from a text file.

The list needs the following format (from comment of @Martin) for example:

1000,000000
10000,000000
25000,000000
50000,000000
100000,000000
250000,000000
500000,000000
750000,000000
1000000,000000
3000000,000000
10000000,000000

How do I set this list as default list for all users on all workstations?


texts:
I have around 1200 tiff files (example below), and corresponding .tfw files and .aux files with sonar data that I'd need to merge into one file and I'm struggling to do so.

I browsed some topics about that problem but doing it manually is just not an option in this case and none of the other solutions helps my case.

Does anyone know if QGIS is capable of merging tiff files and ignoring the white color (i.e. treat it as transparent)?


texts: Modifying coordinate values of MultiPolygon using GeoTools?

texts:
I can already get the vertex coordinates of the polygon, but I want to modify the coordinates and write it back to the source file.
This is part of the code I read the coordinates.

try{
      while (iter.hasNext()) {
          SimpleFeature sf = iter.next();
          Collection<Property> property =((ComplexAttribute)sf).getProperties();
          Iterator<Property> iterP = property.iterator();
          while (iterP.hasNext()) {
              Property pro = iterP.next();
              if (pro.getValue() instanceof MultiPolygon) {
                  gisObjectList.add((T) new GisMultiPolygon ((MultiPolygon)pro.getValue(), (SimpleFeature) sf));
                   Geometry gc = (MultiPolygon) pro.getValue();
                  Coordinate[] points = gc.getCoordinates();
              }
          }
      }
  }

I created a new shapefile from the original shapefile, but there was a runtime error.
error is org.geotools.data.directory.DirectoryDataStore cannot be cast to org.geotools.data.shapefile.ShapefileDataStore.

public static void transShape(File shpFile, File shpFile2) {
    try {
        //源shape文件
        ShapefileDataStore shapeDS = (ShapefileDataStore) new ShapefileDataStoreFactory().createDataStore(shpFile.toURI().toURL());
        //创建目标shape文件对象
        Map<String, Serializable> params = new HashMap<String, Serializable>();
        FileDataStoreFactorySpi factory = new ShapefileDataStoreFactory();
        params.put(ShapefileDataStoreFactory.URLP.key, shpFile2.toURI().toURL());
        ShapefileDataStore ds = (ShapefileDataStore) factory.createNewDataStore(params);
        // 设置属性
        SimpleFeatureSource fs = shapeDS.getFeatureSource(shapeDS.getTypeNames()[0]);
        //下面这行还有其他写法，根据源shape文件的simpleFeatureType可以不用retype，而直接用fs.getSchema设置
        ds.createSchema(SimpleFeatureTypeBuilder.retype(fs.getSchema(), DefaultGeographicCRS.WGS84));
        //设置writer
        FeatureWriter<SimpleFeatureType, SimpleFeature> writer = ds.getFeatureWriter(ds.getTypeNames()[0], Transaction.AUTO_COMMIT);
        //写记录
        SimpleFeatureIterator it = fs.getFeatures().features();
        try {
            while (it.hasNext()) {
                SimpleFeature f = it.next();
                SimpleFeature fNew = writer.next();
                fNew.setAttributes(f.getAttributes());
                writer.write();
            }
        } finally {
            it.close();
        }
        writer.close();
        ds.dispose();
        shapeDS.dispose();
    } catch (Exception e) { e.printStackTrace();    }
}


texts: Cannot input new data into FeatureDataForm

texts:
I am wanting to update the fields on my server using a FeatureDataForm. I have followed the documentation, but it still does not work. It also does not have the buttons in the form that say 'OK' and 'Delete' like the documentation does. Any help will be appreciated.

XAML (FeatureDataForm):

<Border x:Name="FeatureDataFormBorder" Visibility="Collapsed"
                HorizontalAlignment="Right" VerticalAlignment="Top" Margin="10,10,10,0" Width="300" Height="300" >
                <Border.Effect>
                    <DropShadowEffect Color="Black" Direction="-45" BlurRadius="20" Opacity=".75" />
                </Border.Effect>
                <esri:FeatureDataForm x:Name="MyFeatureDataForm"  Foreground="Black" BorderBrush="{StaticResource PanelGradient}" BorderThickness="12"
                                         IsReadOnly="False" LabelPosition="Left" />
            </Border>

XAML(MyMap):

<esri:Map x:Name="MyMap" Background="#FFE3E3E3" WrapAround="True"  MouseRightButtonDown="MyMap_MouseRightButtonDown" Extent="-9834972.92753924,4441899.425293319,-9833977.88119163, 4442762.485358352" Loaded="MyMap_Loaded">
                <esri:Map.Layers>
                    <esri:ArcGISTiledMapServiceLayer ID="MyLayer" Url="http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer" />
                    <esri:FeatureLayer ID="WaterValves" Url="http://localhost:6080/arcgis/rest/services/CalvertCity_Test_2/MapServer/0" Renderer="{StaticResource MySimpleRenderer}"
                                                  OutFields="*" MouseLeftButtonUp="FeatureLayer_MouseLeftButtonUp" MouseLeftButtonDown="FeatureLayer_MouseLeftButtonDown_1" DisableClientCaching="True"  Mode="OnDemand"/>
                    <esri:GraphicsLayer ID="CandidateGraphicsLayer">
                        <esri:GraphicsLayer.MapTip>
                            <Grid>
                                <Rectangle Stroke="Gray"  RadiusX="10" RadiusY="10" Fill="#77FF0000" Margin="0,0,0,5" >
                                    <Rectangle.Effect>
                                        <DropShadowEffect/>
                                    </Rectangle.Effect>
                                </Rectangle>
                                <Rectangle Fill="#DDFFFFFF" Stroke="DarkGray" RadiusX="5" RadiusY="5" Margin="10,10,10,15" />
                                <StackPanel Orientation="Vertical" HorizontalAlignment="Center" Margin="30,20,30,30">
                                    <TextBlock Text="{Binding [Address]}" HorizontalAlignment="Left" Foreground="Black" />
                                </StackPanel>
                            </Grid>
                        </esri:GraphicsLayer.MapTip>
                    </esri:GraphicsLayer>
                    <esri:ArcGISDynamicMapServiceLayer ID="CalvertCity" Url="http://localhost:6080/arcgis/rest/services/CalvertCity_Test_2/MapServer" />
                </esri:Map.Layers>
            </esri:Map>

C#:

private void MyMap_Loaded(object sender, RoutedEventArgs e)
        {
            MyFeatureDataForm.FeatureLayer = MyMap.Layers[1] as ESRI.ArcGIS.Client.FeatureLayer;
        }

        private void FeatureLayer_MouseLeftButtonDown_1(object sender, GraphicMouseButtonEventArgs e)
        {
            e.Graphic.Selected = !e.Graphic.Selected;
            if (e.Graphic.Selected)
                MyDataGrid.ScrollIntoView(e.Graphic, null);
        }

I believe this tells me that the map is editable:


texts:
I'm trying to carry out weather conditions interpolations. I have information of 92 weather stations for 800,000 km2. As you can imagine the extension of land is huge and I need some kind or weather models to predict wind velocity and direction and solar irradiation for points in between the stations via interpolations. This means if I have a weather station in A and another one in C, I can know what is happening in B with weather interpolation models.

Do you know if QGIS does something alike? (maybe there is a plugin to download?) Are you familiar with any models that could help me out?


texts: Should I be able to open personal DB converted to MSSQL as a MSSQL Layer?

texts:
I'm feeling like this should work, but I'm banging my head against the keyboard, the documentaion and the google as I can't find any examples to follow for MSSQL.

I have 3rd party data only in an ESRI Personal DB (.MDB). I converted that to MSSQL because I believe there's no way to open that in QGIS.

The ACCESS to SQL conversion ran fine and I can see the tables and data in the SQL Server DB. I setup a User DSN to the DB and when I add a new MSSQL Layer connection, the Test Connection works.

However when I attempt to Connect to the DB in the Add MSSQL Table Dialog I get the following error:

[Microsoft][ODBC Driver 11 for SQL Server][SQL Server]Invalid object name 'geometry_columns'. QODBC3: Unable to execute statement

Maybe there a schema migration I need to run?


texts:
How to get coords for NW and SE if we have about 4-5 pixels on the image with assigned latit/longit? I'm working on javascript app which would overlay that image on Google Map then.

Thanks.


texts: Orienting line labels to graticules in QGIS

texts:
I have been trying to add labels to a map that are oriented parallel to the graticule lines
In ArcGIS, with Maplex, it is simple link, but in QGIS I cant find any way to do this.


texts: Vector Building Footprint Resources - Tokyo

texts:
I typically know where to look for this data in the United States, but I'm at a loss when looking for building footprint shapefiles abroad. I thought I could find what I was looking for with Open Street Map data, but when I check out http://downloads.cloudmade.com/asia/eastern_asia/japan/tokyo#downloads_breadcrumbs I don't get the files I need.
Anyone know other resources, either official or osm?


texts: Resizing symbols in QGIS Print Composer legend?

texts:
I am trying to map a count data type via graduated symbol against a chloropleth and the legend is too big.

The circles of the graduated symbols are the actual size of their corresponding values in the map, which is good for interpreting the legend but is not good because they overlap and in order for space to be created between them, the legend stretches absurdly along the y axis.

Is there any way to scale the circles back to half the size of the mapped symbols?

I stretched out the legend as you can see below, aesthetically it wasn't ideal.

got the answer!


texts: Add PostGIS layer using PyQGIS standalone application

texts:
I'm trying to add a PostGIS layer to my canvas, but when I debug it I see that QgsVectorLayer does not return anything, however that same code was copied from a plugin that I am making and if the connection and creation of the layer works well. It seems that for this type of applications outside of QGIS it should be done in another way but I do not know it.
This is the code I am using
    def cargarPostGIS2(self):
        uri = QgsDataSourceUri()
        uri.setConnection("localhost", "5432", "catastro", "postgres", "postgres")
        uri.setDataSource("catastro", 'municipios', 'mi_shape')
        vlayer = QgsVectorLayer(uri.uri(False), 'municipios', "postgres")
        if vlayer.isValid():
            # Agregar el layer al registro
            QgsProject.instance().addMapLayer(vlayer)
            self.canvas.setLayers([vlayer])

How resolve this problems?


texts: Grade separation routing using pgRouting

texts:
I am using pgRouting version 2.6.
I am creating a routing topology from a subset of lion data https://www1.nyc.gov/site/planning/data-maps/open-data/dwn-lion.page
I have gone through the process of setting up the network and I am able to do some basic shortest path routing using the pgr_dijkstra function.
However I am running into the routing issue of overpasses/bridges/tunnels where the network is not aware of the grade separations or heights of the roads that are intersecting. As a few resources on the web have mentioned this issue

How to create a routable road network from Ordnance Survey data?
Create a network topology for pgrouting with grade separation
https://github.com/mixedbredie/highways-for-pgrouting

These questions and answers have helped me contextualize my problem but the solutions they came up with do no really apply to the dataset I am using.
I have already figured out the logic on what I need to do but I do not know how to implement it into either the network topology creation or the routing functions
Example
WITH
dijkstra AS (
    SELECT  *
      FROM pgr_dijkstra('SELECT id, source, target, cost_drive AS cost, rcost_drive as reverse_cost
                         FROM edges where driveable = true',
           getnearestnode(40.5993439,-74.0642904),
           getnearestnode(40.71204,-73.9616284),
           --getnearestnode(40.7059044,-74.0078728),
           directed := TRUE)
           )
    SELECT *
    FROM dijkstra left JOIN edges ON (edge = id)
    ORDER BY seq

The query above routes from Staten Island to Williamsburg BK and for the most part it does a good job accurately finding the route however in the picture below the route is on the BQE and then makes an impossible left turn onto Borinquen Place.

I have been able to identify when a turn or route should not be allowed.
In the data there are two columns that deal with gradient. nodelevelf (Level code indicator vertical topology at the start of the street segment. goes from A-Z where z is the highest level) and nodelevelt (Level code indicator vertical topology at the end of the street segment. goes from A-Z where z is the highest level)
Two segments should ONLY be routed together if the nodelevelt of the 1st segment equals the nodelevelf of the 2nd segment.
In the above example the BQE segment (highlighted in red) has a gradient of I for both its nodelevelf and nodelevelt and the next segment which is a left turn onto Borinquen Place has a nodelevelf and nodelevelt of M which is a level higher than the previous BQE segment. This should not happen.
How do I implement such a rule in the routing?
Things I have tried:

pre-processing the data to merge streets with similar attributes but this creates multilinestrings which does not work in the pgrouting functions
run the above SQL query to generate the route then identify the wrong turn link and delete that node from the vertices table (abandoned that idea once other routes were not working because of it)


texts: Web AppBuilder import error: The theme used by the app is in a newer version

texts:
I am using the latest version of the Web AppBuilder SDK (v 2.19 is the latest as of Jan 2021) and I was trying to (unsuccessfully) import an app that was created via ArcGIS Online.  The app uses the dashboard theme and features a very basic web map and three infographic widgets.  I tried importing it via the account option and uploading a zipped up version of the app.  Both methods produced the following message in the SDK logs:
[2021-05-07T15:32:52.965] [ERROR] app - fail to import app. { message: 'Fail to import. The theme used by the app is in a newer version. Name:[HD Test Map]' }
[2021-05-07T15:34:01.535] [ERROR] app - { message: 'Fail to import. The theme used by the app is in a newer version. Name:[HD_Test_Map]' }
To the naked eye, I could not see any discernible differences between the dashboard themes in the SDK and AGO.  Is this a bug?


texts: Appending label to every layer in current map bounds using Leaflet

texts:
I am trying to reduce the amount of permanent tooltip on my leaflet map for performance reasons. I have over 5000 geoJSON Polygons in my non-geographic map. At a specific zoom level, my Polygons should receive a label a.k.a Tooltip.

So my idea is to check each layer if it is in the current map bounds. But getting bounds and center LatLngs of over 5000 layers and checking if that layer is in the map bounds is time-consuming. If I add this in the map.moveend Event my map freeze. In addition, I have to fetch text for my labels from an SQLite Database.

Is there another solution for this problem?

EDIT:

My code so far:

var allLayer = L.geoJson(xhr.response, {style: style}).addTo(mymap);

var layerInCurrrentMapBounds = [];
mymap.on('moveend', function(e) {
    if (mymap.getZoom() >= 5) {
        allLayer.eachLayer(function(layer) {
            var bounds = layer.getBounds();
            if (bounds.isValid()) {
                var latLng = bounds.getCenter();
                if (mymap.getBounds().contains(latLng)) {
                    layerInCurrrentMapBounds.push(layer);
                }
            }
        });

        layerInCurrrentMapBounds.eachLayer(function(layer) {
            layer.bindTooltip("test", {
                permanent: true,
                className: "my-label",
                offset: [0, 0]
            });
        });
    }
});

On the next moveend I have to remove all tooltips. In addition, I have to fetch my tooltip text from a sqlite database. And what if the user starts to drag the map rapidly. Could I make my function async and cancel this every time the user starts to drag again?


texts: How can I show custom titles on my table fields with OpenLayers?

texts:
I have a server with Geoserver, Apache and OpenLayers and it displays the map okay; I can see my data with an OpenLayers.Popup.FramedCloud call and it shows the content of my layer.

I want to know if it's possible to show custom titles for my data instead of the ones in the  table.

Example:

I have COSTRUC field and I want to show Construction

So there is any chance to do it with Openlayers?


texts: ZIP and import to ArcGIS online a personal geodatabase

texts:
This may sound dumb, but how do you zip a personal geodatabase in order to import it to ArcGIS online?

Are there other ways to import the geodatabase, different from the one below?

ArcGIS Online fails to add zipped File Geodatabase?

There is the "share as" option in ArcMap 10.1 but it's been 2 hours since it's packaging the map... it's a small and simple map.


texts: Listing layers with or without geometries using PyQGIS

texts:
I need to list all layers with geometries but not those without.
Can I be sure geometryType = 4 always are tables without geometries?
I cannot find the documentation explaining which number is which type in plain text.
for lyr in QgsProject.instance().mapLayers().values():
   #if lyr.geometryType() != 4: #Do something
   print(lyr, lyr.geometryType())


texts:
I am using ArcPad 10.2 to do some field data collection.  On of the functions I would like to be able to do is associate a photo (taken in the field with the built in camera) to an existing point layer.  My assumptions is that I will need an empty field to store the photo, but I am unsure of what that field needs to be as a data type.

Any help and info would be appreciated.

Thanks...


texts: How do I use polygon centroids within a GeoQuerySet?

texts:
I have two models: neighbourhoods and boroughs, each containing a polygon field.

I am trying to construct a QuerySet that will give me the neighbourhoods within a given borough. I expected to use within or contained:

neighbourhoods_in_borough = Neighbourhood.objects.filter(poly__within=Borough.poly)

However, it is not always the case that a neighbourhood lies cleanly within a borough. Neighbourhood polygons roughly align with a Borough, but may weave in and out slightly. As such, filters like within or intersects either leave out neighbourhoods or include too many surrounding a Borough which happen to slightly overlap.

Because I only want the Neighbourhoods which are predominantly inside of a given Borough, I thought I would test to find Neighbourhoods whose centroids are within a given Borough:

neighbourhoods_in_borough = Neighbourhood.objects.filter(centroid__within=Borough.poly)

However, you can spot the issue here: centroid is not a field on the model, and so the above does not work. I did think adding a preceding query (e.g. Neighbourhood.objects.centroid() ) would make the centroid available to be used in some form, but it fails for the same reason as the example above.

Is it possible to filter based on a polygon's centroid rather than the polygon's shape?


texts:
I have a gpkg with multiple tables, and I am trying to load it with geopandas:

import geopandas as gpd
filename = '~/example.gpkg'
gdf = gpd.read_file(filename)'

This returns the first table/layer in the gpkg. Is there any way to access the other layer, i.e. something like:

filename = '~/example.gpkg|second_table'

I'm guessing I might need to use fiona directly, but couldn't figure it out from the docs


texts: OpenLayers - How to show UK Map with OSGB Projection

texts:
We are using  OpenLayers for our GIS Application.

We use Google and BingVE as source for the BaseMap.  We use Projection in the OpenLayers with following code. We want to show UK Map with OSGB.

if(source=='bingve')
{
    zoomlength=17-minZoomLevel;
    map=new OpenLayers.Map('map',{

        projection:new OpenLayers.Projection('EPSG:900913'),
        displayProjection:new OpenLayers.Projection('EPSG:27700'),
        units:'m'

    });
}

Does this will Automatically show Map in UK format OR we need to explicitly call Transform Method to convert Lag / Lon.


texts:
I'm using Qgis 3.0.1 and trying to polygonize a raster layer. I'm getting the following error:

  Processing algorithm…
  AlgorithmDialog
  Input parameters:
  { 'INPUT' : 'crs=EPSG:4326&dpiMode=7&format=image/png&layers=klim-stat_temperatur:medeltemp_maj_yta&styles=&url=http://...', 'BAND' : 1, 'FIELD' : 'DN', 'EIGHT_CONNECTEDNESS' : False, 'OUTPUT' : 'C:/...' }

  GDAL command:
  cmd.exe /C gdal_polygonize.bat crs=EPSG:4326&dpiMode=7&format=image/png&layers=klim-stat_temperatur:medeltemp_maj_yta&styles=&url=http://... "C:..." -b 1 -f "ESRI Shapefile" None DN
  GDAL command output:

  gdal_polygonize [-8] [-nomask] [-mask filename] raster_file [-b band|mask]

  [-q] [-f ogr_format] out_file [layer] [fieldname]

  'dpiMode' is not recognized as an internal or external command,

  operable program or batch file.

  Invalid parameter - /png

  'layers' is not recognized as an internal or external command,

  operable program or batch file.

  'styles' is not recognized as an internal or external command,

  operable program or batch file.

  'url' is not recognized as an internal or external command,

  operable program or batch file.

  Execution completed in 0.86 seconds
  Results:
  {'OUTPUT': }

  Loading resulting layers
  The following layers were not correctly generated.
  C:/...
  You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

Any ideas of what it could be?


texts: canvas.toBlob(function(blob)) is giving SecurityError: The operation is insecure

texts:
I'm trying to get the image of Map from OpenLayers 4.6.5 and the following error is appearing each and every time.

I'm using cdn link of OpenLayers 4.6.5 for the base Map in JavaScript with Apache Tomcat server in web application (HTTP format) everything is working fine.

But when I start to get the image from OpenLayers Map
it's giving :

  SecurityError: The operation is insecure.

in case of FireFox and

  Uncaught DOMException: Failed to execute 'toBlob' on 'HTMLCanvasElement': Tainted canvases may not be exported.

It's really head eating and I have tried many ways
like CORS, timeout and anything available from other sources:

   canvas.toBlob(function(blob)
            {

             var imgData = new Image();
             imgData.setAttribute('crossOrigin', 'anonymous');
             imgData.src = canvas.toDataURL('png');
           }

I found out something like (I'm not sure) : it's due to security of the network i.e I'm using secure cdn link from OpenLayers :

https://cdnjs.cloudflare.com/ajax/libs/openlayers/4.6.5/ol.js

which is https:// and my application is http:// ---- is this a problem? I'm really tired of this.

Could anyone look into it?


texts:
I have a couple virtual machines with Python (2.7.16 and 3.6.6 - 64 bit) with GDAL installed using GIS internals installers and using this process

However when running gdalinfo --version on the shell, in both I get an error:

C:\Users\K>gdalinfo --version
ERROR 1: Can't load requested DLL: C:\Program Files (x86)\GDAL\gdalplugins\ogr_MSSQLSpatial.dll
126: No se puede encontrar el m¾dulo especificado.

ERROR 1: Can't load requested DLL: C:\Program Files (x86)\GDAL\gdalplugins\ogr_MSSQLSpatial.dll
126: No se puede encontrar el m¾dulo especificado.

GDAL 2.2.3, released 2017/11/20

C:\Users\K>gdalinfo --version
ERROR 1: Can't load requested DLL: C:\Program Files\GDAL\gdalplugins\ogr_MSSQLSpatial.dll
126: No se puede encontrar el m¾dulo especificado.

ERROR 1: Can't load requested DLL: C:\Program Files\GDAL\gdalplugins\ogr_MSSQLSpatial.dll
126: No se puede encontrar el m¾dulo especificado.

GDAL 2.4.0, released 2018/12/14

The dll is located in the gdalplugins folder. Why does it fail in both installations?


texts:
I have created a 1 km buffer around a polygon using Geoprocessing Tools in order to identify local sites which lie within 1 km. However, whilst one local site is showing as being within 1 km of my polygon, when I measure the distance using the Measure Line tool, the closest distance that I make it is 1,064 m, so why is it included in the buffer?

Is it a smoothing error, is the Measure Line tool not as accurate?

I am currently working in EPSG:27700.


texts: Create and show a vector layer from Geoserver using Openlayers3

texts:
It seems the transition from OL2 to OL3 will be harder than I thought. There are significant differences and the online support its still limited (regarding examples etc.)

I am trying to do something relatively simply..to display a WFS layer on my map using OpenLayers3. I have found this example and tried to modify it to my needs.

He presents 3 different ways to create a vector layer from geoserver and I used the last one (JSONP format). That means that I had to enable JSONP in geoserver as described here.

When I try to see the map I get a blank page with the error:

Uncaught TypeError: undefined is not a function

on line 82 which is:

strategy: ol.loadingstrategy.tile(new ol.tilegrid.XYZ({

If I change this line into:

strategy: ol.loadingstrategy.box

then I see the map layer (but not the vector layer) and I get this error:

XMLHttpRequest cannot load http://localhost:8080/geoserver/wfs?service=WFS&version=1.1.0&request=GetFe…4558475.055843057%2C2678334.3618555046%2C4592871.718571385%2CEPSG%3A900913.

No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'null' is therefore not allowed access.

This is my code in the init() function

function init() {

var center = ol.proj.transform([23.709719,37.973333], 'EPSG:4326','EPSG:900913');
var view = new ol.View ({
    center: center,
    zoom: 12,
});

var newLayer = new ol.layer.Tile({
    source: new ol.source.OSM()
});

// WFS layer from Geoserver

loadFeatures = function(response) {
formatWFS = new ol.format.WFS(),
sourceVector.addFeatures(formatWFS.readFeatures(response));
};

sourceVector = new ol.source.Vector({
loader: function(extent) {
    $.ajax('http://localhost:8080/geoserver/wfs',{
        type: 'GET',
        data: {
            service: 'WFS',
            version: '1.1.0',
            request: 'GetFeature',
            typename: 'dSpatialAnalysis:categoriesdata',
            srsname: 'EPSG:900913',
            outputFormat: 'application/json',
            bbox: extent.join(',') + ',EPSG:900913'
            },
        }).done(loadFeatures);
    },
    strategy: ol.loadingstrategy.bbox
});

layerVector = new ol.layer.Vector({
    source: sourceVector,
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'rgba(0, 0, 255, 1.0)',
            width: 2
        })
    })
});

map = new ol.Map({
    target:'map',
    renderer:'canvas',
    view: view,
    layers: [newLayer, layerVector],
});

}

I have downloaded the library and call it locally:

<!-- Load Openlayers - openlayers 2.13 cause issue with geoserver-->
<script src="libs/OpenLayers-3/ol.js"></script>
<script src="https://code.jquery.com/jquery-1.11.2.min.js"></script> <!-- add jquery -->

EDIT

Removing the port from the URL makes the request correctly and I can see in Firebug the response. But I still have issues displaying the points on the map.


texts: ArcGIS javascript API Update Records of a Feature Layer

texts:
I'm new to ArcGIS javascript API, now I want to create a feature layer and then add an array of attributes or simply update the array of attributes of a specific field. I only need to replace the original values array as my own array. Is that possible?
There's a method: attr(name,value), Can I use this?
Can any one give me some hints or code? Thanks in advance!


texts:
I'm actually trying to calculate a grid point position in a SRTM grid by row and column count.
So I got a position by lat/lon, then calculating the row and column in my hgt-file (this post helped me) and reading the height. After that I want to calculate the position of the grid point at the specific row and column. For that I'm calculating the distance from the grid point to the origin of the grid and then I get the lat/lon of the grid point.

An example:

Lat: 47.743718 = 47º 44' 37.3848"
Lon: 09.260396 = 09º 15' 37.4256"

=> I have to use the file "N47E009.hgt"

47º 44' 37.3848" => (44 * 60)" + 37.3848" = 2677.3848"
09º 15' 37.4256" => (15 * 60)" + 37.4256" =  937.4256"

GridRow = floor(2677.3848" / 3") // =892
GridCol = floor(0937.4256" / 3") // =312

fileRow = 1201 - GridRow
fileCol = (1201 * (fileRow - 1)) + GridCol

//Reading height from hgt file

Calculating position of grid point relative to origin (47N09E):

x = GridCol * 3 * 30.87 //3 arc secs and 1" = 30.87m
z = GridRow * (3 * 30.87 * cos(47)

=> x = 28894m   z = 56339m

If I'm calculating the distance with lat/lon with a tool I'm getting this results:

=> x = 19470m   z = 82700m

So where is my error in reasoning?


texts: Loading Vector data into OpenLayers3 from Geoserver

texts:
I'm using Geoserver to serve my spatial data which is being displayed using OpenLayers 3.

I've so far been able to load my data as TilesWMS from my geoserver, but in order to select the data I need to be able to load it in as a Vector.

I'm using this code as my import

var vector = new ol.layer.Vector({
    source: new ol.source.GeoJSON({
        projection: 'EPSG:4326',
        url: 'http://localhost:8080/geoserver/LiDAR/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=LiDAR:lidar_projects_completed&maxFeatures=50&outputFormat=application/json'
    }),
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'rgba(255, 255, 255, 2.0)',
            width: 2
        })
    })
});

map.addLayer(vector);

Loading the URL into a web browser returns the multipart feature, so it leads me to think I need to handle it differently then just 'adding the layer'

Any suggestions?


texts:
I want to make a choropleth map from Census data, where my table structure is similar to the below, but I want to create a graduated style for the population based on a particular Country of birth (e.g. for Greece):

Region     Country of Birth         Population
1          Australia                52
1          Iraq                     22
1          Greece                   13
2          Australia                16
2          Iraq                     3
2          Greece                   21

I've tried doing this via an expression, however, have been unsuccessful so far.


texts:
I have two rasters: Raster A and B. Raster A has has only values of 0. The white area in Raster B are No Data values. I want to merge them to one raster (Raster C). So basically I want all pixels in Raster B with a No Data value to have a value of 0.

I tried to do this with the raster calculator in QGIS by adding Raster A to Raster B, but it only returns again Raster B.

I am looking for solutions with QGIS, GDAL, or Grass GIS.


texts: How can I programatically create a (tiff) raster file for a DEM to reference a constant altitude?

texts:
I'm trying to create a DEM to go with a simple shapefile (one polygon, four points) in python.

I used fiona to create the shapefile and now need to create the DEM as a tiff file to denote the elevation of the points in the shapefile. I'm stumped. Everything I google looks sooo complicated. How hard can this be?

Some things I don't understand yet include:

what is the relation of the dimensions (in pixels) of the tiff file to the shapefile extents? does it matter?
what kind of meta data is also stored in the tiff file? e.g. is the CRS stored in the tiff file?

I will be using the output with ArcScene Desktop 10.4.1, but would really like to create everything programatically. I guess if the solution involves arcpy or arcgisscripting, that is not a big problem, but geopandas and fiona sound like the way to go, right?


texts: Use QGIS 3.2 for TMS such as leaflet 'TMS:true' tag

texts:
Forgive me for using a loose title, I wanted to add a XYZ tile (http://this/is/url/{z}/{x}/{y}.png) with QGIS 3.2 XYZ tiles. Unfortunately, I had some issues showed in the image.

The TMS layer works queit well with Leaflet.js with 'tms: true'

var lyr = L.tileLayer('http://this/is/url/{z}/{x}/{y}.png', {tms: true, opacity: 0.8, attribution: ""});

You can see that link for more imformation about the TMS Layer.


texts:
I have 80 LANDSAT 8 Images. I can able to create composite image using Arcgis and then Pansharpen it to produce 15m multi-spectral image as a input for classification process. My requirement is to extract only roads vector(Polygon) data from LANDSAT 8 Images using ecognition developer.
After segmentation, I have tried NDVI,NDWI ruleset in eCognition to unclassify vegetation and water areas(so that I can get road) but it doesn't seem possible to classify road data completely.

Is there any specific ruleset available for extracting roads?

I have licensed Arcgis 10.2.1, Ecognition developer 9 as of now.


texts: Building foot print detection from raster imagery in QGIS?

texts:
I read the article about Building RECognition tool (BREC) from  GEM foundation visit www.globalquakemodel.org/resources/use-and-share/tools-apps/

Tried to work on Building foot print extraction in QGIS with this plugin BREC4GEM visit http://www.ticinumaerospace.com/foss/

But, I was unable to load the plugin in QGIS 2.14 Essen, tried in Mac OS X EI Capitain, Windows 8 and OSGeo Live 10.  using this procedure. plugin via `Settings->Options->System->Plugin paths->Add->"folder path/BREC"'.  I solved my problem by going to the folder of plugins ~/.qgis2/python/pluginsand pasting the folder 'BREC' in it, in mac OS X & Windows 8.

Did anyone work on this plugin "BREC4GEM" before OR there are other plugins to extract building foot print from raster imagery ?

I also read the 'User guide: Geospatial tools for building footprint and homogeneous zone extraction from imagery' and want to implement this for Sentinel 2 image data set.   visit https://www.globalquakemodel.org/resources/publications/technical-reports/user-guide-geospatial-tools-building-footprint-and/

Can anyone suggest other raster imagery from Open Source projects for sole purpose of building foot print extraction ?


texts: Statistical validity percentage check on raster

texts:
I would like to check if my set of NDVI images with min/max values from -1 +1 is good for at least 80%.
I took my cue from this other topic: here. Considering that the nodata is -32768
But the values I get on multiple images are all the same.
with rasterio.open(NDVI[0], 'r') as src:
    meta = src.meta.copy()
    meta.update({"nodata":-32768})
    data = src.read(1)
    arr=data[data>-1.0]
pct_valid = 100 * (arr != meta['nodata']).sum() / (meta['width'] * meta['height'])

print('%.4f'%pct_valid)

OUTPUT:
-6.87526

How should I interpret it for my image?
To select only those images that are 80% good instead?
if pct_valid > 80.0:
    write(raster)


texts:
I am having problems opening the OpenLayers plugin (any of them).  It seems to be a Python error. I have searched around in previous questions to find something, but none of the solutions I came across fixed the issue.

I have posted the error message below:

An error has occurred while executing Python code:

QgsCsException: forward transform of (-140537588568161501184.000000, -99802055649825423360.000000) PROJ.4:  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 +to  +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs Error: latitude or longitude exceeded limits
Traceback (most recent call last):
  File "C:/Users/Russells/.qgis2/python/plugins\openlayers_plugin\weblayers\weblayer.py", line 100, in addLayer
    self._addLayerCallback(self)
  File "C:/Users/Russells/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py", line 158, in addLayer
    self.setMapCrs(coordRefSys)
  File "C:/Users/Russells/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py", line 206, in setMapCrs
    extMap = coordTrans.transform(extMap, QgsCoordinateTransform.ForwardTransform)
QgsCsException: forward transform of
(-140537588568161501184.000000, -99802055649825423360.000000)
PROJ.4:  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 +to  +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs
Error: latitude or longitude exceeded limits

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.18.0 Las Palmas, 0332f5a

Python Path:
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python
C:/Users/Russells/.qgis2/python
C:/Users/Russells/.qgis2/python/plugins
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\bin\python27.zip
C:\PROGRA~1\QGIS2~1.18\apps\Python27\DLLs
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\plat-win
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGIS2~1.18\bin
C:\PROGRA~1\QGIS2~1.18\apps\Python27
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\PIL
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/Russells/.qgis2//python
F:/School/History/Mapping/AtlantaCampaign


texts: Non-destructive improvement of vector layer (shapefile) representing a road network?

texts:
I have a QGIS project with some vector layers backed by shapefiles. One of these layers contains all the roads, and I would like to "improve" it:

It is organized in a way that each vector feature goes between consecutive intersections, but I would like to join every same-name segment into one continuous feature;
The road network is incomplete (data is quite old), so I would like to add my own segments to it;
Occasionally there might be changed roads (unpaved roads that were "straightened" before being paved), that I would like to replace in the original vector layer.

I am not sure what workflow should I follow so that I take most advantage from the original layer, without modifying it, but being able to edit it - both automatically when convenient, or manually when necessary.

Some strategies I thought about:

Creating a new layer and copying the data, selectively, to another layer, from which I could remove or edit old data, and add new entities;
Creating a copy of that layer right away, losing reference with the original one, and editing it exclusively, much like a repository fork, "making it my own".

I ask because I am quite new to GIS in general and QGIS in special, and I am not aware of the good practices yet, so I am afraid of following the wrong approach and spending a lot of work that I would eventually have to discard and start over.


texts: Writing multiple layers to GeoPackage using writeOGR() in R?

texts:
I'm trying to write multiple layers to the same GeoPackage in R, but I get an error, Creation of output file failed. I have tried to search for documentation of reading and writing to .gpkg files with RGDAL, e.g. to figure out whether writeOGR() actually supports multiple layers, with little success. Is this even possible, if so, how to do it? Minimal working example:

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

data(wrld_simpl)

norway <- wrld_simpl[wrld_simpl$NAME == "Norway", ]
sweden <- wrld_simpl[wrld_simpl$NAME == "Sweden", ]

file <- tempfile("scandinavia", fileext = c(".gpkg"))

writeOGR(norway, dsn = file, layer = "norway", driver = "GPKG")
writeOGR(sweden, dsn = file, layer = "sweden", driver = "GPKG")

ogrListLayers(file)

There is apparently an ogr2ogr shell command that does the trick (hat tip mdsumner), which I can wrap in an R function. However, it would be neat if writeOGR() and/or st_write() in the sf package had this built in. I reckon it depends on GDAL's layer_options, but there doesn't seem to be an append-type option for GPKG in GDAL.

I could write a simple wrapper function for st_write() but native support in sf or rgdal would be better.


texts: how to read/write XML file using Openlayers?

texts:
I am new to OpenLayers. I am trying to load data from a datasource and show it on a OSM map, then also add the functions such as tagging and making new features as well as editing their attributes. My questions are:

what's the best (efficient in terms of speed) solution for data storage? should i store data in DBMS such as MySQL, or how is to store points in XML (GML) file and read and write the file?? How many users can load and save XML files at the same time?
If I want to do it with XML, how can I do it? any good posts that explain this?


texts: Heatmap Processing Exception in QGIS 3.2

texts:
I am working through training exercises with the Alaska training data using QGIS 3.2 on Windows 8.  To create a Heatmap, I load the airports shapefile.

I cannot find the Heatmap plugin using the Plugins Manager for QGIS 3.2 so I try using Processing Tools/Interpolation/Heatmap(KDE) and set the radius to 300000.  I get the following error:

Traceback (most recent call last):
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\Heatmap.py", line 215, in processAlgorithm
self.tr('Could not create destination layer'))
_core.QgsProcessingException: Could not create destination layer

Execution failed after 0.09 seconds

I am assuming the Interpolation/Heatmap(KDE) is not the correct tool to use.  Is there a way to get the Heatmap plugin?


texts: How can I programmatically display the "Create Features" dockable window in ArcMap 10?

texts:
When I start editing in ArcMap 10, by default the "Create Features" window is automatically opened and docked. Stopping the edit session hides this window.

When I use the IEditor:StartEditing method this does not display the "Create Features" window however if it is already displayed the IEditor:StopEditing method hides it.

How can I Start Editing programmatically and display the "Create Features" window?


texts: Cannot be processed raster statistics for polygon processing on the processing toolbox QGIS 2.18.1

texts:
Using the QGIS 2.18.1 application to process satellite image data and shapefiles.
Raster statistics for polygon processing on the processing toolbox with coordinate reference system Selected CRS (EPSG: 32651, WGS 84 / UTM zone 51N) experienced an error:

The following layers were not correctly generated.

Statistics

You can check the log messages to find more information about the execution of the algorithm


texts: Geoserver global services and geowebcahe

texts:
I have disabled the "global setting" in my Geoserver to require the workspace name in all URL endpoint.

Unfortunately, GeoWebCache does not work anymore.
It can return the tiles allready generated but cannot generate new tile with geoserver.

It say "400: Problem communicating with GeoServer"
For example :

http://localhost:8080/geoserver/adm_ws/gwc/service/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=adm3&[...]

Do you have any idea to solve this problem ?
Geoserver V2.13.1
GeoWebCache V1.13.1


texts: Metric system not working properly while creating buffer in QGIS

texts:
I am creating 10 meter buffer on line geometry shapefiles.
But the output buffer is too large.

CRS of my shapefile is EPSG:24383, which is in metric system.
I am creating buffer using Vector > Geoprocessing Tools > Buffer

The units of buffer are shown correctly in meters. But the output buffer created is too large. May be 100 or 1000 times large. Or maybe more than that, I have no way to find it out.

I have tried MMQGIS plug-in also, but no different result.

Is there anything which I am missing out here.


texts: Potree Converter - adding to output folder later - duplicate names issue

texts:
I use the Potree Converter (Windows application) to bulk convert .LAZ files and generate them into a (shared) output folder.  This works great, as each output file/folder generated is uniquely named.

If I want to convert more .LAZ files in the future and save them to the same output folder, some of the generated output files get allocated the names for files that already exist in the output folder.  I have checked the file sizes of the duplicate files and they are different - indicating they are indeed different files.  I do not want to overwrite any files unless they are the same file (i.e. contain the same data).

The only way around this issue I can see is to convert ALL the .LAZ files again - something that is obviously not practical - considering I have TBs worth of source .LAZ files.

Could someone tell me if it is possible to "add" to the output folder ensuring the latest files generated have names unique to the output folder?

As a very basic example, the shared directory contains the following files:

r.hrc
r.laz
r0.laz

If I want to add to this directory later by converting some latest source data, chances are the above files will be re-created but for the latest source data.

I do not wish to have to store the latest files in a new directory.


texts: How to preserve formulas in XLS cells when working in QGIS?

texts:
Opening layers directly from XLS spreadsheets is a very useful (if not well documented) feature of QGIS, but I can't seem to be able to properly import cells where the contents of the cell is the result of a formula.
The events are:

Set up an XLS spreadsheet where one column is the calculated result of others (e.g the sum of the other columns), the 'sum' column
Import spreadsheet as a table into Qgis by using 'add layer' and selecting the XLS file and clicking open.
The table loads up with all the correct values, including the calculated values in the 'sum' column
Save the table
Return to the XLS spreadsheet and your 'sum' column no longer has any formulas, it's just the integer values.

Is there anything I can do in the importing or saving of edits that will preserve the column as a formula, not just the results?


texts: How can I use QgsRasterDataProvider.write from Python in QGIS?

texts:
I would like to change a raster cell value from Python in QGIS using QgsRasterDataProvider.write method. It takes the following arguments  (from QGIS API):

virtual bool QgsRasterDataProvider::write   (
    void *      data,
    int     band,
    int     width,
    int     height,
    int     xOffset,
    int     yOffset
)

How can I construct the pointer (data) with the new cell value?


texts: Is it possible to run a search cursor on a multipart feature class?

texts:
I have a multipart feature class of cities.  I want to run a search cursor on a cities feature class where:

one city is selected and a select layer by location is run to select all parcels that are within a certain distance X;
Values from the selected parcels are then summed together;
move on to the next city record

Currently, I am  having trouble getting my select by location statement correct.  Instead of only selecting the parcels that are within distance X of one city, it's selecting all the parcels that are within distance X of all the cities in the feature class.

    arcpy.MakeFeatureLayer_management(city, c_layer)
    cCursor = arcpy.SearchCursor(c_layer)
    print "Created search cursor on c_layer"

arcpy.MakeFeatureLayer_management(parcel, p_layer)

for row in cCursor:
    geoName = row.getValue("GEO_NAME")
    print geoName
    geoID = row.getValue("GEO_ID")
    print geoID

    arcpy.SelectLayerByLocation_management(p_layer, "WITHIN_A_DISTANCE", c_layer, search_radius, "NEW_SELECTION")
    print "Completed select by location on p_layer"
    y = arcpy.GetCount_management(p_layer)
    print y


texts:
I have a table with id and date
Users are not really IT-literate so anything involving a join might be too complicated for them, but they love Excel. They can easily change their spreadsheets but doing it on QGIS may cause errors.
Is there any easy way of updating the table? Ideally a user would just import their CSV with ids and dates, and the corresponding ids would get the dates in the CSV


texts: Can QGIS print nested or overlapping polygons in order?

texts:
I use MapInfo-Discover and have calculated Catchment areas for analysis. MI won’t print overlapping polygons in one table in any set order; dividing the data into 100’s of tables isn’t a solution.

Can QGIS do this?


texts: Creating a vector layer from loaded layer from GeoPackage

texts:
Firstly, I load layers from gpkg (GeoPackage). Secondly, I need to create a vector layer from one loaded layer which is named stops.
Second point I cant do it, because I need to mark a column, that one is x field and the other is y field and assign EPSG code.
This part of code doesnt work and I don't know if it's possible to do this way (specifically uri variable):

gpkg_name = os.path.splitext(os.path.basename(path2))[0]
path_to_gpkg = os.path.join(os.path.dirname(path2), gpkg_name + '.gpkg')
conn = ogr.Open(path_to_gpkg)

for i in conn:
    layer=iface.addVectorLayer(path_to_gpkg + "|layername=" + i.GetName(), i.GetName(), 'ogr')
    print(i.GetName())
    if i.GetName()== 'stops':
        uri = (path_to_gpkg + "|layername=" + i.GetName() + ",&xField=stop_lon" + ",&yField=stop_lat" + ",&crs=epsg:4326").format(i)
        print(uri)
        name = i.GetName()
        #os.path.splitext(os.path.basename(i))[0]
        layer = QgsVectorLayer(uri, name, 'delimitedtext')
        print(layer)
        QgsProject.instance().addMapLayer(layer)

There is a whole code:
from qgis.core import *
from osgeo import ogr
import os.path
import shutil

path2=r'C:\Users\Martin\Desktop\gpkg' #rewrite path

if not os.path.exists(path2):
    os.mkdir(path2)

files = []
    # r=root, d=directories, f = files
for r, d, f in os.walk(path2):
        for file in f:
            exten = os.path.splitext(os.path.basename(file))[1]
            if exten == '.txt':
                files.append(os.path.join(r, file))

options = QgsVectorFileWriter.SaveVectorOptions()
firstt = True
for file in files:
    uri = 'file:///{}?delimiter=,'.format(file)
    name = os.path.splitext(os.path.basename(file))[0]
    layer = QgsVectorLayer(uri, name, 'delimitedtext')
    if firstt == True:
        options = QgsVectorFileWriter.SaveVectorOptions()
        options.driverName = 'GPKG'
        options.layerName = "_".join(layer.name().split(' '))
        error_message = QgsVectorFileWriter.writeAsVectorFormat(layer,path2,options)
        firstt=False
    else:
        options = QgsVectorFileWriter.SaveVectorOptions()
        options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
        options.driverName = 'GPKG'
        options.layerName = "_".join(layer.name().split(' '))
        #print(options.layerName)
        error_message = QgsVectorFileWriter.writeAsVectorFormat(layer,path2,options)
        if error_message:
            print(layer.name(), error_message)

gpkg_name = os.path.splitext(os.path.basename(path2))[0]
path_to_gpkg = os.path.join(os.path.dirname(path2), gpkg_name + '.gpkg')
conn = ogr.Open(path_to_gpkg)

for i in conn:
    layer=iface.addVectorLayer(path_to_gpkg + "|layername=" + i.GetName(), i.GetName(), 'ogr')
    print(i.GetName())
    if i.GetName()== 'stops':
        uri = (path_to_gpkg + "|layername=" + i.GetName() + ",&xField=stop_lon" + ",&yField=stop_lat" + ",&crs=epsg:4326").format(i)
        print(uri)
        name = i.GetName()
        #os.path.splitext(os.path.basename(i))[0]
        layer = QgsVectorLayer(uri, name, 'delimitedtext')
        print(layer)
        QgsProject.instance().addMapLayer(layer)


texts: Debugging AttributeError: 'module' object has no attribute 'Parameter' from Python Toolbox of ArcPy?

texts:
I have a big .pyt script (python toolbox) and now i trying split it to many files (1 file - 1 tool).

In single .pyt file everything works perfectly, but when file is splited i get this message:
AttributeError: 'module' object has no attribute 'Parameter'.

Structure of files:

My Catalog:
-- toolbox.pyt
-- toolpackage:
---- configurator.py
---- __init__.py

toolbox.pyt:

# This Python file uses the following encoding: utf-8

import arcpy

from toolpackage.configurator import ToolboxConfigurator

class Toolbox(object):
  def __init__(self):
    """Define the toolbox (the name of the toolbox is the name of the .pyt file)."""
    self.label = "label"
    self.alias = "Tools"

    # List of tool classes associated with this toolbox
    self.tools = [ToolboxConfigurator]

configurator.py:

# This Python file uses the following encoding: utf-8

import arcpy

class ToolboxConfigurator(object):

  def __init__(self):
    """Define the tool (tool name is the name of the class)."""
    self.label              = u"config"
    self.description        = u""
    self.canRunInBackground = False

  def getParameterInfo(self):
    """Define parameter definitions"""
    new_config_file   = arcpy.Parameter(
      displayName     = u"?",
      name            = "new_config_file",
      datatype        = "GPBoolean",
      parameterType   = "Optional",
      direction       = "Input")
    parameters = [new_config_file]
    return parameters

  def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True

  def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal validation
    is performed.  This method is called whenever a parameter has been changed."""
    return

  def updateMessages(self, parameters):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return

  def execute(self, parameters, messages):
    """The source code of the tool."""

__init__.py is clear.

My error:

Traceback (most recent call last):
  File "C:\tools_v4\toolpackage\configurator.py", line 15, in getParameterInfo
    new_config_file   = arcpy.Parameter(
AttributeError: 'module' object has no attribute 'Parameter'


texts:
I am trying to print a map that I have made in QGIS print Composer. When I try to print the map or export it to a PDF (or any other format) as a raster all I get is a blank page. When I try to print not as a raster (unchecking the "print as raster option") it will print but obviosly not the way I want it to. eg. Layers arent transparent nor do other bits look right. I don't understand why this would happen. Why wont do I have a blank page when trying to print as a raster?

I am running QGIS 2.0 on a UBUNTU 13.10/64bit/16GiB ram


texts: ArcGIS Conversion Error While Converting ED50 Datum to WGS84 Datum

texts:
With ArcGIS for Desktop, I tried a datum conversion of Transverse Mercator ED50 Datum to Transverse Mercator WGS84 Datum.

Normally, the difference between two projections is approximately 180 meters but when the conversion is completed with ArcGIS' "Project" tool, the shift of the two features becomes 80 meters, not 180 meters.

On a projected data frame, ArcGIS shows the exact same place for two datums, which is an expected issue, but in fact the WGS84 one is not at the correct location.

After completing this process with a different software I opened two features on the same data frame with an ED50 datum.

Features are not on the exact locations and there is approximately 100 meters visualization error (it means, the datum conversion error of the ArcGIS is 100 meters):

A data frame without a projection and datum shows the real difference between the datums, which is approximately 180 m:

Is this a known issue with ArcGIS or is there a solution for converting ED50 to WGS84 with "Project" tool of ArcGIS correctly?


texts: Polygons have holes after pgr_pointsAsPolygon

texts:
I calculate a 500m Isochrones with pgr_drivingDistance and after that use pgr_pointsAsPolygon to create the Isochrones out of the resulting nodes.

SELECT * FROM pgr_pointsAsPolygon('SELECT id,x,y FROM Calc_nodes')

Sadly in some case the Dijkstra algorithm create holes in the resulting poylgons. I know the reason for this is the nature of the network and how the algorithm works.

I just need function to check if there are holes and then get rid of them. Any suggestions?


texts:
i have a game that described a planet and is a sphere.

on the client side, the world is separated to tiles (images) and i mapped it follow by this and this, so now when the user is zooming / dragging the world he's see the right location follow by the middle-center geo coordinate and global pixel.

now, on my server side, i need to put data on my planet (in other words, to draw it on the tiles i mentioned above - but it's doesn't matter for this discussion).

i have array of data items with same dimensions (width/height) that should be populate on the planet.
from my own understanding, i need to find the geo coordinate of each data-item and then use the function geoToGlobalPX() to get it's location against the whole 2d planet's map, and finally 'put' it on the planet.

my problem is to find this geo coordinate, since i can't run over the planet from pixel(0, 0) to pixel(PLANET_SIZE, PLANET_SIZE) and each iteration grow the current-position in item-data-width because the range of location is decreasing as much you're going up/down from the equator, so on the top latitude (~85) i'll have more less data-items than on the equator.

i tried to read on the articles above about meters-per-pixel but i still don't really understood how i can integrate it on my system.

i hope you understood my and someone will have some advice for me.

thanks.

Shlomi.


texts: Interpolation of point data gives no data values in QGIS 2.18

texts:
I am for the moment, quite stumped. I have a layer of points to which i have retrieved values from an underlying raster. From these values, I am trying to create a TIN, using the Interpolation tool in QGIS. The tool runs with no errors, but the result is a raster of no data. I have tried:

Using QGIS 2.18 (Las Palmas de G.C:) and running Windows 10

Setting the regional settings to using a period and a comma for the decimal marker
Different tools, including IDW and squared TIN in the Interpolation dialog. This includes SAGA Trinagulation tool (this returns an error - tool not found)
Different projections for the points and the raster
Opening the points and running the analysis in a separate project
Changing the directories to avoid spacces in directory names
Changing the output file type (.asc and .tif)
Running the interpolation in an older version of QGIS (2.8)

Nothing has worked so far. The only diference was using squared TIN, that produced values of 0 instead of No data.

EDIT: here's a screencap showing the points and a sample of the attribute table. Second column is the data used, it's the same as the first with negative nubmers set to 0.

Any ideas?


texts: PostGIS Geography Buffer producing different areas?

texts:
I'm trying to create circles of equal area in different places, but depending on the location of the centroid different areas are calculated. For example

SELECT ST_Area(ST_Buffer('POINT(0 0)'::geography, 100000))/1000000

results in an area of 31,291 km², whereas

SELECT ST_Area(ST_Buffer('POINT(30 30)'::geography, 100000))/1000000

results in an area of 31,233 km².

Is this the correct approach to take for this problem, and if so, is there a method that would allow the creation of circles of equal area?


texts: Accessing MPS Atlas from Python with ArcPy?

texts:
I'm having trouble finding documentation on accessing MPS Atlas features using python in ArcGIS 10.

I need to be able to modify dataframe manager properties, swap map sheets, etc.

Is it part of arcpy.mapping?

From what I can see, there is no documentation of it on the ESRI Resource Center.


texts:
I have set up on Ubuntu a tile server by following this tutorial

I followed the instructions for downloading and configuring osm-bright, but now I want to use the openstreetmap-carto stylesheet.

In another question (mapnik rendering with OSM Carto style) I read that after importing osm data in the postgis database with osm2pgsql (with the use of the carto.style file), then I have to convert mml to xml via carto (carto project.mml > osm.xml)

For the osm-bright to work, before converting the project.mml to .xml, first I have to edit the configure.py file and then run the "preprocessor", e.g. run make.py.

Do I have to do something similar for openstreetmap-carto?

In the scripts folder there are 3 python scripts. Should I run them all, or just one them, or none of them?


texts: Representing json data on a map with GeoTools (Java)

texts:
I'm working on a Web Map Service in Java. Practically this WMS allows me to represent the map of open data type (shp, kml and gml). I'm going to start working even with data in json format, and for this I am using the GeoJSON libraries GeoTools. My intention is to convert a json file in shp format and parse it. I found several functions of https://gis.stackexchange.com/ but none of it works for me.

I'm using GeoTools 14.0.

Assuming that my WMS arrivals a request like this:

request = GetMap & Layers = hotel, school & Style = Red% 20Point, Black% 20Point

At this point, my WMS picks from this URL the things that I want to represent, namely school and hotel and especially how I want to represent them. Then I do a search for these layers within a directory which I agreed. Inside this folder, I check if a file by that name, and then do a check on its extension. For example:

public void parsing(String layers) {
    FileReader f=new FileReader();
    f.openFile("C:/Users/...../Servlet/Data");

    File[] lista = f.getFiles();

    for(File i:lista){
        String p=i.getAbsolutePath();
        if(i.isDirectory()){
            for(File fil:i.listFiles()){
                String g=fil.getAbsolutePath();
                checkFile(g,layers);
            }
        }else checkFile(p,layers);
        System.out.println("Layers: "+layers +" trovato");
    }                                                                                }

the checkfile function controls me the kind of extension:

public void checkFile(String p,String layers){
        if(p.contains(layers)){
            if(p.contains("json"))
                    try {
                        parsingJSON(p);
                        found=true;
                    } catch (IOException e) {
                        e.printStackTrace();
                    }
}

My intention is to make the even of open data parsing type "json" and I want to parser this dataset for example, but I should write a function that allows me to parser any json file . school.json

I thought about doing it this way, what do you think?

   public FeatureCollection parsingJSON(String p) throws URISyntaxException, SAXException, IOException {
        File geojson = new File(p);
        Map<String, Serializable> params = new HashMap<String, Serializable>();
        params.put("url", geojson.toURI().toURL());
        try {
            DataStore dataStore = DataStoreFinder.getDataStore(params);
            if (dataStore == null) {
            System.out.println("Error dataStore == null");
            }
            // read featureCollection
            Features handler = new Features();
            List<SimpleFeature> features = handler.getFeatures();

            collection = new DefaultFeatureCollection();
            for(SimpleFeature s:features){
                collection.add(s);
            }
        } catch (IOException e) {
        e.printStackTrace();
        }
        return collection;
    }

Unfortunately the line "DataStore dataStore = DataStoreFinder.getDataStore(params);" I generate a terminal error.


texts: Showing map inside polygon in layout QGIS

texts:
I have layout 1 with four different map presentations. Each of these
four has a name (ex: mapa_dv, mapa_dr, mapa_lf, mapa_rh, ...).

I have layout 2, where I want show "map_dv" inside an irregular polygon (Add Polygon can be irregular polygon).

How do I draw the "map_dv" of layout 1 within this irregular polygon of layout 2?


texts: Is is possible to use arcpy statement for multiple fields? e.s. arcpy.Intersect_analysis()

texts:
I am new in both coding and GIS. Now I am trying something new for me, which is using the for loop to execute some of the arcpy statement. For example, the Intersect_analysis statement.

Following is the code.

------------------------------------------------------
'''
interpretation of the Intersect_analysis statement:
Identity_analysis(in_features, identity_features, out_feature_class, {join_attributes}, {cluster_tolerance}, {relationship})

the one I copy from Arcmap when I did the statement using Arctoolbox:
arcpy.Intersect_analysis("HAWKER_CENTRE #","E:/university/NUS/SEM 2/courses/UD5628 - Sustainable Urban Design and Development/GIS/HAWKER_CENTRE_I.shp","ALL","#","POINT")
'''

# Import and setting worksapce
from arcpy import env

env.workspace = "E:\university\NUS\SEM 2\courses\UD5628 - Sustainable Urban Design and Development\GIS"

#List of the names of shp files to execute
fields = ["SUPERMARKET", "WET MARKET", "SHOPPING MALL", "HAWKER_CENTRE", "FOOD CENTRE", "CONVENIENCE STORE"]

#creating for loop to execute the statement for all the shp files

for i in fields:
    arcpy.Intersect_analysis(i "#", "E:/university/NUS/SEM 2/courses/UD5628 - Sustainable Urban Design and Development/GIS/i_I.shp", "ALL","#","POINT" )


texts: Problem using ArcGIS API for Silverlight

texts: I've install Visual Studio 2010 SP1, Microsoft Expression Blend 4, and ArcGIS API 3.0 for Silverlight. When I create a new project using a standard map application in Visual Studio, I see the error in design tab. When I create a project in Microsoft blend 4, I don't see any template and the appliction shows an error. How can I fix the problem?


texts: Stretched raster symbology not working in ArcGIS for Desktop?

texts:
I have some DEM rasters (from the same data source so I think they're uniform in all file specifications) that I'm symbolizing using the same stretched color ramp.  I set the stretch type to min-max, then edit high/low values to cover the low - high values for these rasters (0 - 128 in this case).  For most of the rasters this works perfectly to create a nice seamless map, however for 3 of the rasters it is still not using the same symbology resulting in big ugly seams.

Any ideas on what may be causing this?

I'm using ArcMap 10.2.2.


texts: Splitting shapefile per feature in Python using GDAL?

texts:
is it possible to split a shapefile per feature in python? (best would be a solution where i can temporarily save the resulting vector objects to the memory instead to disk).

The reason: I want to use the gdal rasterizeLayer function with several different subsets of shapefile. The function requires a osgeo.ogr.Layer object.

mkay, i tried a little bit around and it could work as following. You can get the geometry of gdal layer objects per feature as following.

#  Load shape into gdal
shapefile=str(vectorPath)
layer_source = ogr.Open(shapefile)
lyr = layer_source.GetLayer(0)
for i in range(0,lyr.GetFeatureCount()):
     feat = lyr.GetFeature(i)
     ge = feat.geometry()

Now i just need to know how to create a osgeo.ogr.layer object based on this geometry.

For clarification. I need a function in plain ogr/gdal code! This appears to be of some interest to other people as well and i still want a solution without any secondary modules (although any solution coming from here will be used in a free available qgis plugin).


texts: PyScripter crashes at "import processing"

texts:
I set up PyScripter latest version, 2.7) to work with QGIS (2.2, OSGeo4W install) following http://chelahmy.blogspot.ca/2011/02/setting-up-pyscripter-for-quantum-gis.html. It generally works well. I'm running Windows 7.

Here's my code:

import sys
import os
sys.path.append("C:\\OSGeo4W\\apps\\qgis\\python")
sys.path.append("C:\\OSGeo4W\\apps\\qgis\\python\\plugins")
from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt4 import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *

from processing.core.Processing import Processing

I need the sys.path.append to import qgis properly (rather than setting them in the PYTHONPATH). However, when I run the line

from processing.core.Processing import Processing

PyScripter simply shuts down, no error message, window just closes. Has anybody else experienced this? Same thing happens with other imports, like

import processing

I can import all the other plugin libraries (GdalTools, fTools, db_manager) with no problems - but I don't need those, I need the proceesing library to access all the nice OTB tools.


texts:
I have a series of contour lines with relative elevations (Z value). I am trying to align the labels so that they are all on the same imaginary line.
It's a process that I usually do with Autocad manually. Could someone tell me if there is an easier way in QGIS? The only result I can obtain is this:


texts: Create Python script to update field based on change in another field ArcGIS

texts:
I am trying to update a field "TPN" based on another field "PASS_NUM". In my script "TPN" is a newly added number field and thus blank, "PASS_NUM" is also a number field but it contains the wrong numbers. I want to assign values to "TPN" when "PASS_NUM" is no longer equal the the value in the row above it. The feature class table is sorted by "PASS_NUM". See example of attribute table I want to achieve here:

Is there any way to check if the value in the row above is equal to the next value (some sort of onchange event)?

The PASS_NUM field does not increment by 1 each time (i.e. it may skip 7 as in the example) but I want to TPN to only go up one.

Is this possible using the UpdateCursor in arcpy?


texts: How to color symbols by data in MapInfo?

texts: In MapInfo 10.0, how can I quickly colorize point symbols according to a field in the data table, without creating a new map? My data is just longitude, latitude and time, and I'd like to colorize the points according to the time field. This seems really basic but I can't find any way to do it.


texts: image.getDownloadURL download pixel grid dimension limit

texts:
I am trying to use image.getDownloadURL(downloadParam) to get the download link of an area of interest, it works fine with small areas, but when I pick a larger area, I will encounter an error like this: Pixel grid dimensions (48569x17664) must be less than or equal to 10000. Is the 10000 pixel dimension limit set by GEE on the server side? Is there any way to increase this limit? I am aware of the Export.image.toDrive function, to which you can set the maxPixels, I am using image.getDownloadURL(downloadParam) because I am building an EE app that allows users download raster data for a large area.
Thanks in advance for any helpful suggestions!
My code example is here: https://code.earthengine.google.com/84de75665f99d0cb78136bfe8045bfd2
//import SRTM data
var image = ee.Image("USGS/SRTMGL1_003");

//Define small and large AOI
var Small_AOI =
    ee.Geometry.Polygon(
        [[[-83.36097825867773, 36.44238302457335],
          [-83.36097825867773, 36.17346476392684],
          [-82.68806566102148, 36.17346476392684],
          [-82.68806566102148, 36.44238302457335]]], null, false);

var Large_AOI =
    ee.Geometry.Polygon(
        [[[-85.40031541688086, 37.58251126010928],
          [-85.40031541688086, 33.88278318496734],
          [-79.02824510438086, 33.88278318496734],
          [-79.02824510438086, 37.58251126010928]]], null, false);
//print(Map.getBounds())
Map.centerObject(Small_AOI)
Map.addLayer(Small_AOI)

//get download link for Small_AOI and Large_AOI
var url_small = image.getDownloadURL({
  image: image,
  //crs: ee.String(crs_textbox.getValue()),
  scale: image.projection().nominalScale(),
  region: Small_AOI,
  //maxPixels: 2e10,
})
print(url_small)

var url_large = image.getDownloadURL({
  image: image,
  //crs: ee.String(crs_textbox.getValue()),
  scale: image.projection().nominalScale(),
  region: Large_AOI,
  //maxPixels: 2e10,
})
print(url_large)


texts: Image disappeared after changing CRS in QGIS?

texts:
I am using QGIS 2.10.1
Why my image disappeared after changing the CRS?

My procedure:
Project>Project properties>CRS
Enable 'on the fly' CRS transformation

After I choose the desired CRS, my image is disappeared.


texts: Combine several raster files with heterogeneous bands

texts:
I am attempting to (a) combine several non-overlapping pop density raster files, (b) create a velox object from this raster, (c) and then extract the mean population value for each overlapping polygon from a shapefile I provide that covers the big raster.

I'm stuck on (a) with the error (R):

  Warning 6: gdalbuildvrt does not support heterogeneous band characteristics.

I don't know enough about raster files to know whether I need to somehow standardize the input rasters or pass an additional parameter to one of these functions.

library(gdalUtils)
library(rgdal)

# rasters
# 42: https://data.humdata.org/dataset/southasia_as42-high-resolution-population-density-maps
# 43: https://data.humdata.org/dataset/southasia_as43-high-resolution-population-density-maps
# 44: https://data.humdata.org/dataset/southasia_as434-high-resolution-population-density-maps
# 47: https://data.humdata.org/dataset/southasia_as47-high-resolution-population-density-maps

# zip with all, ~600MB: https://www.dropbox.com/s/26tc77upgpr51sr/rasterExample.zip?dl=0

# example https://stackoverflow.com/a/39554006/841405
  all_my_rasts <- list.files(pattern="\\.tif$")
  e <- extent(61.87125, 97.73931, 5.91125, 37.27042)
  template <- raster(e)
  projection(template) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
  writeRaster(template, file="MyBigNastyRasty.tif",
              format="GTiff",
              overwrite=TRUE)
  mosaic_rasters(gdalfile=all_my_rasts,
                 dst_dataset="MyBigNastyRasty.tif",
                 of="GTiff")
  gdalinfo("MyBigNastyRasty.tif")

Update 1:

I verified that all rasters have type "FLT8S". Removing rest of old update.

Update 2:

Building on @Jeffrey_Evans answer below, here's my first attempt at breaking up the processing. Basically, extract 2 rasters and rbind, then rbind that result to the third raster, and then to the forth, in the process removing things we don't need.

Here's where I stand at the end of step 1.

Is there a better way to approach this chunking?

# step 1
  e1 <- exact_extract(raster(r[1]), p)
  e2 <- exact_extract(raster(r[2]), p)
  results12 <- list()
  for(i in 1:nrow(p)) {
    results12[[i]] <- na.omit(rbind(e1[[i]], e2[[i]]))
  }
  remove(e1)
  remove(e2)
  gc()

# step 2
  e3 <- exact_extract(raster(r[3]), p)
  results123 <- list()
  for(i in 1:nrow(p)) {
    results123[[i]] <- na.omit(rbind(results12[[i]], e3[[i]]))
  }
  remove(e3)
  remove(results12)
  gc()

# step 3
  e4 <- exact_extract(raster(r[4]), p)
  results1234 <- list()
  for(i in 1:nrow(p)) {
    results1234[[i]] <- na.omit(rbind(results123[[i]], e4[[i]]))
  }
  remove(e4)
  remove(results123)
  gc()

  lapply(results1234, FUN=function(x) mean(x[,1]))
  p$rmean <- unlist(lapply(results1234, FUN=function(x) mean(x[,1])))


texts: Decimal Degrees on bottom right do not change

texts:
I added a shapefile into a new ArcMap 10.4.1 .mxd to check the coordinate system. It says Canada_Lambert_Conformal_Conic as the coordinate system.
When I do the mouse-over, the decimal degrees in the bottom right do not change. It just stays at this number:

How come the decimal degrees do not change when I move the mouse?


texts: Weighted mean of raster by polygons in R?

texts:
I have two rasters, raster_a and raster_b, and a shapefile, world_shapefile. Both rasters are of the world extent and overlay each other.

What I want to do is to create a weighted average of raster_a, where the weights are from raster_b, by each polygon in my world_shapefile.

I guess it is quite straightforward and goes something like: read in shapefile and rasters, create the weighed mean by polygons.

EDIT:
My code currently looks as follows:

require(raster)
require(sp)
require(gdal)

setwd("C:/Users/XXX/Desktop/data"

polys <- readOGR(dsn=get_wd(), layer = "world_shapefile"
raster.a <- raster(raster_a.txt)
raster.b <- raster(raster_b.txt)

plot(polys)
plot(raster.a)
plot(raster.b)

Everything looks perfect up to this point and then I get stuck...


texts: How to deactivate my custom tool by pressing the Escape key using PyQGIS?

texts:
I'm developing a QGIS plugin and I have created a custom tool. Once the process is over, if the user presses the ESCAPE key, I want my custom tool not to trigger mouse events on the map anymore.

How can I do that using python programming?

I have tried:

event.ignore()


texts: R: Clustering observations together based on distance and time

texts:
In R, I have a dataframe with roughly 3 million observations, with the columns being longitude, latitude and time respectively. My goal is to form clusters (using a custom distance function), and then form a single dataframe containing the observation from each cluster with the earliest time value. Two observations should be in the same cluster if the absolute difference in their time values is less than 100, and their spatial distance is less than 4km.
Firstly, here is some sample data:
library("geosphere") # contains the distm function, which we'll use to calculate distance given long-lat.
# Here is a reproducible dataset.
Long1<-c(149.1250,149.0774,149.1250,148.5352,149.0994, 149.0301, 149.0883, 149.0663, 149.1097)
Lat1<-c(-21.9876, -21.9171, -21.9876, -22.4645, -21.9206, -22.0115, -21.9188, -21.9249, -21.9546)
Time<-c(1000,500,754,129,6050,1908,6109,245,6049)
# Create our dataset
DatPoints<-cbind(Long1,Lat1,Time)

# Create a spatial distance matrix
DataMat<-as.matrix(cbind(Long1,Lat1))
Dist_Mat<-distm(DataMat,DataMat,fun=distHaversine)/1000 # creates a spatial distance matrix in km

From the distance matrix above, we can see that observation 5 and 7 have a spatial distance (1.1637km) that is less than 4km, and their time values differ by less than 100. So, observations 5 and 7 should be in the same cluster. Further, observations 5 and 9 have a spatial distance (3.931458) less than 4km, and their time values differ by less than 100, and so observations 5 and 9 should also be in the same cluster. So observations 5, 7, and 9 should all be in the same cluster, despite observations 7 and 9 being further than 4km apart.
I am wary of using any raster techniques because some of the observations (in my original file as well as the example data above) have the same lat-long coordinates. I have seen many similar examples on stack exchange, (such as Clustering spatial data in R?) and they all seem to suggest the use of the hclust function, using a "single" method. Indeed, if I was only interested in clustering using the spatial distance (with 4km as the cut off distance), then I think I can just do something like
chc<-hclust(as.dist(Dist_Mat,diag=TRUE),method="single")
chc.d40<-cutree(chc,h=4)

But I am unable to generalize this approach to the case where we need to consider a third attribute such as time (as opposed to just lat-long).
How can I achieve the goals that I specified at the beginning of this question?


texts: QGIS Python - Error: Feature 0 for attribute update not found

texts:
I'm new to QGIS python development. I'm trying to create new features in an existing vector layer and then to display the Feature Form so that the user can fill the attribute values.

The problem is, after the code is executed (after I fill the features atribute values and click Ok), an error is shown: Feature 0 for attribute update not found.

After that, I run again the code to create another feature with different attribute values. No error appears, but I notice that the attribute values from the firts created feature are exact like the attribute values from the second created feature.

Here is the code:

mylayer = QgsProject.instance().mapLayersByName('my layer name')[0]
mylayer.startEditing()
feat = QgsFeature ()
feat.setFields(mylayer.fields())

//here is the code that creates the geometry//

feat.setGeometry(geom)
self.iface.openFeatureForm(mylayer, feat, False)
mylayer.addFeature(feat)
mylayer.commitChanges()

What am I doing wrong? I don't know if it helps, but mylayer is an empty layer.


texts:
I extensively use ArcPy, recently needed to apply a specific color ramp on a single band raster and export it. Arcpy does not allow me to apply color ramp and save the raster with the applied color ramp. However, QGIS has the option to do that.

I created a new color ramp in QGIS 2.18.20 and name it to "Yield", which is 28th (defaultColorRampNames[27]) item in the color ramp list.
Is there any PyQGIS codes that can automatically apply the color ramp to the raster. I have the codes to save the raster after applying the color ramp but struggling to find a way to apply that.

Followed the codes from Using QGIS GUI predefined colorramps in Python console / PyQGIS? - but first couples of lines are working not the last three lines. The following codes are working -

#accessing the raster layer
vl = iface.activeLayer()

## get default style
myStyle = QgsStyleV2().defaultStyle()

## get a list of default color ramps [u'Blues', u'BrBG', u'BuGn'....]
defaultColorRampNames = myStyle.colorRampNames()

## setting ramp to "Yield" which was recently created
ramp = myStyle.colorRamp(defaultColorRampNames[27])

But, the codes provided below seems like for a shapefile, not sure what to use for the raster layer

# set up an empty categorized renderer and assign the color ramp
renderer = QgsCategorizedSymbolRendererV2(field, [])
renderer.setSourceColorRamp(ramp)
vl.setRendererV2(renderer)


texts: Want to delete multiple duplicate fields

texts: I have a feature class which contains duplicate fields like a,b,c,a1,b1,c1. I want to delete these duplicate fields. The number of fields is more than 100. So instead of manual deletion I am looking for some automation.


texts:
What is the method used by Raster2pgsql to create overview? I couldn't find details in the docs. Nearest neighbor isn't a very good choice for the dataset I'm working on, but that seem to be the one used, just from visual inspection. I was hoping to be able to coerce something else.
Does it use existing overviews? That would allow me to create them with whatever I like with gdaladdo then import those in PostGIS.


texts: Resolving trade off between smoothing via polynomial fitting and simplify-like algorithms (do not eliminate vertices)

texts:
I am struggling for a few weeks to use Python to smooth various linestrings and bumped into a conundrum. There seems to be a trade-off between using techniques that apply parametric curve fitting and those that apply Non-parametric algorithms like shapely simplify. I have attempted to use the following parametric strategies:
I used numpy polyfit and polyval
coefs_xy = np.polyfit(xcoords, ycoords, 4)
fitted_XY = np.polyval(coefs_xy, xcoords)

that fit polynomial functions to x and y coordinates. And used scipy UnivariateSpline that does a similar process.
spl = UnivariateSpline(xaxis, yaxis, s=.01, k=5)

They work fine if linestrings are well-behaved, i.e., if an x coordinate has only one corresponding y coordinate. This is violated when a linestring is a vertical arch, for example. In this case, the algorithm interprets two y coordinates as a statistical error around a mean and tries to find a 'mean' y coordinate between the two points. Or, in the scipy case, it starts zigzagging between y coordinates. Of course, the smoothing goes to the drain. I have even attempted to invert x with y coordinates and pick the one that has the smallest deviation wrt original linestring. However, there is always an exception that causes problems.
Alternatively, I have applied shapely simplify. It does solve the problem above. It may smooth corners locally, but looking to the whole line, I can still see vertices, especially when you have an almost straight line then a small angle. This is the code I used:
def interpolateBySimplify(lineCoords, crs):
    lineToBeSimplified = LineString(lineCoords)
    denseLine = np.linspace(0, lineToBeSimplified.length, 1000)
    tempVertList = []
    for vert in denseLine:
        tempPoint = lineToBeSimplified.interpolate(vert)
        tempVertList.append(tempPoint)
    interpolatedLine = []
    for i in range(len(tempVertList)):
        interpolatedLine.append(tempVertList[i].coords[0])
    lineToBeSimplified = LineString(interpolatedLine)
    gdf = gpd.GeoDataFrame(columns =['id', 'geometry'], crs=crs)
    new_row = {'id':0, 'geometry':lineToBeSimplified}
    gdf = gdf.append(new_row, ignore_index=True)
    tolerance = 3
    simplified = gdf.simplify(tolerance, preserve_topology=False)
    simplified.buffer(200, join_style=1).buffer(-200, join_style=1)
    finalLine_interpolated = {'id': [1], 'geometry':[LineString(lineToBeSimplified)]}
    return finalLine_interpolated, interpolatedLine

Even with a very large buffer (this is UTM coordinates), I still can see vertices like those in the figure below
Have you come across this problem?
It may be the case I am missing some detail in my approach.


texts: Elevation profile 10 km each side of a line

texts:
How can I get an elevation profile for a band of terrain?

The highest elevation within 10 km (on each side of the defined line) should be taken into account.

I hope my question is clear.
Thank you very much in advance.