texts: Displaying KML Extended Data in OpenLayers?

texts:
I have just got a problem displaying extended data from a KML file in OpenLayers 2. I hope you can help me with this issue. I have already searched for a solution, but I do not find a real solution for it. I have to mention that I am not a programmer and just modify free available scripts and codes.

I have a standard KML file like:

<?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document id="root_doc">

<Schema name="pr14" id="pr14">
    <SimpleField name="Description" type="string"></SimpleField>
    <SimpleField name="Nr." type="int"></SimpleField>
    <SimpleField name="Landkreis" type="string"></SimpleField>
    <SimpleField name="Lage (Koordinate X)" type="float"></SimpleField>
    <SimpleField name="Lage (Koordinate Y)" type="float"></SimpleField>
    <SimpleField name="Kapazität" type="string"></SimpleField>
</Schema>
<Folder><name>P+R pr14</name>
  <Placemark>
    <name>P+R Altenbach</name>
    <description>8011018</description>
    <ExtendedData><SchemaData schemaUrl="#pr14">
        <SimpleData name="Nr.">1</SimpleData>
        <SimpleData name="Landkreis">LL</SimpleData>
        <SimpleData name="Lage (Koordinate X)">51.356068</SimpleData>
        <SimpleData name="Lage (Koordinate Y)">12.681333</SimpleData>
        <SimpleData name="Kapazität">0</SimpleData>
    </SchemaData></ExtendedData>

      <Point><coordinates>12.681333,51.356068</coordinates></Point>
  </Placemark>

</Folder>
</Document></kml>

Now, I would like to display beside of the "name" and "description" field which work fine also the attributes from Extended Data fields.

 var kmllayer2 = new OpenLayers.Layer.Vector("P+R Pl&auml;tze", {
                maxResolution: map.getResolutionForZoom(9),
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "pr14_1.kml",
                format: new OpenLayers.Format.KML({
                    extractStyles: true,
                    extractAttributes: true,
                    maxDepth: 2

To display the name in my popup I use '+feature.attributes.name+', but I cannot sort out how to modify this to get also the extended data values like 'Kapazität" for example.

var controls2 = {
          selector2: new OpenLayers.Control.SelectFeature(kmllayer2, { onSelect: createPopup2, onUnselect: destroyPopup2 })
        };

        function createPopup2(feature) {
          feature.popup = new OpenLayers.Popup.FramedCloud("pop",
              feature.geometry.getBounds().getCenterLonLat(),
              null,
              '

'+feature.attributes.name+'

',
              null,
              true,
              function() { controls2['selector2'].unselectAll(); }
          );

          //feature.popup.closeOnMove = true;
          map.addPopup(feature.popup);
        }

        function destroyPopup2(feature) {
          feature.popup.destroy();
          feature.popup = null;
        }

         map.addControl(new OpenLayers.Control.LayerSwitcher({'ascending':false}));

        map.addControl(controls2['selector2']);
        controls2['selector2'].activate();

Could you give me a hint?


texts: Georeferencing two images at same time in QGIS?

texts:
Has anyone ever tried to georeference two images at the same time on QGIS?

In my application I don't atually care much about the positioning itself, but in the scale and orientation.

I have GeoTIFFs of aerial views of buildings. Unfortunately these images do not have a complete picture of the buildings (low flight drone). I would like to "construct" the complete picture of a building use various parts. These parts are already georeferenced, but due to small errors in the GPS positioning, we do see discontinuation failures on the seams. I don't care about the positioning itself during this construction, as I can later position the complete picture when done.

I'm not sure if I'm using the correct tool, but in QGIS there is the georeferencer. It seems to be able to rotate and move images as I want, but only one image at a time. This isn't ideal for me as I'm trying to assemble what is essentially a puzzle.


texts: No output layer while using GRASS r.viewshed in QGIS

texts:
I have a DEM (in EPSG 3857; same as project EPSG) and try to create a viewshed with GRASS r.viewshed. Function seems to run fine. However, when it is over there is a message that no output layer (tif-file) has been created with. Protocol says

{'output': <QgsProcessingOutputLayerDefinition {'sink':[file path]/view.tif, 'createOptions': {'fileEncoding': 'System'}}>}

Working with QGIS 3.12.0 on Win 10 (x64).

A possible hint: I tried also the Visibility Analysis Pluging. When I tried to create Viewpoints with this tool, I get the following erroer-message:

  Raster pixels are irregular in shape (probably due to incorrect
  projection)!

However, my DEM is in EPSG 3857, a projected CRS - you can't use a geographic CRS for Viewshed analysis, if my information is correct.

Any idea what went wrong and how to fix it?


texts: Creating polygons from open ended line (riverine) feature

texts:
I have a dataset of rivers which I'd like to build into polygons of areas enclosed by rivers.
The rivers are continuous lines but there are open ends.

Is this possible? It's OK if the calculated boundary between the open aspects of the rivers is less precise, and it doesn't matter if the process causes gaps or forms a homogeneous polygon, either is OK.
Currently I'm drawing them by hand using the trace tool. I tried grouping by river/feature name but this wasn't successful. I also tried minimum bounding geometry, both as lines and converted to points, and if I select the points manually then it's fine, but I have no way to group them as again the name isn't helpful.

Is this possible, or should I carry on tracing?


texts:
I am trying to get an overlay (intersection) of two geodataframes of polygons, but it returns me a ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

It's weird as far as I did an intersection some lines before without issue.

clip = gpd.overlay(result_join, concave_hull, how='intersection')
# ^^This one works well
clip.crs = {'init': 'epsg:3857'}

parc = gpd.read_file(parcellaire)
parc = gpd.GeoDataFrame(parc).drop(['name'], axis=1)
parc = parc.set_geometry('geometry')
parc = parc.explode()
parc = parc.buffer(0)
parc.crs = {'init': 'epsg:4326'}
parc = parc.to_crs({'init': 'epsg:3857'})
parc = gpd.GeoDataFrame(parc).rename(columns={0:'geometry'}).set_geometry('geometry')

zonage = gpd.overlay(clip, parc, how='intersection') # This one crashes

Full traceback :

File "/Tim/Google Drive/Fruition Sciences/Dev/Geomatic/Scripts/Zonage/Zonage_standalone.py", line 278, in <module>
    zonage(shape)
  File "/Tim/Google Drive/Fruition Sciences/Dev/Geomatic/Scripts/Zonage/Zonage_standalone.py", line 225, in zonage
    zonage = gpd.overlay(clip, parc, how='intersection') # Intersection du zonage par le parcellaire
  File "/usr/local/lib/python2.7/dist-packages/geopandas/tools/overlay.py", line 175, in overlay
    prop2.drop(df2._geometry_column_name)])
  File "/usr/local/lib/python2.7/dist-packages/pandas/core/reshape/concat.py", line 213, in concat
    return op.get_result()
  File "/usr/local/lib/python2.7/dist-packages/pandas/core/reshape/concat.py", line 372, in get_result
    name = com._consensus_name_attr(self.objs)
  File "/usr/local/lib/python2.7/dist-packages/pandas/core/common.py", line 130, in _consensus_name_attr
    if obj.name != name:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

To reproduce it ---> files

Command: python .../Zonage_standalone.py .../shpdetest 10 .../block.shp


texts: What is the process for identifying a 20 year flood zone?

texts:
I am currently working on an assignment where I need to model flood risk areas. Firstly I was asked to identify potential flood risk zones based on the following criteria:

identify areas that I would consider to have steep slopes. I think I have completed this task. I used the Slope tool here.
Identify a 20 year flood zone and base this on an estimated distance from rivers. Then repeat the task for a 100 year flood return. A guide value would be not to go beyond 2 kilometers for the 100 year and less for the 20 year. I have no idea how to even start this. I have googled different things have come up with extensions that can be used but we never went over any such extensions in the college practicals.


texts: How do I show GWC cached tiles offline in Openlayers 3.5?

texts:
I want to compare performance of my OL app using tiles from GWC and compared to loading them locally from the HD.

Usually I load local tiles that were generated using GDAL with the following code :

var trails= new ol.layer.Tile({
title: 'Bike Trails',
type: 'base',
visible: false,
source: new ol.source.XYZ({
    tileUrlFunction: function(coordinate) {
        if (coordinate == null) {
            return "";
        }
        var z = coordinate[0];
        var x = coordinate[1];
        var y = -coordinate[2] - 1; //3.7.0

        return 'tiles/trailes/' + z + '/' + x + '/' + y + '.png';
    }
})

});

but GWC tiles are generated in different structure / numbers.

GWC struct is :

EPSG4226_Zoom/NNN_NNN/NNNNNN_NNNNNN.png

below the root folder there are many folders named like 076_042 and below them the tiles are named 12 digits with underscore in the middile like :

019470_010907.png

Any idea how to adapt the function to support GWC local tiles ?


texts: gdalwarp malfunctioning for GEOS projection

texts:
I'm trying to reproject a flat lat/lon image into geostationary projection.
The input data is 1441x1441 pixels and covers the region -90 -> 90 lat and lon.
The output data should be 3712x3712 pixels in GEOS projection centred upon 0.0 degrees longitude (i.e: Meteosat/SEVIRI full disk).
I run gdalwarp with this command on a land/sea mask input image:
gdalwarp -t_srs '+proj=geos +a=6378169.0 +b=6356583.8 +lon_0=0.0 +h=35785831.0' -te -5568748.276 -5568748.276 5568748.276 5568748.276 -tr 3000.403165817260742 3000.403165817260742 -multi -r cubic  tmp_in.tif tmp_out.tif

And I get this output (land/sea boundaries added in red):

As you can see, the actual georeferencing looks OK, but parts of the image are missing: Such as Greenland, the western parts of South America and also Antarctica. These are all present in the input data file, but for some reason gdalwarp is cutting them off in the output.
Any ideas what the heck is happening here? I've tried various combinations of -te, -ts and the -t_srs projection without luck.


texts:
I have an orthomosaic, made by drone on with an AI Detector marked all visible plants. The orthomosaic was computed which GCPs collected by emlid reach rs+ with RTK.
As result i got a shape point layer with all plants. The orthomosaic and the resulting point shape layer are in CRS EPSG4326.
From the planing office i got also a  line shape in CRS UTM 33N. This line shape layer contains all information and attributes from each plot of this field. It was used for seeding.
Each plot has two rows of corn. Out of this line shape layer i computed polygons to count the points in each polygon. The process works good so far.
But there is one problem. The line shape fits very good from east to west and is located in the middle of the plots on the walking way. But in the south in runs over the edge of the plots. In the image you can see it. The line runs over the interrupted corn line. The should be in de middle of the wy between the polts and also directly in front of the interrupted rows.
How can i move the polygon layer? I have tried it moving the shape, but when i move the problem occours on the north.
I also tried to reproject all data in UTM 33 but the result was the same.
Is there a probelm with the CRS?
Is there a way to change the polygon grid so the each polygon is also changed relatively?
I hope someone can help[!


texts: Multiple Spatial Point Patterns, create "average" heatmap

texts:
I am using R's spatstat package to analyze spatial patterns of objects across plant leaves.

Let's say I have 10 leaves, each with corresponding X,Y coordinates of the objects and a window in the shape of a leaf. Each leaf is differently sized and each image containing the leaf may have varying dimensions.

I would like to take the separate point patterns, and create one "average" point pattern, with a corresponding heat map. Due to the different sizes, differently shaped/sized windows, is there a way to create a composite pattern from this data?

I have analyzed all the patterns, etc. separately, but would like to be able to view a composite visual pattern. Is there a machine learning exercise I could use?

Apologies if this does not make sense, this is not my area of expertise.

Example tesselations for each leaf:


texts:
I have a string read from a csv file:
'LINESTRING (-79.494722 43.777724, -79.494805 43.777636, -79.49503799999999 43.777581, -79.495277 43.777524)'
I want to convert this string object to
shapely.geometry.linestring.LineString object.
Thousands of lines are in that csv file. Is there a way to convert them?


texts: Why are map tiles being repeated? Using MapServer and Openlayers

texts:
Problem I have: map tiles from my layer are being repeated.

Shapefile is projection 3059.
I have specified 3059 is in the mapfile and also when creating map object.

I know OpenLayers.Layer.WMS is recommended, but when I use it I only see empty map. When I use OpenLayers.Layer.MapServer, map is displayed (and tiles are repeated), so I'm using that for now.

My layer (EPSG:3059) looks like this

Versions:
Openlayers 2.13.1
MapServer 6.0.3

WMS Server mapfile:

MAP
    NAME country_map
    STATUS ON
    SIZE 800 500
    EXTENT 313187.4 146161.56 766469.92 470367.64
    UNITS METERS
    SHAPEPATH "shapefiles"
    IMAGECOLOR 155 155 155

    WEB
        TEMPLATE "template.html"
        IMAGEPATH "C:/MS4W/ms4w/apache/htdocs/gis/output/"
        IMAGEURL "/gis/output/"

        METADATA
            "wms_title"             "MapServer + Openlayers"
            "wms_onlineresource"    "http://localhost:1025/cgi-bin/mapserv.exe?map=C:\MS4W\ms4w\Apache\htdocs\gis\mapfile2.map"
            "wms_enable_request"    "*"
            "wms_srs"               "EPSG:3059"
        END
    END

    PROJECTION
        "init=epsg:3059"
    END

# example layer
    LAYER
    METADATA
        "wms_title"             "lat_rob_region"
    END
    NAME "lat_rob_region"
    DATA lat_rob_region
    STATUS ON
    TYPE POLYGON
        CLASS
            NAME "lat_rob_region"
            STYLE
                COLOR 0 255 0
            END
        END
    END
END

Javascript:

var lon = 25;
var lat = 57;
var zoom = 5;
var map, layer1;

function init()
{
    map = new OpenLayers.Map( 'map',{displayProjection: new OpenLayers.Projection("EPSG:3059")} );

    layer1 = new OpenLayers.Layer.MapServer( "Country map",
            "http://localhost:1025/cgi-bin/mapserv.exe?map=C:/MS4W/ms4w/Apache/htdocs/gis/mapfile2.map", {layers: 'lat_rob_region'},
            {gutter: 15},{srs: 'EPSG:3059'});
    map.addLayer(layer1);

    map.setCenter(new OpenLayers.LonLat(lon, lat), zoom);
    map.addControl( new OpenLayers.Control.LayerSwitcher() );
}


texts: Export CSV from QGIS without thousands separator

texts:
I need to export a GeoPackage file from QGIS 3.2 as CSV file.
The GeoPackage contains two fields: UID as Integer and S_SCH as String (yes, it really is a string).
It looks like:

UID S_SCH
654321 0920056
1256897 0920058

When I export this as CSV, QGIS adds thousand separators to the Integer field:

UID S_SCH
654.321,0920056
1.256.897,0920058

How can I prevent QGIS from doing this? I need it as plain Integer without any separators, just like in the original GeoPackage.

Note this is about exporting a CSV from QGIS, not about importing as all the questions I found so far. Exporting a .csvt won't help here as long as there are still separators in the integer field.


texts: Exporting Sentinel-2A imagery from QGIS to ERDAS Imagine 2013

texts:
I have downloaded a Sentinel-2 image and loaded the 10m bands from the .xml file in QGIS which works fine. I simply want to open the image in Erdas Imagine for use in accuracy assessment in a true colour composite.

In QGIS the bands are as follows: Band 1 (red) 2 (green) 3 (blue) which looks fine.

I then right clicked > saved as a GeoTIFF and opened in Erdas, however when I use the same band combination the colours looks different, as below.

Does anyone know a solution to display the image as seen in QGIS?


texts:
I have a point layer (shapefile) insert to PostGIS. I'm using MapServer and OpenLayers to display layer on map. Now, I want get lon/lat of all features of the layer.

Hope somebody can help me.


texts: Adding data from a SQL Server Spatial database to ArcMap or ArcEngine application with ArcObjects for better performance

texts:
I have an ArcSDE geodatabase stored in MS SQL Server 2008 R2.  The datasets in this geodatabase range from 10 record to 7 million records and there is a mixture of tables and spatial datasets (SQL Geometry).  All datasets are currently related via relationship classes.

It is nearly impossible to work with this data in ArcMap or via ArcObjects because the datasets are so large.  ArcMap hangs when querying or getting related records.  I tried query layers in ArcGIS 10 that query the SQL tables directly, but because these requery the database with every re-draw of the map or identify or attribute table change, the performance was still extremely slow.

Is there a way to query the SQL tables directly once and load the data into ArcMap as a layer?


texts: Georeferencing system of rivers using XY To Line tool?

texts: I have a State Water Hydro layer from our GIS county office that does not have georeferenced attributes (It is sitting on top of a State Plane county basemap). Looking at the XY to Line tool I need Start/End XY fields to run the tool. Problem is that there are no such attributes, and I'm looking at an entire county's worth of streams & rivers. Is there a way to determine these XY values without manually going through every stream feature? Also, I'm trying to suggest to my department why having all of our dept. specific features georeferenced is useful for future analysis, but I'm a bit green in the when it comes to justifying the ends to some of our work. Any suggestions on top of the XY to Line tool help would be great.


texts: How to calculate the Vertical Distance to Channel Network in ArcGis?

texts:
In Saga there is the tool "Vertical Distance to Channel Network". Is there a way to produce the same data in ArcGIs?
With con("flowacc" > xx, 1) I was able to create the Channel Network.
I am guessing that I have to subtract the channel networks elevation from the original elevation, but I am not sure how to do so.


texts:
My ajax request for a geojson works fine as a static web app

  $.getJSON("myfile.geojson", function (data) {
    theaters.addData(data);
    map.addLayer(theaterLayer);
  });

I am trying to wrap the project now around flask.  All of my js/geojson files are in my static folder.

However when I try to pass static folder url to geojson inside flask (with map.js also residing in static/js) like such:

  $.getJSON("/static/js/myfile.geojson", function (data) {
    theaters.addData(data);
    map.addLayer(theaterLayer);
  });

I get 'Cannot read property 'addLayer' of undefined'.

This is the only file I load from the js.  What is the proper way of referencing local files inside my map.js in a flask environment?


texts:
I found a topic on Esri ArcGIS site about using ArcGIS 3D Analyst for designing and placement of surveillance cameras, but there isn't any clear explanation to guide designer.

Is there any idea or guidance to use ArcGIS in this subject?


texts: gdalwarp `-te -5.8 41 10 51.5 -ts 1980 0` not working with mercator?

texts:
I was working with native lat/log projection on etopo1:

#donwload:
curl -L -C - 'http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/georeferenced_tiff/ETOPO1_Ice_g_geotiff.zip' -o ./ETOPO1.zip
unzip -n ./ETOPO1.zip '*.tif'
#commands:
gdal_translate -projwin -5.8 51.5 10 41 ../data/noaa/ETOPO1_Ice_g_geotiff.tif cropXL.tmp.tif
gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:4326 cropXL.tmp.tif reproj.tmp.tif
gdalwarp  -te -5.8 41 10 51.5 -ts 1980 0   reproj.tmp.tif resized.tmp.tif

all working fine.

I now wish to work in mercator. Also, I made a tiny change to reproject to mercator -t_srs EPSG:3857, but it doesn't work:

gdal_translate -projwin -5.8 51.5 10 41 ../data/noaa/ETOPO1_Ice_g_geotiff.tif cropXL.tmp.tif
gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:3857 cropXL.tmp.tif reproj.tmp.tif
gdalwarp -te -5.8 41 10 51.5  -ts 1980 0   reproj.tmp.tif resized.tmp.tif

The reprojection (2) seems fine. The last command (3) runs but resized.tmp.tif is now flat, all black. No error message. But the values x,y,z seems non-sense.

What am I doing wrong ?


texts: SelectLayerByLocation works in ArcMap but not in standalone script

texts:
I have a piece of code that selects a layer based on the data frame extent.
the code works when I use the code within ArcMap.  However, as soon as I try run this as a stand alone script, nothing gets selected. So some empty Excel files get created.  I have no understanding of why running this as a standalone script would be any different.

import os
import arcpy
import sys

#create polygon for DF extent
def getDataFramePoly(mxd):
    mxd = mxd
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    extent = df.extent
    array = arcpy.Array()
    array.add(extent.lowerLeft)
    array.add(extent.lowerRight)
    array.add(extent.upperRight)
    array.add(extent.upperLeft)
    array.add(extent.lowerLeft)
    polygon = arcpy.Polygon(array)
    return polygon

#get relevant layers function
def getLayers(layers, datasources):
    analysis_layers = []
    layers = layers
    for lyr in layers:
        if lyr.supports('DATASOURCE') and str(lyr.dataSource)in datasources:
            analysis_layers.append(lyr)
    return analysis_layers

mxd = arcpy.mapping.MapDocument(r"G:\mxdName.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0] #First Dataframe
extent = df.extent

#select layer by DF
df_extent = "DF_EXTENT"
polygon = getDataFramePoly(mxd)
arcpy.MakeFeatureLayer_management(polygon,df_extent)

#feature classes that should be included in analysis
dataSources = [r'datasourcePath' r'datapayt']

analysis_layers = getLayers(arcpy.mapping.ListLayers(mxd,"",df), dataSources)
excel_out_path = r'C:\Users\Bbrokling\Desktop\Test'

#Select visible features
i = 0
for lyr in analysis_layers:
    arcpy.SelectLayerByLocation_management(lyr,"INTERSECT",df_extent)
    desc = arcpy.Describe(lyr)
    if not desc.FIDSet  == '':
        print 'Layer has a selection'
    if lyr.symbologyType == u"UNIQUE_VALUES":
        desc = arcpy.Describe(lyr)
        symField = lyr.symbology.valueField
        symClassValues = lyr.symbology.classValues
        otherValues = lyr.symbology.showOtherValues
        if not otherValues:
            for x in desc.fields:
                if x.name == symField:
                    print x.type
                    fieldType = x.type
            if fieldType in ['SmallInteger', 'Double', 'Integer']:
                newclassvalues = tuple([float(x) for x in symClassValues])
            else:
                newclassvalues = tuple([x.encode('UTF8') for x in symClassValues])
            newQuery = "{} in {}".format(symField, newclassvalues)
            arcpy.SelectLayerByAttribute_management(lyr,"SUBSET_SELECTION",newQuery)
    out = r'%s\%s.xls' % (excel_out_path,lyr.name)
    desc = arcpy.Describe(lyr)
    if not desc.FIDSet  == '':
            print 'Layer has a selection'
    table_view_Name =  "tableView%s" % str(i)
    arcpy.MakeTableView_management(lyr, table_view_Name)
    arcpy.TableToExcel_conversion(table_view_Name, out)
    i += 1


texts: OpenLayers load WFS vector layer: how to check whether all features are returned from server?

texts:
My app takes several seconds to return all features of a WFS layer from server, I am wondering whether there is a way to check whether all features are returned in coding, so that we can add a loading-mask...

In other words, how do you check whether a WFS request is successful or not?

Code sample:

var newEventLayer = new OpenLayers.Layer.Vector
(
    'test'
    , { "strategies":
            [new OpenLayers.Strategy.BBOX()
            , new OpenLayers.Strategy.Refresh({ "force": true })
            ]
        , "protocol": new OpenLayers.Protocol.WFS(
        {   "version": '1.0.0'
            , "srsName": 'EPSG:4326'
            , "url": '/geoserver/ows?service=WFS'
            , "featureType": 'XXX'
            , "featureNS": 'XXXX'
            , "geometryName": 'GEOM'
        })
    }
);

map.addLayer(newEventLayer);


texts: Openlayers 3 screws up rendering of vector data on high zoom

texts:
I'm working on a web app that needs to show a lot of detail.

It's also loading a geoJSON file as background, but I have problems with it when zooming in. It's like the vector renderer doesn't find the off-screen point, and just renders some horizontal line.

See the two images, where the first is zoomed in to the deepest level, and the second is zoomed out one level (so shows the features correctly).

I've tried to enlarge the render buffer of the vector layer, but to no avail, even not with crazy values like 20000 (instead of the default 100).

Does anyone have an idea on how to let OL3 find the off-screen point without changing the maximum zoom level? Preferably also without changing the geoJSON but perhaps a specification of how far points on a line may be apart might work.


texts:
Still learning OpenLayers and I think my problem may be a lack of understanding the terminology, however.

Is it possible to draw polygons in a different colour within one vector layer, or do all polygons get their colour from the colour set within the layer?

I ask because I have a map which draws the borders of several (up to 10) different countries around the world and I'd rather have one layer with different coloured polygons rather than have several different layers.

I have the following code:

var polyCountry = format.read("POLYGON((" + polyString + "))");
Zip15021.addFeatures([polyCountry])

Which works perfectly, but can't see any documentation anywhere on how to style the polygon.


texts:
I would like to prepare a confusion matrix from  decision tree module in R.

Can anyone help me with the code that is needed for calculating user, producer and overall accuracy along with kappa values in R?

My confusion matrix looks like this

                       Leaf-off Trees Leaf_on Trees
  Leaf-off Trees             41            10
  Leaf-on Trees               1            46


texts: Convert SHAPE_AREA units via the REST FeatureServer API

texts:
All similar questions seem to involve desktop software but I'm strictly using the arcgis API like so:
https://services2.arcgis.com/API_KEY_HERE/arcgis/rest/services/Category/FeatureServer/0/query?where=&objectIds=&time=&geometry=-62.67407%2C41.76575&geometryType=esriGeometryPoint&inSR=4326&outFields=*&outSR=4326&f=pgeojson&token=

I change the Output Spatial Reference to 4326, which converted the coordinates successfully in the response, but the SHAPE_AREA is an unknown unit. My goal is to get the SHAPE_AREA in sqm or sqft.
Right now that query returns something like: "Shape__Area" : 210.56640625
I'm assuming I need to project/convert it somehow like I did with the coordinates. Any advice?
The metadata for the Geoservice says this:

Dataset Source: GIS Dataset. Source Format: GIS SDE Feature Class

and this:


texts:
I'm trying to calculate center line from polygon, so i'm trying to do  with postgis function st_straightskeleton but it is rising the error of function not found.

postgis version: POSTGIS="2.4.3 r16312" PGSQL="100" GEOS="3.5.1-CAPI-1.9.1 r4246" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" LIBPROTOBUF="1.2.1" (core procs from "2.4.2 r16113" need upgrade) RASTER (raster procs from "2.4.2 r16113" need upgrade)

i have tried this example but getting error

SELECT ST_StraightSkeleton(ST_GeomFromText('POLYGON (( 190 190, 10 190, 10 10, 190 10, 190 20, 160 30, 60 30, 60 130, 190 140, 190 190 ))'));

Error

ERROR:  function st_straightskeleton(geometry) does not exist
LINE 1: SELECT ST_StraightSkeleton(ST_GeomFromText('POLYGON (( 190 1...


texts: Coordinates transformation of four columns in qgis

texts:
I have four columns with grid coordinates X_min Y_min X_max Y_max. As these coordinates are in ESPG 3857, i want to convert each of the columns to 4326 in QGIS.

Can anyone help?


texts: QGIS: GPS Tracking with improved Altitude

texts:
I am preparing for archeological fieldwork in Sudan. I recently came up with RTK equipment and software using RTKlib to export the real time corrected data (cm level) to QGIS via a simulated NMEA device.

However I checked the NMEA data flow and it includes elevations down to mm; in QGIS only one decimal is given. Is there any way to improve this and have the altitude information shown in cm as well?

Also QGIS applies geoid corrections but I wish to use ellipsoidal elevation. And finally right now I have to include the elevation data manually in the datamask for every point while saving. Is there a way QGIS can save the elevation in points and polylines?


texts: How to speed up an st_intersects() query?

texts:
I have a "huge" spatial database: 72GB, 200 millions documents with geometries (points)
I have a GIST index on the geometries

the database:

CREATE TABLE resources
(
id bigint NOT NULL,
accuracy varchar(3),
owner varchar(100),
title varchar(500),
posted timestamp,
nbtags int,
url text,
description text,
rawtags text,
tags text,
geom geometry NOT NULL,
CONSTRAINT resources_pkey PRIMARY KEY (id)
);
 INDEX resources_index_rep ON resources USING GIST (geom);

CREATE RULE resources_rule AS ON INSERT TO resources
WHERE EXISTS(SELECT 1 FROM resources
            WHERE id =NEW.id)
DO INSTEAD NOTHING;

Of course I want to query this database. Nevertheless, it's very very slow ...
(I did analyzes and Vacuums)

for example with area of France:

EXPLAIN ANALYZE
SELECT count(*)
FROM resources
WHERE st_intersects(st_geomfromtext('POLYGON((-5.134723 41.364166,-5.134723 51.09111,9.562222 51.09111,9.562222 41.364166,-5.134723 41.364166))',4326), geom)

QUERY PLAN
Aggregate  (cost=14411182.45..14411182.46 rows=1 width=0) (actual time=17368911.185..17368911.186 rows=1 loops=1)
 ->  Bitmap Heap Scan on resources  (cost=890493.27..14391392.83 rows=7915848 width=0) (actual time=2627029.235..17357552.440 rows=21381758 loops=1)
    Recheck Cond: ('0103000020E6100000010000000500000077137CD3F48914C05628D2FD9CAE444077137CD3F48914C0562B137EA98B4940F52EDE8FDB1F2340562B137EA98B4940F52EDE8FDB1F23405628D2FD9CAE444077137CD3F48914C05628D2FD9CAE4440'::geometry && geom)
    Filter: _st_intersects('0103000020E6100000010000000500000077137CD3F48914C05628D2FD9CAE444077137CD3F48914C0562B137EA98B4940F52EDE8FDB1F2340562B137EA98B4940F52EDE8FDB1F23405628D2FD9CAE444077137CD3F48914C05628D2FD9CAE4440'::geometry, geom)
    ->  Bitmap Index Scan on resources_index_rep  (cost=0.00..888514.31 rows=23747545 width=0) (actual time=2626090.609..2626090.609 rows=21381788 loops=1)
          Index Cond: ('0103000020E6100000010000000500000077137CD3F48914C05628D2FD9CAE444077137CD3F48914C0562B137EA98B4940F52EDE8FDB1F2340562B137EA98B4940F52EDE8FDB1F23405628D2FD9CAE444077137CD3F48914C05628D2FD9CAE4440'::geometry && geom)
Total runtime: 17368946.473 ms
7 row(s)

it took almost 5 hours.

and with Switzerland:

EXPLAIN ANALYZE
SELECT count(*)
FROM resources
WHERE st_intersects(st_geomfromtext('POLYGON((5.96611 45.829437,5.96611 47.806938,10.488913 47.806938,10.488913 45.829437,5.96611 45.829437))',4326), geom)

QUERY PLAN
Aggregate  (cost=5623425.68..5623425.69 rows=1 width=0) (actual time=6798842.064..6798842.064 rows=1 loops=1)
->  Bitmap Heap Scan on resources  (cost=88982.29..5621452.65 rows=789212 width=0) (actual time=783383.496..6797250.482 rows=2539057 loops=1)
    Recheck Cond: ('0103000020E61000000100000005000000AF5A99F04BDD1740D28BDAFD2AEA4640AF5A99F04BDD174028F38FBE49E74740B22D03CE52FA244028F38FBE49E74740B22D03CE52FA2440D28BDAFD2AEA4640AF5A99F04BDD1740D28BDAFD2AEA4640'::geometry && geom)
    Filter: _st_intersects('0103000020E61000000100000005000000AF5A99F04BDD1740D28BDAFD2AEA4640AF5A99F04BDD174028F38FBE49E74740B22D03CE52FA244028F38FBE49E74740B22D03CE52FA2440D28BDAFD2AEA4640AF5A99F04BDD1740D28BDAFD2AEA4640'::geometry, geom)
    ->  Bitmap Index Scan on resources_index_rep  (cost=0.00..88784.99 rows=2367636 width=0) (actual time=782664.847..782664.847 rows=2539063 loops=1)
          Index Cond: ('0103000020E61000000100000005000000AF5A99F04BDD1740D28BDAFD2AEA4640AF5A99F04BDD174028F38FBE49E74740B22D03CE52FA244028F38FBE49E74740B22D03CE52FA2440D28BDAFD2AEA4640AF5A99F04BDD1740D28BDAFD2AEA4640'::geometry && geom)
Total runtime: 6798852.955 ms
7 row(s)

it took almost 2 hours.

I'm looking for how to improve the performance?

About our server:

- CPU: 4
- RAM: 16 GB
- OS: Ubuntu 12.04 LTS 64bit


texts: ArcMap stops displaying feature classes when zoomed to 1:846 or larger scale

texts:
There is a specific zoom-in point at which all of the layers in data view stop drawing and only the basemap is visible.  This includes geodatabase feature classes, shapefiles, raster images, everything but the basemap disappears.  The specific zoom-in point is 1:846.  At 1:847, everything displays properly.  At 1:846 (or larger) scale, everything disappears but the basemap.  This has been ongoing for a couple years and we've worked around it but it is a problem for us because we create detailed barn and building layouts for farms, and sometimes the farms are quite small and call for a large scale to show the detail.

It is happening to multiple users and we are all using ArcMap 10.6.1

Any ideas on how we can fix it?


texts:
I would like to think with all the animations out there thsi would be relatively simple but it hasn't been.

I would like to animate the movement of a stream (polyline) from one end to another, with all its forks and tributaries (smaller polylines).  I have been reduced  to trying the Time Stamp property in ArcGIS 10.  I have a number of questions related to this, such as why is it that the Time Stamp property always defaults to hours, or that if I try to alter the time property to a different time step, it causes the one segment I have been working with to disappear?

If anyone has worked with the time stamp property in ArcGIS 10 successfully for a whole layer I would love to hear from you.

Dan Brenner
Eureka, CA


texts: GeoDjango: order PointField by distance?

texts:
I've got a GeoDjango model as follows:

class City(models.Model):
    id = models.IntegerField(primary_key=True)
    county = models.ForeignKey(County)
    location = models.PointField(null=True, blank=True)

How do I get all the cities within a particular bounding box, ordered by distance? I've got the relevant latlngs as strings, but I'm not really sure how to construct the bounding box, or the central point, or the query :)

EDITED: have now worked out how to construct the bounding box and get places - but can anyone help me sort by distance?

swLat = request.GET.get('swLat')
swLng = request.GET.get('swLng')
neLat = request.GET.get('neLat')
neLng = request.GET.get('neLng')
centreLat = request.GET.get('centreLat')
centreLng = request.GET.get('centreLng')
centre = fromstr('POINT(%s %s)' % (centreLat, centreLng))
boundaries = fromstr('POLYGON((%s %s,%s %s,%s %s,%s %s,%s %s))' % \
       (swLng, swLat, swLng, neLat, neLng, neLat, neLng, swLat, swLng, swLat))
cities = City.objects.filter(location__within=boundaries)
# how to order cities by distance from centre?

Thanks for helping out a beginner!


texts: Radiance and Raw Image difference after Radiometric calibration

texts:
I have done Radiometric calibration on Envi 5.3 and got the output of my input Landsat 5 image.
What is the difference between the raw image and the radiance image and why is this value of radiance called as
TOA radiance?

Figure 1: Radiance Image

Figure 2: Raw Image


texts:
Is it possible to retrieve osm data by using overpass api, for a specific boundary area? I have a custom area in poly format (can change it to geojson) and I would like to retrieve the ways edited by specific users.
Is it possible to find the same as above, by other means?


texts:
I would need to predict wind speed and direction at a wind farm station, using an artificial neural network. Unfortunately, the only historical weather data available are the one at some weather stations close to my target point. Thus, I will be able to predict the short term wind speed and direction at those sites and then I will need to interpolate those predicted value in order to find the forecast for the target site.

Can you please suggest me how to do it?

Should I use ArcGIS?

I read about the use of Inverse distance weighting or Kocriging method but I don't know how to use them.


texts: Looking for guidance with beginner use cases

texts:
Background

I'm more or less completely new to GIS. At work we use a program called Microvision which does stuff like work out drive distance (in addition to straight line distance) from which you can then create distance bands. Plus has some other mapping features like creating layered ranges. However it has got to be one of the worst programs I have ever had the displeasure of using (not in terms of features, but just in the way it has been coded).

So I'm looking at open source alternatives instead. My problem is that I really have no idea where to start. I downloaded SAGA but quickly realised how far out of my depth I am. The features SAGA boasts are cool, but in some respects it is may more powerful than I need. I also don't understand how to import a base map into SAGA. I found some demo maps which are cool, but I have no idea how to create those maps in the first place. Someone on another forum suggested Quantum GIS, so I've going to give that a go very soon.

My main needs/desires are below, if you can suggest something to match them, I'd be really grateful:

A prebuilt map of the world (or least Australia), or some process to follow to create one. This is so I can then overlay our own maps and views over top of it.
Something with a lot of customisation and preferably a good API. An ability to manipulate and transform data within the program would be great (even if it requires grunt work coding to do). Microvision is one of the most rigid pieces of software I've ever used and doesn't let you customise anything.
Typically my data will be at a postcode level, so the ability to map postcodes to latitudes and longitudes would be awesome. Something that comes with postcode catchments, or even census collection districts for Australia would be AMAZING.
Drive time, aka googlemaps style. This isn't essential, but it would really help us for some of the work we do. Of course I can do approximations using straight line distance, but Australia isn't the flattest of countries so something that can look at drive distances would be fantastic.

Help appreciated :)


texts: Having new points inserted into shapefile automatically display using ArcPy?

texts:
Displaying points from a table is easy - you just right click the table and choose "Display XY Data".

But what is the comparable method for shapefiles?  When you right click a shapefile the "Display XY Data" option is nowhere to be found.  The "Add XY Coordinates" tool is not appropriate because it just adds two new columns "POINT_X" and "POINT_Y" to the table, not to the map.

My shapefile's attributes table is as follows.

The first 4 points were pre-existing... then I ran a script that added the last two from an excel table.  So now they're in my shapefile's attributes table, but not showing on the map.

Relevant portion of My script:

#Insert NEW points into Target Shapefile
x = 0
with arcpy.da.InsertCursor(Target_Shapefile, ['CSSP_ID', 'Lat', 'Lon']) as curs:
    while x < total_new:
        curID = new[x]
        newcoord = latLonDi[curID]
        newlat = newcoord[0]
        newlon = newcoord[1]
        curs.insertRow((curID, newlat, newlon))
        x+=1

How do I display these last two points so that I have all six showing on the map?  This seems like it should be a straightforward task for Arc.  Surely there's a relatively simple process akin to "Display XY Data" for tables.


texts:
I am using EnterpriseDB to run Postgis.

when I run "shp2pgsql" on got an error message.

~ macbook$ shp2pgsql
-bash: shp2pgsql: command not found
~ macbook$ which shp2pgsql

I read that I should find it in this directory, but It seems that I am missing that file.

~ macbook$ ls /Library/PostgreSQL/9.3/bin/

clusterdb       pg_controldata      pgbench
createdb        pg_ctl          pgsql2shp
createlang      pg_dump         pltcl_delmod
createuser      pg_dumpall      pltcl_listmod
dropdb          pg_isready      pltcl_loadmod
droplang        pg_receivexlog      postgres
dropuser        pg_resetxlog        postmaster
ecpg            pg_restore      psql
initdb          pg_standby      raster2pgsql
oid2name        pg_test_fsync       reindexdb
pg_archivecleanup   pg_test_timing      shp2pgsql
pg_basebackup       pg_upgrade      vacuumdb
pg_config       pg_xlogdump     vacuumlo

when I run this comma
/Library/PostgreSQL/9.3/bin/shp2pgsql.sh Outlines.shp city | psql -d City

It says:
-bash: psql: command not found
-bash: shp2pgsql.sh: command not found

How can I run this command ?
should I go for shp2pgsql-gui loader ?


texts:
Does anyone know if there is a shapefile available which moves Alaska and Hawaii closer to the continental USA?

I have attached a photo which might help.


texts: Separating points with identical coordinates using ArcPy?

texts:
I need to move apart points with identical coordinates.

I have a point layer where there can be as high as 4 points on each other (same coordinate). I want to separate these points from each other (with a defined distance in x and y), using a pyscript. Right now, only the script's first part (Part 1) works, but then for some reason the script stops moving points at the second "cursor loop".
Is there a better way to solve the task?

#------------------------------------------------------
#Enviroments
    import arcpy
    arcpy.env.overwriteOutput = True
    workspace = r"*"
    arcpy.env.workspace = workspace
    fc = "*"
    fields = ['SHAPE@X', 'SHAPE@Y']

    #Copy feature class:
    arcpy.CopyFeatures_management(fc, fc + "Copy")
    #Set FC to Copied data
    fcCopy = fc + "Copy"

    edit = arcpy.da.Editor(workspace)

    ----------------------Part 1----------------------

    Xlist1 = []
    Ylist1 = []
    edit.startEditing(False, False)

    edit.startOperation()
    with arcpy.da.UpdateCursor(fcCopy, fields) as cursor:
        for row in cursor:
            if row[0] in Xlist1 and row[1] in Ylist1:
                row[0] = row[0] + 10.0 #Move the point 10 meters in X
            if row[0] not in Xlist1 and row[1] not in Ylist1:
                Xlist1.append(row[0])
                Ylist1.append(row[1])
            cursor.updateRow(row)
        print(Xlist1, Ylist1)
        del row
        del cursor
    edit.stopOperation()
    edit.stopEditing(True)

    print("Script Part Done 1")

#----------------------Part 2 ----------------------

Xlist = []
Ylist = []

    edit.startEditing(False, False)
    edit.startOperation()
    with arcpy.da.UpdateCursor(fcCopy, fields) as cursor:
        for row in cursor:
            if row[0] in Xlist and row[1] in Ylist:
                row[1] = row[1] + 10.0   #Move the point 10 meters in Y
            if row[0] not in Xlist and row[1] not in Ylist:
                Xlist.append(row[0])
                Ylist.append(row[1])
            cursor.updateRow(row)
        del row
        del cursor
    edit.stopOperation()
    edit.stopEditing(True)

    print("Script Part Done 2")

From

To:


texts: How do you convert SRTM HGT elevations from WGS84 into meters above sea level?

texts:
I am building a JavaScript script that reads binary HGT files - as generated by the SRTM - and converts the data into bitmaps viewable in your browser.

You can try out the script and view the source code of a preliminary version on GitHub:

https://github.com/jaanga/terrain-plus/tree/gh-pages/cookbook/read-hgt-files

The app seems to be reading the elevation data correctly, but the numbers are huge ( over 20>). See the demo which shows the first 50 heights of the file. In this case, the demo shows the southern tip of Point Reyes near San Francisco CA.

The SRTM 'Quickstart Guide' says:

  Heights are in meters referenced to the WGS84/EGM96 geoid

I cannot find a definition or explanation for what this exactly means. Nor have I found any JavaScript code that could be used to convert the 16 bit integer HGT data into something like meters of elevation above or below sea level.


texts: OpenLayers take user input of XY and create a point

texts:
I'm looking to take input from a form (X and Y WGS84 coordinate) and then draw a point on my map when the user clicks execute.

I could not locate something like this in the OpenLayers examples.


texts: Can QGIS2 refresh layers without the need to reclassify each layer individually?

texts:
I have a set of layers which I'd like to have refreshed each time I open QGIS2.2.0 (Windows 7).  Each of 14 layers (polygons) is linked via a join to a csv file with no geometry. The layer displays a graduated fill based on data in the csv file. The CSV file is in the layer tree as well. The original csv file is update daily (using an R script).  So my procedure is to run the R script which updates the csv file, then open QGIS.  The csv files in qgis are in fact updated - i can see the new data in the attribute table, but the layer's gradient fill is  not updated, and I have to go through each layer properties and reclassify to see the new data displayed.

I thought there might be a setting, but I can't find one.  Or perhaps there is a better way to organise this project - maybe a database would be a better option, but I'm not clear how to set this up so I can retain the same level of automation.  Thanks for you help!


texts: Convert a collection of raw Landsat 8 scenes to top-of-atmosphere reflectance

texts:
This link , use ee.Algorithms.Landsat.TOA(input)  to convert a raw Landsat scene to top-of-atmosphere reflectance (TOA). I am trying to convert a collection of raw Landsat 8 scenes (LANDSAT/LC8
) to top-of-atmosphere reflectance (TOA):

//Filter USGS Landsat 8 Raw Scenes base on a point location (geometry)
//then Apply TOA algorithm to the raw collection:
var rawcollection = LC8.filterBounds(geometry);
var TOA = ee.Algorithms.Landsat.TOA(rawcollection);

Earth Engine reply that this is  Invalid type. Expected: Image.
So how to apply ee.Algorithms.Landsat.TOA to an ImageCollection?


texts: How can I control the zoom level of an OpenLayers layer in QGIS?

texts:
When you go to a web map, the map server serves a set of tiles that control the detail displayed and the font size of feature labels based on the zoom level.

I am using the QGIS OpenLayers plugin to add a basemap in a project. The plugin must be choosing some zoom level to request for display. However, the zoom level chosen seems inappropriate to the area I am trying to create a print layout for. There are too many features labelled and the font sizes are too small (in some cases to the point of illegibility). See for example the difference between the zoom level chosen in the OpenLayers Overview pane (lower left), which seems reasonable for the area of interest, and the zoom level that appears in the main map (right).

How can the zoom level be adjusted when using the QGIS OpenLayers plugin?


texts: How to divide a shapefile based on a 500x500m grid clip areas based on a grid?

texts:
I want to divide a shapefile based on a 500x500m grid into a series of individual shapefile using the clip function, each only having the contents of one of the 500x500 areas.

May I ask how to do this?

The suggestions I found are primarily for using a single shapefile as a clip shape to cut through multiple layers. But there is no method to do this.

As shown below, I have a few building and a grid of 200m:

I want to use each tile in the grid to clip the building layer and save the clipped part as a new layer:


texts:
For the past couple weeks I have been trying to view GOES17 data from netCDF files (converted to geotif) in QGIS 3.10 but cannot get the projection to work correctly.

I have attempted numerous methods but most recently tried the procedure of the top answer in this post:
Converting NetCDF dataset array to GeoTiff using rasterio Python

When loading the .tif into QGIS, it appears in the wrong place relative the OpenStreetMap. I have tried several projections including EPSG:3857, which is what appears in the bottom right in QGIS when the OSM is loaded. They are all wrong.

I have also tried this answer:
How do I add projection to this NetCDF file? (Satellite)

When attempting the reproject function I get an error.

xds3857 = xds.rio.reproject("epsg:3857")

Error:

DimensionError: x dimension not found. 'set_spatial_dims()' can address this.

xds:

<xarray.Dataset>
Dimensions:                                 (number_of_LZA_bounds: 2, number_of_SZA_bounds: 2, number_of_image_bounds: 2, number_of_time_bounds: 2, x: 1086, y: 1086)
Coordinates:
    t                                       datetime64[ns] 2020-02-03T19:05:05.476645888
  * y                                       (y) float32 0.1519 ... -0.15190002
  * x                                       (x) float32 -0.1519 ... 0.15190002
    goes_imager_projection                  int32 -2147483647
    y_image                                 float32 0.0
    x_image                                 float32 0.0
    retrieval_local_zenith_angle            float32 85.0
    quantitative_local_zenith_angle         float32 70.0
    solar_zenith_angle                      float32 180.0
    time                                    int32 -2147483647
    spatial_ref                             int64 0

The issue continues to persist after doing the suggestion.

xds.rio.set_spatial_dims("x","y",inplace=True)


texts: I am using the template picker from the ESRI JS API.  As far as I see it, the look and feel of the picker, cannot be customized very much at all.  The rows and columns can be controlled, and in the API reference it mentions a style parameter, which does not seem to be able to do very much.  What i would like to do, is to display the features in the Template picker in more of a list style.  As you see in the first image (to the right), the patch is shown on the top of the image with the description below that.  This style takes up a lot of real estate on the webpage.  I would prefer to horizontally align the patch and the description, like you would see in a typical legend.  is this possible, or does anyone know of another method to use a template picker like menu to activate the Editor in JS?


texts:
I have two shape files with lines: one has 41k+ lines and the other has 606 lines. Lines from these layers often overlap. I want to use the second one as a mask to clip lines from the first one so I have a layer with around 40,5k with overlapping features deleted.
I have tried "Difference" and "Symmetrical difference" operations in QGIS but they didn't work for me. How can I do it then?


texts: Installing the QGIS C++ API for Visual Studio?

texts:
I've recently being trying to get the C++ API working in Visual Studio, but there really not allot of information out there.

I've found these code examples https://github.com/ranchodeluxe/qgis_cpp_plugin_cookbook

But I just don't know how to install the API for it to work.


texts: Represent time with many maps and standard graduation

texts:
I want to represent trends over time.

So I have 4 diferent maps (e.g. 95' 00' 05' 10') but my problem is in graduation.

Is there any form to make the graduation with standars breaks and colors (e.g. in all of the maps 30-60 to be red, 60 to 90 to be blue etc)

I read many post like that one:
 https://www.gislounge.com/time-and-gis/
but i couldnt found my answer.


texts:
I have a set of sensors at fixed locations across a landscape. Each is connected to an Arduino Uno board which is coupled to a WiFi board. New sensors get added occasionally and old ones stop working. I want to display the data gathered by each sensor, in (nearly) real-time, in QGIS, as a map of points with attributes displayed as (for example) colours of point symbols. The data value is read by each sensor every 15 minutes or so and thus I want to upload these values to QGIS and update my QGIS vector map of the sensor values at this interval.

The data themselves are very simple - a single real numerical value which I could convert and present as integer, ASCII text string, HTML, XML or pretty much anything else because I can convert them at the sensor end. Since I know the location of each sensor I can also include values of lat/long (for example) in the data each sensor URL presents to the world.

What's the very simplest way to get my data from the sensors onto a QGIS map which periodically gets the data and displays them?


texts: How to insert geometries from shapefile into existing table?

texts:
I'm running into a problem trying to figure out how I should insert a GIS coverage polygon into an existing table. I'm primarily using QGIS and the PostGIS Manager plugin to handle it. After some trial and error, I believe I have my table set up correctly to accept the polygon object that I want to link to my existing table. Below is the schema, pruned down a bit -- the actual table has more traditional SQL columns that are not applicable to this.

However, effectively, I want to enable a spatial query against each of the current rows by using a polygon -- i.e., query against the rows by determining if a point is contained within the polygon.

CREATE TABLE "distribution"
(
  location character varying(5) NOT NULL,
  unit smallint NOT NULL,
  "unitName" character varying(15) NOT NULL,
  coverage geometry,
  CONSTRAINT "distribution_pkey" PRIMARY KEY (location, unit),
  CONSTRAINT enforce_dims_coverage CHECK (st_ndims(coverage) = 2),
  CONSTRAINT enforce_geotype_coverage CHECK (geometrytype(coverage) = 'POLYGON'::text OR coverage IS NULL),
  CONSTRAINT enforce_srid_coverage CHECK (st_srid(coverage) = 4326)
)

Given the above table and a ESRI shapefile, how do I actually go about inserting the shapefile into the appropriate row's coverage field?


texts: Tiling from complex, multi-step GDAL VRT?

texts:
I've got 1002 ~100 MB (5000'x5000') .tif files in a state plane (transverse merc) projection comprising a county-wide, 1-ft aerial dataset. After two days, gdal_merge.py is only at 93 of 1002. And when it's done, I'll still need to apply a mask/cutline, then reproject to EPSG:3857 before I can even start tiling. So I'm trying to find a quicker approach..

I'm relatively familiar with the VRT format and this idea of lazy evaluation, but after several rounds of testing, I'm not convinced the workflow I'm pursing is more efficient.

Here's the workflow I'm asking about (skipping the mask/cutline for now):

gdalbuildvrt --> merge the 1002 .tif files -------------> outputs mosaic VRT
gdalwarp -----> transform to EPSG:3857 ------------> outputs web merc VRT
..start tiling against the 2nd-generation VRT..

The tiling script I wrote in Python (which may be a weak-link in its own respects) is definitely capable of reading the VRT. Small-area windows render in a few seconds, still longer than I'd like. But large area windows take quite a long time. As I'm writing on this, I'm waiting on a sample z=12 tile to render.. I bet it's already eaten 20+ minutes.

Ultimately, my goal is to create 256x256 tiles in the OSM scheme, like you'd get out of Mapnik, but using a Python/GDAL-only approach. ..I've deemed this necessary because I can't get Mapnik to properly consume my 2nd-generation VRT.

My python script is creating the OSM grid properly---I sniped some code from a Mapnik utility (generate_tiles.;y), and I'm getting proper cell coordinates for my raster. But I question whether or not my ReadRaster() / WriteRaster() implementation is efficient.. I'll paste in that part; values are hard-coded for my sample z=12 tile..

def getTileData():
    gdal.AllRegister()
    dataset = gdal.Open('C:/xGIS/Aerial/3857/2009_1ft_od_3857.vrt', GA_ReadOnly)
    bands = dataset.RasterCount
    band_type = dataset.GetRasterBand(1).DataType

    memDriver = gdal.GetDriverByName('MEM')
    memDS = memDriver.Create("TEMP", 256, 256, bands, band_type, [])
    memDS.SetGeoTransform( [94490, (256.0/(121069-94490)), 0, 96198, 0, 256.0/(96198-94490)] )

    for b in range(1, bands+1):
        s_band = dataset.GetRasterBand( b )
        t_band = memDS.GetRasterBand( b )
        data = s_band.ReadRaster( 94490, 69618, 121069, 96198, 256, 256, band_type )
        t_band.WriteRaster( 0, 0, 256, 256, data )

    createDriver = gdal.GetDriverByName('PNG')
    createDS = createDriver.CreateCopy("C:/xGIS/RichlandCounty/Aerial/SCRich09_1ft/3857/tiles/snap_2.png", memDS)
    createDS = None
    memDS = None

Does anyone have any insights or reactions that might help me get results quicker?


texts: QGIS Clipping Raster ERROR Cutline polygon is invalid

texts:
I want to clip a raster tif image based on a shapefile. However,the error as Fig.2 popuped when I use the Clipper Toolbox.What's the reason?I can make it in ArcGIS by Clip function in Toolbox.Thanks.

              Fig.1

             Fig.2


texts: Calculating unique ID value for a field in QGIS using python

texts:
I am trying to calculate a unique ID field for a vector layer in qgis, using python. I found some code to calculate attribute values, and tried to increment the value in the expression, but it doesn't work - it calculates all feature values as 0. I'm guessing the intial value I set for the expression is what it is reading each time - it doesn't seem to change it inside the loop. Anyone have any ideas? Here is my code:

QgsMapLayerRegistry.instance().addMapLayer(dsslvExplodeLyr)

idx = dsslvExplodeLyr.fieldNameIndex('UniqueID')

eVal = '0'
e = QgsExpression(eVal)
e.prepare(dsslvExplodeLyr.pendingFields())

dsslvExplodeLyr.startEditing()

for f in dsslvExplodeLyr.getFeatures():
    num = int(eVal)
    num = num + 1
    eVal = str(num)
    f[idx] = e.evaluate(f)
    dsslvExplodeLyr.updateFeature(f)

dsslvExplodeLyr.commitChanges()

I should mention I've also tried '$rownum' for the expression as that works in the field calculator gui but it also just gives a value of 0.


texts: Does hibernate spatial support st_geometry created by ArcSDE?

texts:
We are using ArcSDE 10.0 to export data to the Oracle database, and the spatial property are saved using the SDE.st_geometry format.

Now in our J2EE application we tried to manage these data by hibernate spatial.

And hs support Oracle Spatial out of box, but it have no documents about the st_geometry for ArcSDE.

So I wonder if we can use it for st_geometry ?  If not, is it possible to extend the functions of hs?

Update:

This is my dialect, and I register two functions overlaps and length:

public class Oracle11gSTGeometrySpatialDialect extends Oracle10gDialect {
    public Oracle11gSTGeometrySpatialDialect() {
        super();
        registerFunction("overlaps", new STFunction(new BooleanType(), "sde.st_overlaps(?1,sde.st_geometry(?2,3))"));
        registerFunction("length", new SQLFunctionTemplate(new DoubleType(), "sde.st_length(?1)"));
    }
}

And the query:

Query q=em.createQuery("from HBDS h where length(h.shape) > -1");

However the generated sql is invalid:

select hbds0_.objectId as objectId1_0_, st_astext(hbds0_.shape) as shape4_0_ from TEST_BDS hbds0_ where sde.st_length(st_astext(hbds0_.shape))>-1

While the valid sql should be:

select hbds0_.objectId as objectId1_0_, st_astext(hbds0_.shape) as shape4_0_ from TEST_BDS hbds0_ where sde.st_length(hbds0_.shape)>-1

It seems that this is caused by the ColumnTransformer annotation.

I am surprised that your example works as expected?

Now I can use the native query:

Query q = em.createNativeQuery("select st_astext(shape) from TEST_BDS h where sde.st_overlaps(h.shape, st_geometry(?,3)) = 1");
q.setParameter(1, JTS.defaultGeometryFactory.createPoint(new Coordinate(x, y)).toText());

It works, but there are some cons:

1 I have to write the select fields manually since I can not use select * which will be a problem once I have a lot of fileds in the entity.

2 I can not count the row count for the specified query.


texts: Displaying Feature Info for a user defined polygon

texts:
To retrieve feature info for a polygon selection, I have the OpenLayers code below:

        var featureTypes  = ["nig_admin1", "nig_admin2", "nig_admin3", "nig_landscan"];
        var populationCtrl = new OpenLayers.Control.GetFeature({
            protocol: OpenLayers.Protocol.WFS.fromWMSLayer(nig_landscan, {
                url: "http://localhost:8081/geoserver/wfs",
                featureType: featureTypes
            }),
            hover: false
        });

        populationCtrl.events.register("featureselected", this, function(e){
            select.addFeatures([e.feature]);
        });

I'd like to view the features retrieved in a pop up, or something similar. I've tried adding a pop up code within the featureselected event, but I can't get anything to show. Any pointers on how I can proceed will be highly appreciated.

nig_landscan in the code above refers to the corresponding WMS layer, which is a layer group in GeoServer.

My environment setup: GeoServer 2.1.3, PostGIS 1.5, OpenLayers 2.11.


texts: Calculating length downstream of edges in a network

texts:
I have a fairly complex and large river system that I would like to calculate length downstream between any two points within the river network. I have a calculation between any two points in the network in distance, but not downstream distance.

The river system was generated from aerial photography. It is connected as a network (it successfully produced distance between all points of interest in the network using network analyst). It also has river order. We are imposing flow on the river system now.

I made a trial geometric network with a simple river line I drew. I assigned a hydroID and generated From/To Nodes along the line. It populates the field "LENGTHDOWN" with 0s when I run "Calculate Length Downstream for Edges". Not sure what's up.

What is a workflow for being able to use the "calculate length downstream for edges" in the archydro toolbox using a river line generated not from a DEM?


texts: Sources for Geo-referenced data for Tiger populations in India or other countries

texts: I'm currently an Undergrad student in a GIS class, and I'm trying to conduct some research on tiger populations and urbanization. Does anyone know where I might be able to find data for tiger populations?


texts: How to select polygons containing points?

texts: I have a dataset that contains more than 10,000 polygon features. I also have a point layer with 300 points, some of which intersect the polygon features. I want to make a query and create a dataset of those polygons which contain at least one point. How can I do this?


texts: How to add my WFS data to ArcGIS operations dashboard?

texts:
I'm trying to build an operations dashboard in ArcGIS Online.
After adding my WFS layer to a web map I had to add some custom parameters to actually show the real-time data (points and metadata). When I added the map to my dashboard I wanted to add a chart, but the following notification came up:

<The map doesn't contain layers that can be used>.

Is there a way to fix this problem?
Parameters are srsNAME for dutch coördinate system and box to narrow down the amount of points.
Here's the WFS if it might help
https://waterwebservices.rijkswaterstaat.nl/services/distributielaagWFS/distributielaag_dbo?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&TYPENAME=locatiesmetlaatstewaarneming&srsName=EPSG:28992&bbox=140000,400000,180000,430000,urn:ogc:def:crs:EPSG:28992&Maxfeatures=50


texts: GeoServer; Blocking GeoServer-users' access to specific folders on host drives

texts:
Is it possible to adjust GeoServer in a way to block GeoServer-users' access to specific folders in the data directory? (I mean the actual folders on the computer drive where the shapefiles are physically stored)

I ask this question because of the following problem:

I have several users logging to the running GeoServer. When loading shapefiles (using "add store") they must not be able to access the shapefiles which do not belong to them. It seems that GeoServer do not work with the security rules of the host network (i.e. NTFS in this case).

I am not a security expert. I am familiar with the subject as much I usually need it.


texts:
I am trying to build an image pyramid data store in geoserver 2.14 running windows 10.

I have built the image pyramid as per these instructions: https://www.slideshare.net/geosolutions/creating-a-pyramid-with-gdal-reltile-for-serving-with-geoserver

My pyramid properties file looks like this:
Name=satellite

Levels=10.0,10.0 20.0,20.0
LevelsDirs = 0 1
LevelsNum=2
Envelope2D = 1084000.0,4722000.0 2092000.0,6217280.0

File structure is

c:\new\pyramid.properties

C:\new\0\0.dbf
C:\new\0\0.fix
C:\new\0\0.prj
C:\new\0\0.properties
C:\new\0\0.qix
C:\new\0\0.shp
C:\new\0\0.shx
C:\new\0\0.satellite_01_21.tif.....

c:\new\pyramid.properties
C:\new\1\1.dbf
C:\new\1\1.fix
C:\new\1\1.prj
C:\new\1\1.properties
C:\new\1\1.qix
C:\new\1\1.shp
C:\new\1\1.shx
C:\new\1\1.satellite_01_10.tif.....

Both folders 0 and 1, can read in as an imagemosaic, however I cannot get image pyramid to connect to anything

Everytime I hit connect I get the following error:

Could not list layers for this store, an error occurred retrieving them:
Failed to create reader from file://C:\new\pyramid.properties and hints
Hints: REPOSITORY = org.geoserver.catalog.CatalogRepository@12f7942
EXECUTOR_SERVICE = java.util.concurrent.ThreadPoolExecutor@182b120[Running,
pool size = 0, active threads = 0, queued tasks = 0, completed tasks = 0]
System defaults: GRID_COVERAGE_FACTORY = GridCoverageFactory TILE_ENCODING =
null FILTER_FACTORY = FilterFactoryImpl COMPARISON_TOLERANCE = 1.0E-8
FEATURE_FACTORY = org.geotools.feature.LenientFeatureFactoryImpl@1ec8458
STYLE_FACTORY = StyleFactoryImpl FORCE_LONGITUDE_FIRST_AXIS_ORDER = true
LENIENT_DATUM_SHIFT = true FORCE_AXIS_ORDER_HONORING = http

My research on this folder suggests a permissions issue with tomcat but I have been unable to find any installation of tomcat on my computer.


texts: Reprojecting EPSG:5514 to EPSG:3857 for use in Leaflet

texts:
I have one JPG (170 MB) and JGW file of a small area projected in EPSG:5514 (Image provider told me.)

JGW file:

                 0.1250000
                 0.0000000
                 0.0000000
                -0.1250000
           -584310.1875000
          -1213026.9375000

I am trying to create a custom map tiles for Leaflet in EPSG:3857.

I have tried to convert source JPG to TIFF

gdal_translate -of Gtiff -a_srs EPSG:5514 -co COMPRESS=JPEG -co JPEG_QUALITY=100 -co PHOTOMETRIC=YCBCR source.jpg output.tiff

Then I have used gdalwarp for reprojection

gdalwarp -of Gtiff -s_srs EPSG:5514 -t_srs EPSG:3857 -co COMPRESS=JPEG -co JPEG_QUALITY=100 -co PHOTOMETRIC=YCBCR output.tiff output_3857.tiff

Then I have transformed it to JPEG and changed the image size for zoom level

gdal_translate -of JPEG -outsize 16384 17700 output_3857.tiff output_zoom6.jpg

As a last step I generated the tiles

gdal2tiles -p raster -w all output_zoom6.jpg tiles

But if I use these tiles in Leaflet and I try to add a marker to the map, marker is in the wrong place...

Marker should be near the center of the image... If I use OpenStreet map tiles, marker is in the correct position, so my tiles are on the wrong place.

Any ideas what I am doing wrong?


texts:
What is the best way to access data I uploaded to mapbox with a click on tile.

I'm working on a web site using Node.js. I need to display a local map, click on the tile for a property and return some of the data in the db I uploaded to continue processing in my node web app.

I thought the best option was to use Mapbox.js to display the map but how do I get the data from the click? Or is there a better option then Mapbox.js.


texts: Flex Viewer: How do I return a subset of the input data using the geoprocessing widget?

texts:
I have a set of point data (representing crashes) that has been published as a feature service and which I have pulled into a Flex (3.5) application. I want the users to be able to submit a query and create a subset of the original data - for example, to access all the crashes that occurred on a Tuesday, or which were fatal, etc.

To do this I have created a Geoprocessing widget based on a very simple model (Get expression, Select by attribute; copy features) that has been published as a GP service.

However, I haven't been able to get the GP widget to work properly. Once run, it adds the output 'Query' layer to the layer list, but there are no features within that layer (ie, nothing on the map).

Geoprocessing Widget config file:

<?xml version="1.0"?>
<configuration label="Select by Attribute">
    <taskurl>http://ags-dev2:6080/arcgis/rest/services/CAS_Crashes/SelbyAttr/GPServer/SelbyAttr</taskurl>
    <shareresults>true</shareresults>
    <inputparams>
        <param name="Expression"
               label="Expression"
               required="false"
               type="string"
               visible="true">
        </param>
    </inputparams>
    <outputparams>
        <param name="Output"
               label="Query"
               type="featurerecordset"
               visible="true"
               geometrytype="point">
            <simplerenderer>
                <simplemarkersymbol color="0x2D9C49" width="4"/>
            </simplerenderer>
        </param>
    </outputparams>
</configuration>

Does anyone have any ideas about what I might need to change to bring the features through to the output layer?

I'm using ArcGIS Desktop 10.2, ArcGIS Server 10.2, Flex 3.5.

Thanks!


texts:
I wanted to use the r.viewshed tool in QGIS. Everytime I want to start the tool there is a error (german):

Dieser Algorithmus kann nicht ausgeführt werden run :-(
GRASS GIS 7 ist anscheinend nicht richtig auf dem System installiert und konfiguriert. Bitte vor der Ausführung von GRASS GIS 7-Algorithmen installieren.

Do I have to install GRASS 7 to use the tools?

I am using a Mac with 10.13 and QGIS 2.18


texts: How do PostGIS functions handle buffers and external storage?

texts:
I see an ever increasing list of new functions in PostGIS, some of which involve GEOS (e.g. ST_ClusterKMeans). Some functions (e.g. those in pgrouting) rely on other libraries (e.g. BGL).

My impression is that many of these underlying libraries (often in C/C++) do not handle buffers management between the memory and secondary memories/storage/disks.

So do the PostGIS functions on top of them work on large data sets that cannot be stored in the physical (or virtual) memory?

If so, where do these buffer management capabilities come from (from the point view of implementation)?


texts: Calculating length of road network within a polygon

texts:
I have a road network stored in a gdb and some polygons in another gdb.

I want to populate an attribute for each polygon with the total length (in km) of the road network that it overlaps.

Any idea on how to do that?


texts: Why does WGS84 refer to so many different things?

texts:
Sounds like a stupid question, but I've been struggling with this.

Here's a diagram of the structure for EPSG 4326 - the geographic CRS known as WGS 84.

Here we have WGS 84 (a CRS), World Geodetic System 1984 (a datum) and WGS 84 (an ellipsoid). As a newcomer to GIS this is very confusing, and I'm wondering why it was named like this.

So what's the most common usage of the term WGS 84 - the ellipsoid, the datum, or the CRS?


texts: How to convert a bil image into grid in ArcGIS?

texts:
I've downloaded RFE data from FEWSNET. It's a compressed file which includes .bil files.

Now i want to convert this ".bil image" into a grid to use the precipitation data in an area.

I'm using ArcGIS.

How do i do it?


texts: GEOS and Shapely installation on Windows

texts:
When trying to install Shapely on my Windows 64bit computer, I cannot get the GEOS library to work.

So far, I have run the OSGeo4W installer from which I installed GDAL (I believe the geos library is included in that package). After that, I checked and I have geos_c.dll on my C:\OSGeo4W\bin directory, but either I have missed some configuration steps or the library does not work.

I need Shapely to work, so I also ran pip install shapely after  installing GDAL, and it apparently worked (although it could not find the C library for GEOS).

In my code, I can import Shapely, but when I try to use it, I get an error telling me "geos.dll" is not found.


texts: OL3 How to select features of a WMS layer

texts:
I want to show on the map the features of the saved layers from a GeoServer server, but I don't know how to do this...

Like this, but with the dragbox interaction...

I have this code:

var select = new ol.interaction.Select();
mapa.addInteraction(select);
var selectedFeatures = select.getFeatures();
var dragBox = new ol.interaction.DragBox({
        condition: ol.events.condition.always,
        style: new ol.style.Style({
            stroke: new ol.style.Stroke({
                color: [0, 0, 255, 1]
            })
        })
    });

    mapa.addInteraction(dragBox);

    dragBox.on('boxend', function(e) {
        if (capaActiva != null) {
            var info = [];
            var extent = dragBox.getGeometry().getExtent();
            var source = capaActiva.getSource();
            console.log(source);
            var url = source.getGetFeatureInfoUrl(extent, viewResolution, 'EPSG:3857',{'INFO_FORMAT': 'text/plain'});
            console.log(url);

        } else {
            Ext.MessageBox.show({
                title: 'Error',
                msg: 'Check a layer',
                buttonText: {
                    ok: 'Aceptar',
                }
            });
        }
    });


texts: Find the closest line vertex from the linestring using PostGIS

texts: How to find the closest line vertex from a linestring? ST_ClosestPoint will give a point which is not from the linestring.


texts: Are Rwanda and Burundi part of "Eastern Africa countries" or "East African countries". Some sources are differing.


texts: QGIS draw line segments with "geometry generator" based on segments angle

texts: I have quite a hard problem to solve here and I don't know if its even possible with the geometry functions in QGIS. I have a line geometry with equidistant vertices, that is, same length segments. Based on a treshhold angle of a segment and its subsequent segment, I'd like to generate new offset line segments - if the angle is > 10 but < 180, I would have all the 'hard' right turns and would like the newly created line segments for these two segments to be offset to the left, if the angle is > 180 and < 350 then I'd have the 'hard' left turns and would like the newly created segments offset to the right side. That is, the segments with 'hard turns' would have a offset line segment at the outer side of the original two line segments. Like in the example picture below! The math for the angle, I think, are in this example: https://stackoverflow.com/questions/4294638/how-to-calculate-angle-between-two-direction-vectors-that-form-a-closed-open-sha


texts: Address locator not using modified .cls file

texts:
I have made a customization to the US Streets locator style for our application.  We've changed "SR" to map to "STATE ROAD" instead of "HC" in the us_addr.cls file in the C:\Program Files (x86)\ArcGIS directory.  When I run StanEdit, it maps "SR" correctly.  When I do a test in ArcMap using the Geocoding Standardize Addresses tool, picking the the US Streets with City, State, Zip and AltName, it works correctly also.

I've created a locator in ArcCatalog using the US Streets with City, State, Zip and AltName style and registered it on my local server (which is the same machine I did the above testing on).  When I run the batch geocoding tool using that locator or use the Find Address tool, it maps "SR" to "HC".

So I'm thinking that I'm not editing the correct us_addr.cls file. But I have no idea where the locator is pulling from. Is there any documentation anywhere (besides the Geocoding Rule Base Developer Guide - been all thru that) that explains how to modify the locator syles?

We are running ArcGIS 9.3.1 SP 2.


texts:
I have a problem with setting up a qgis web client on remote server. I set it up on localhost and it runs smoothly, but when I done exactly the same on remote server it is just not working, even a default maps. Helloworld example is not showing up layers, and my project is working only when i exclude about 20 layers, and by the way all styles are crashed. I could not find an issue, and exacly the same projects and many more runs great on localhost.

It is link to remote server: 185.25.148.67/qgis-web-client/site/index.html

EDIT:

I noticed that on localhost server sends POST requests and on remote server it is always a GET request. I think that could be the issue, since request to display all layers in GET method would be too long. But i still couldn't figure out how to enforce POST method.

2nd EDIT:

It is POST on remote server at the end, but it still doesn't work, here is link to the project with more layers:

185.25.148.67/qgis-web-client/site/qgiswebclient.html?map=/home/web/qgis-web-client/projects/rzgow/start.qgs


texts: Is point on land (Natural Earth shapefile)

texts:
I am attempting to use the second option in Check if a point falls within a multipolygon with Python to check if a given point is on land or water.  The shape file is from Natural Earth: https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-land/

Here is my code:

#!/usr/bin/python3

import shapefile
from shapely.geometry import Point # Point class
from shapely.geometry import shape # shape() is a function to convert geo objects through the interface

#point = (-80,40) # an x,y tuple
point = (40,-80)
myshp = open('./natural_earth/ne_10m_coastline.shp','rb')
mydbf = open('./natural_earth/ne_10m_coastline.dbf','rb')
shp = shapefile.Reader(shp=myshp, dbf=mydbf) #open the shapefile
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()
for i in range(len(all_shapes)):
    boundary = all_shapes[i] # get a boundary polygon
    if Point(point).within(shape(boundary)): # see if the point is in the polygon
        print('land')
       myshp.close()
mydbf.close()

I used this code on the USG Tiger data set to tell what state an county a point was in.
I have tried reversing the order of the coordinates.  This code never finds the point over land. Can someone diagnose this?


texts:
I have created an ImageMosaic layer with GeoServer. In the application that uses GeoServer's WMS the layer has whitespaces at the edges of .sid as seen:

However when using GeoServer's layer preview I see that the layer is fine.
Update. In layer preview I was viewing not the same layer as in first image. In layer preview was the layer that was made from .tif files (I used gdaltranslate on .sid files). So when attempting to see a layer preview (of same .sid layer as in first image) I see white image only! I am leaving second image for reference.

Note: black lines - .sid edges, blue square - area seen in first image, red - area where whitespace is

Query of a layer preview request:

service=WMS&VERSION=1.1.0&REQUEST=GetMap&layers=X_AROUND_75_25_2_C&FORMAT=image%2Fpng&STYLES=&SRS=EPSG:3346&bbox=575000.0,6025000.0,585000.0,6035000.0&width=768&height=768

Query of one of GetMap requests that browser makes:

SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=X_AROUND_75_25_2_C&FORMAT=image%2Fpng8&STYLES=&SRS=EPSG%3A3346&BBOX=581202.1631283718%2C6032956.741882952%2C582556.8325043772%2C6034311.411258957&WIDTH=512&HEIGHT=512&TRANSPARENT=true&TILED=true


texts: Rule-Based: filter test is not working on QGIS 3.4

texts:
I made one Atlas Map based on this exercise (https://anitagraser.com/2016/10/15/qgis-atlas-tutorial-material-design/) and I had good results working with QGIS 2.18.

But when I use QGIS 3.4 I can not get the same results.

As you can see in this picture bellow in QGIS 2.18 I have this notification when I run the filter test:

This is what happens when I run the same filter code in QGIS 3.4:

  Filter expression parsing error:

  $atlasfeatureid function is not known

  $map function is not known

So in this case - What do I need to do?


texts: How can I use the SFCGAL version of ST_Union?

texts:
I'm trying to combine a bunch of geometry boundaries to get the linework of their borders, but I'm getting an error with the default, GEOS-based ST_Union:

SELECT ST_Union(ST_Boundary(geom)) FROM my_table;

gives

ERROR: GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING (391229 4.79279e+06, 391444 4.79279e+06) and LINESTRING (391444 4.79279e+06, 391229 4.79279e+06) at 391437.43387050979 4792789.729486906

According to the PostGIS documentation, SFCGAL has a version of ST_Union that I can use:

  Some SFCGAL functions replace standard ones (ST_Intersects, ST_Intersection, ST_Difference, ST_Union, ST_Area and ST_Distance), to switch between standard functions and SFCGAL function use:

SET postgis.backend = sfcgal;

However, this has no effect:

SET postgis.backend = sfcgal;
SHOW postgis.backend;
SELECT ST_Union(ST_Boundary(geom)) FROM my_table;

0 Rows Affected

postgis_backend
----------------
sfcgal

ERROR:  GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING (391229 4.79279e+06, 391444 4.79279e+06) and LINESTRING (391444 4.79279e+06, 391229 4.79279e+06) at 391437.43387050979 4792789.729486906

Based on the fact I'm still getting a GEOS error, I'm presuming it's not actually using the SFCGAL function.

The array version of ST_Union produces the same result:

SELECT ST_Union(ARRAY_AGG(ST_Boundary(geom))) FROM my_table;

How can I get PostGIS to invoke the SFCGAL version of ST_Union?

Note that I do have SFCGAL installed:

SELECT postgis_full_version();

POSTGIS="2.5.1 r17027" [EXTENSION] PGSQL="100" GEOS="3.7.0-CAPI-1.11.0 673b9939" SFCGAL="1.2.2" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 1.11.4, released 2016/01/25" LIBXML="2.9.1" LIBJSON="0.11" RASTER

SELECT * FROM pg_catalog.pg_extension WHERE extname LIKE '%sfcgal%';

extname        | extversion | extconfig
---------------+------------+-----------
postgis_sfcgal | 2.5.1      | NULL

PostgreSQL Version:

SELECT version();

PostgreSQL 10.6 on x86_64-pc-linux-gnu, compiled by gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-28), 64-bit

 Aside: I don't know if the SFCGAL function will help with my error, and that's not really the topic of this question. I just want to try using to see if it does or not.


texts: Where does QGIS save the processing log for a specific project?

texts:
After updating to current stable version (3.4.12 Madeira at time of writing), when opening a project, the processing log as described in https://docs.qgis.org/3.4/en/docs/training_manual/processing/log.html is not accessible any more through the previously visible "Processing" tab. I would like to access the processing history of a specific project to review it.

This already happened to me when I updated from 3.4.1 to 3.4.2, although I wasn't certain about it back then, because I hadn't paid attention specifically to it, since I didn't expect to lose the processing history.

(I handed in my project with an incomplete history file. Nobody noticed because they didn't run it, but, you know, it's bad and not contributing to reproducible science. For my next project, if I update to a newer version of QGIS, I'd like to be able to provide a full processing history :D )


texts: Programmatically export to geotiff in different scales from ArcGIS Desktop?

texts:
I would like to programatically export geotiff files in different scales (e.g. 1:25000, 1:50000 etc) using arcgis10.

As input I have some shape files and a DEM file as well as a corresponding style.

I tried to understand the different possibilities by reading the online-documentation but I am unsure what Api and Calls to use.
C#, C++ and python are all ok but I would like to create a fully automated solution. I think a console application might be the best fit but executing a python script from within arcgis might also be ok.

I have the following questions:

Should I use the .Net API or is a python script sufficient in this scenario?
What are the commands / methods I have to call? A code example would be great but pointers to the right documentation is also nice.


texts: ArcGIS 10.0/arcpy: How to create a polygon geometry from inner and outer ring point arrays?

texts:
I have two arrays, each containing a closed ring of points, one for an outer ring and the other for an inner ring, the latter running counterclockwise, of course.  I am using the arcpy.Polygon() factory function to build the polygon.  My question:  how do I input the points to Polygon()?

Documentation gives the general syntax as Polygon(inputs, spatial reference, has_Z, has_M).  The inputs part is what is sketchy.  The inputs supposedly can be point objects or arrays.  My problem is that I have a polygon with a doughnut hole.  I can build a simple polygon with a single point array as input ... no problem; however, in what way do I add the two sets of points?  In the "Reading Geometries" help page it says that parts (in this case rings) are separated by a null point object.  For polygons with holes, the outer ring goes first.  So, I tried adding the outer array then the inner array to a container array.  That didn't work:  no hole.  I tried the same thing with a null point object between them.  No joy:  same result.  Then, I tried putting both arrays into a Python list as the input.  That threw an error.  The sample scripts in the help files never deal with writing polygons with holes.  So, I really don't know what format to use for the inputs.  Does anyone have any insight into this one?  I would surely appreciate the help.  Thanks, all.


texts:
I am using QGIS 3.10 and would like to visualise my data using a custom svg whose source code is:
<svg enable-background="new 0 0 175 375" height="375" viewBox="0 0 175 375" width="175" xmlns="http://www.w3.org/2000/svg"> <path d="m 40.445708,15.5 c 0,0 30,85 -15,165 -36.897041,65.59474 -20.0000005,115 -10,135 C 44.514588,373.63777 150,365 170.4457,305.5 190.99887,245.68723 189.06791,170.62949 40.445708,15.5 Z" fill="param(fill)" fill-opacity="param(fill-opacity)" stroke="param(outline)" stroke-opacity="param(outline-opacity)" stroke-width="param(outline-width)"/> </svg>.
In QGIS the .svg is displayed correctly and can be adjusted as I wish:

In the print composer, this .svg is distorted into a square:

How can I prevent this?


texts: Problems using the 'Split' tool with ArcGIS

texts:
I'm running into difficulty when using the split tool with ArcGIS 10.

I am trying to split lines using a grid that I generated using the Fishnet command. Each gridcell has a unique text value as a name (FID converted to a string), which is my Split Field. However, I can't tell whether it is the size of my grid (300000 cells), or if it is due to no line being within a cell that causes the error.

I tried breaking my grid cell into smaller pieces of ~50000 cells, but I still got an error occuring. Then I tried selecting only cells that were intersected by the lines, and just splitting using those cells, and still got an error. The strange part about the error messages, were that the analysis would run for about 30 minutes or so, producing output, and then an error message would appear on screen, but no explanation - no error number or no further details.

What I have no resorted to doing is splitting the grid into 500 smaller grids, and then splitting each subgrid on itself. Then I will loop through every cell, and clip the lines to each cell. The problem with this approach is that it is taking a long time (~20 hours so far).

Are there any ways I can improve this process?


texts: ArcMap hangs after clicking on print.  I get the dreaded spining globe icon as if it is working, but does nothing.  I eventually go to task manager to end the task.


texts: How do I let the links of the OpenLayers attributions open up in a new tab?

texts:
Our webserver has set X-Frame-Options to 'sameorigin'. So if I want to open the attribution links in an OpenStreetMap, the requests get blocked. So I want to let the links be opened in a new tab. How do I change the target-property of the a-elements to '_blank' while still using the OpenLayers API?
I want do change it with the OpenLayers API and not "from outside" with JQuery etc.
enter image description here


texts: CRS with Latitude and Longitude the same way round as Google

texts:
I have a set of data with Lat/Lon co-ordinates taken from Google Earth, I've converted the co-ordinates into WKT using a short formula in a spreadsheet, converted to a csv and imported. Unfortunately Google presents its points as Lat/Lon whereas all the CRSs I usually use only seem to understand Lon/Lat. I thought there were conventions for this sort of thing so that GIS systems could talk to one another!

Anyway, does anyone know of a CRS that recognises Lat/Lon, in that order?

Thought it might be useful to give an example - co-ordinates come to me looking like this "50.698162325700466,-4.674852266907692, 50.69849212863936,-4.674183055758476" my formula converts them to this "MULTIPOINT ((50.698162325700466 -4.674852266907692) (50.69849212863936 -4.674183055758476)). The points are all in England, so you can see they are in Lat/Lon.
I've just thought that +axis=neu should work as a custom CRS, but it doen't seem to be translating correctly. It's turning the points around as required, but also shifting them about.


texts: Getting feature class names within dataset with vb.net and ArcObjects?

texts:
I am totally new to arcobjects. I want to get the feature class names withın a dataset and add them to a listbox.  I was able to get the dataset names but i do not know how to iterate through the feature class names within these datasets.

Public Sub GetWorkspace()

    Dim pEnumGxObj As ESRI.ArcGIS.Catalog.IEnumGxObject
    Dim pFCFilter As New ESRI.ArcGIS.Catalog.GxFilterWorkspaces
    Dim pGxDatabase As ESRI.ArcGIS.Catalog.IGxDatabase2
    Dim pGxDialog As New ESRI.ArcGIS.CatalogUI.GxDialog
    Dim pGxObject As ESRI.ArcGIS.Catalog.IGxObject
    Dim pWorkspace As ESRI.ArcGIS.Geodatabase.IWorkspace

    pGxDialog.ObjectFilter = pFCFilter

    If Not pGxDialog.DoModalOpen(My.ArcMap.Application.hWnd, pEnumGxObj) Then Exit Sub

    pGxObject = pEnumGxObj.Next

    If pGxObject Is Nothing Then Exit Sub
    Try
        If TypeOf pGxObject Is ESRI.ArcGIS.Catalog.IGxDatabase2 Then
            pGxDatabase = pGxObject
            pWorkspace = pGxDatabase.Workspace
            System.Windows.Forms.MessageBox.Show(pWorkspace.PathName)
        Else
            System.Windows.Forms.MessageBox.Show("This is likely a shapefile workspace")
        End If

        Dim pEnumDSName = pWorkspace.DatasetNames(esriDatasetType.esriDTAny)

        Dim pSdeDSName = pEnumDSName.Next
        While Not pSdeDSName Is Nothing
            ListBox1.Items.Add(pSdeDSName.Name)
            pSdeDSName = pEnumDSName.Next
        End While

    Catch ex As Exception
        System.Windows.Forms.MessageBox.Show(ex.ToString, "GetWorkspace")

    End Try
End Sub


texts: Automated spline polygons from four points in QGIS?

texts:
I'm looking for a way to use QGIS for tree surveys which involves plotting four points to demarcate the NESW extents of a tree canopy and then generating a polygon that passes through each of these - I know that AutoCAD offers this functionality but I was hoping to achieve it through QGIS. I can create a convex hull easily enough but this is only straight lines connecting the points.

Is there a way to generate a curved polygon (or line) that could be automated over large datasets?


texts: List of data file formats supported in QGIS including read/write support

texts:
Is there a list of data file formats supported in QGIS, including information about read/write support? What I'm looking for is something like:

Geopackage: r/w
Esri Shapefile: r/w
CSV: r
KML: r
etc.

The documentation has a section "Exploring Data Formats and Fields".
However, there is no simple overview over read/write support and you have to read through a lot of text that does not always mention explicitly which of the formats are editable.
In fact, the vector layer export dialog offers a list of file formats to which you can save, but a) this includes only writable formats (not readable as well) and b) it is not available outside of QGIS as a comprehensive list.


texts:
I need help converting a download from Geofabrik to Illustrator. It had worked for me earlier, but I can't remember what I did, and now it doesn't.

I downloaded a country from Geofabrik (because city wasn't available) and now I want to extract three different cities (roads only) for an art project
its an 11GB .osm file
I Go to: Vector (menu) → OpenStreetMap → Import Topology from XML.
Then I get a .osm.db file
Then I do Vector (menu) → OpenStreetMap → Export Topology to SpatiaLite.... (and select the .osm.db file)
I select poly lines and my highway/road/lane tags
But at this state QGIS goes into a freeze (I have 12GB RAM).
What now?

Will try those options.
Meanwhile, the export seems to have happened, but the Map Canvas is blank. Anything I should know here?

My newly exported file is a new layer, has 1943641 feature counts etc, but MapCanvas is blank. How do I see the map, zoom to my area, and create a New Print Composer?


texts: Setting up leaflet map in R using if/else if/else filtering

texts:
I'm trying to put together a map for work. Right now, the script looks as such (The shape file data can be downloaded here):

library(tigris)
library(mapview)
library(leaflet)
library(sf)
library(dplyr)
library(shiny)
library(shinyWidgets)

nv <- st_read('./testdata.shp') %>% mutate(super_region = case_when(region %in% c('1','2') ~ 'A', !region %in% c('1','2') ~ 'B')) %>% mutate(turf = as.character(turf)) %>% select(super_region,region,turf,areas,sales) %>% st_sf()

ui <- fluidPage(
  fluidRow(
    column(8,
           width = 10, offset = 1,
           tags$h3("Select Area"),
           panel(
             selectizeGroupUI(
               id = "filters",
               params = list(
                 SR = list(inputId = "super_region", title = "Super Region:"),
                 Reg = list(inputId = "region", title = "Region:"),
                 Turf = list(inputId = "turf", title = "Turf"),
                 Areas = list(inputId = "areas", title = "Areas:")
               ))
           ),
           leafletOutput("test")
    )
  )
)

server <- function(input, output, session) {

  res_mod <- callModule(
    module = selectizeGroupServer,
    id = "filters",
    data = nv,
    vars = c('super_region', 'region', 'turf', 'areas')
  )

map_data <- reactive({

if(is.null(input$super_region) & is.null(input$region) & is.null(input$turf) & is.null(input$areas)) {

  md <- res_mod() %>% st_sf() %>% st_buffer(0) %>% group_by(super_region) %>% summarize(geometry = st_union(geometry))

}

else if(!is.null(input$super_region)) {

    md <- res_mod() %>% filter(super_region %in% input$super_region) %>% st_sf() %>% st_buffer(0) %>% group_by(region) %>% summarize(geometry = st_union(geometry))

}

else if(!is.null(input$region)) {

    md <- res_mod() %>% filter(region %in% input$region) %>% st_sf() %>% st_buffer(0) %>% group_by(turf) %>% summarize(geometry = st_union(geometry))

}

  else if(!is.null(input$turf)) {

    md <- res_mod() %>% filter(turf %in% input$turf) %>% st_sf() %>% st_buffer(0) %>% group_by(areas) %>% summarize(geometry = st_union(geometry))

  }

else {

  md <- res_mod() %>% filter(areas %in% input$areas) %>% st_sf() %>% st_buffer(0) %>% group_by(areas) %>% summarize(geometry = st_union(geometry))

}

})

output$test <- renderLeaflet({

  res <- map_data()

  mapview(res)@map

})

}

shinyApp(ui,server)

While I thought I had set it up correctly, it is not working as intended, and it appears my issue is with the combination if, else if and else statements I put together, but I'm not sure what I'm doing wrong. Here is what I would expect it to look like:

When the map first loads and no filters are selected:

When only super_region is selected. In this example, I used super_region B

When region is then selected (super_region can be blank or have a value). In this example, I used region 3:

**When turf is then selected (super_region, and region can be blank or have a value, it doesn't matter). In this example, I used turf 3A:

**And finally, when areas is selected (super_region, region, and turf can be blank or have a value, it doesn't matter. In this example, I used areas 1167 is selected:

Does this make sense? So basically, the default setup is showing just the super_regions. Then for each ensuing level, it groups by the one below it.


texts: Vector Layers are not lining up and are scaled differently

texts: Hey I am a first time user of GIS and I am working for a transportation map of North Carolina for a project at my work. So what I have been doing is first placing the Tiger County Boundary by using the add vector layer tool and selecting the directory setting and then select tiger. I then do the same thing for a layer that contains all of the roads. They do not line up and the road layer is significantly larger than the county boundaries. I tried to mess with the CRS, but it did not seem to change anything. I am just taking a guess, but I am assuming that it is a scaling issue? Anyway please help


texts: Pan sharpening with surface reflectance data

texts:
The surface reflectance data for landsat 8 is only available for bands 1-7 (via EPSA bulk download). I want to pan sharpen surface reflectance data. But I don't have a surface reflectance band 8 (pan band).

Do you think I can use for pan sharpening the surface reflectance data (band 1-7) with the raw DN pan band, as I don't have a surfcace reflectance panchromatic band?