texts: Python: shapefile.Reader, how to set up .iloc?

texts:
I would like to set up range of reading by shapefile.Reader. I couldn't find any suitable documentation.
I would like to read rows 0-100. In Pandas or GeoPandas I use df.iloc[0:100], but what should I use here?
Code:
with shapefile.Reader(shapefile_path) as shp:
    total_rows = shp.numRecords
    for row_num, row in enumerate(shp.iterRecords()):
        print(row)


texts:
I would like to use Java GDAL bindings to subset a GeoTIFF and write the result in VRT format.

I have found several old posts, one of them states that this functionality is yet to be implemented in Java GDAL bindings:

http://lists.osgeo.org/pipermail/gdal-dev/2014-June/039337.html
https://stackoverflow.com/questions/25522833/how-to-write-a-output-file-in-java-using-gdal-library

What is the situation in 2017?

Did anyone try to save VRT from Java using GDAL bindings?


texts: Is there a web service to get the area of a city?

texts:
I'm looking for an api/service/db that contains the surface area of every city in the world.

For example, http://en.wikipedia.org/wiki/London, the box on the right contains that information.

I've searched but was unable to found anything useful. Google doesn't seems to have any API for that and I don't want to parse Wikipedia. Even paid services are fine!


texts: Converting ArcGIS date value to string reading in mm/dd/yyyy format using ArcPy?

texts:
I am writing an arcpy tool script in ArcGIS 10.0.  In a fGDB feature class, I have a Date type field.  I gather that the values in the Date field are stored as long integer dates, but they display in mm/dd/yyy format in the table.  What I want to do is capture the integer date value as a string showing the date in mm/dd/yyyy format.  How do I go about this?

There are many examples of how to convert a date/time string to a date object/value but not vice versa.  I am appending features from a feature class with many fields, including a date field, into a standard-schema fc (which I can't change) that has few fields.  To save a lot of the technical attributes in the input, I must dump it all into a catch-all, user-controlled text field that is defined in the target.  The date data from the input is one of those pieces of data I need to save, hence the need to convert the date to a string.

There is likely a fairly easy solution I'm missing, but at the moment I am at a loss.  I appreciate any help that comes my way.


texts: Writing list of feature classes with fields to one csv file per feature class?

texts:
I have a database with multiple feature classes. Now I extracted a list of feature classes and a list of fields. I figured out how to write the list of feature classes and the corresponding fields to csv. But the csv file gets really long. Instead I would like to write each feature class name with the corresponding field names to a separate csv file.
My code looks like this right now:

def csvwriter (csvfile, mylist):
    myfile = open(csvfile,'ab')
    wr = csv.writer(myfile)
    wr.writerow(mylist)

datasets = arcpy.ListDatasets()
for ds in datasets:
    dslist = []
    dslist.append(ds)
    csvwriter (csvfile1,dslist)
    for fc in arcpy.ListFeatureClasses("*","",ds):
        fclist = []
        fclist.append(fc)
        csvwriter(csvfile1,fclist)
        for field in arcpy.ListFields(fc):
            fieldlist = []
            fieldlist.append(field.name)
            csvwriter(csvfile1,fieldlist)


texts: Error occured decoding the epsg code - 3857 (Web Mercator)

texts:
I'm trying to add the following wms service to QGIS using EPSG:3857 and get the following error:

Service (running on L/A/P/P geoserver):

http://66.228.57.208/geoserver/wms?tiled=true or non-cached http://66.228.57.208/geoserver/wms

Layer:

Carto_Basemap_WGS84WM layer group (presented in 3857 which contains ~ 8 layers all in 3857)

Error:

Response: <?xml version="1.0" encoding="UTF-8"?><ServiceExceptionReport     version="1.3.0" xmlns="http://www.opengis.net/ogc"     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/ogc http://66.228.57.208:80/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd">   <ServiceException code="InvalidCRS">
  Error occurred decoding the espg code urn:x-ogc:def:crs:EPSG:3857
Database failure while creating a &apos;CoordinateReferenceSystem&apos; object for code &quot;EPSG:3857&quot;.
Connection is closed

Help!


texts: Is there any method to create WMTS GetCapabilities response in xml using java libraries?

texts: I am creating WMTS server side. The problem is that I cannot find any solution where I will be able to generate GetCapabilities xml response in java using some libraries like geotools or similar. Can somebody help me with this?


texts: Error exporting .geojson file with polygon, line and point data. Attached is the screenshot of the error for reference. Can anyone suggest solution for this kind of error


texts: Plot a world heat map with mollweide projection in R

texts:
I have a data set with global observations with lat/lon information. There are 288 longitudes from -179.375 to 179.375 with spacing 1.25, 180 latitudes from -89.5 to 89.5 with spacing 1. There is one observation at each lat/lon combination. The total number of observations is 288*180=51840.

I want to make a heat map with world map embedded as well. Here is the code to make a rectangular map.

library(fields)
library(maps)
library(maptools)

set.seed(1)
unique.lat <- seq(-89.5, 89.5, by = 1)
unique.lon <- seq(-179.375, 179.375, by = 1.25)
lat <- rep(unique.lat, each = length(unique.lon))
lon <- rep(unique.lon, length(unique.lat))
y <- lat + abs(lon) + rnorm(length(lat))  # generate some data for plotting

dat <- data.frame(lat = lat, lon = lon, y = y)
data.mat <- matrix(dat$y, nrow = length(unique.lon), ncol = length(unique.lat), byrow = F)
surface <- list(x = unique.lon, y = unique.lat, z = data.mat)
plot.surface(surface, type = "I", axes = FALSE, xlab = "Longitude", ylab = "Latitude", legend.cex = 0.8, legend.mar = 7)
map("world", add = T)
map.axes()

What I really want is a heat map using mollweide projection. I learned how to make a simple world map with

data(wrld_simpl)
plot(spTransform(wrld_simpl, CRS = ("+proj=moll +lon_0=0 +x_0=0 +y_0=0")))

How do I somehow combine these two plots and make a heat map on the mollweide-projected map?


texts:
In geoserver only geometry type of SQL server table is detected. Which means SLD style can't be applied ogc filter tag. But none of my SQL server field beside geometry are detected by geoserver.
How can I do that ?

Create table script

CREATE TABLE My_GIS_Table
(
 LOCALNAME varchar(50),
 Coord  geometry not null,

 constraint PK_My_GIS_Table primary key (LOCALNAME)
)

SLD created with uDig from a shapefile.

<sld:Rule>
                    <sld:Name>MY_CITY_NAME</sld:Name>
                    <ogc:Filter>
                        <ogc:PropertyIsEqualTo>
                            <ogc:PropertyName>LOCALNAME</ogc:PropertyName>
                            <ogc:Literal>MY_CITY_NAME</ogc:Literal>
                        </ogc:PropertyIsEqualTo>
                    </ogc:Filter>
                    <sld:PolygonSymbolizer>
                        <sld:Fill>
                            <sld:CssParameter name="fill">#E31A1C</sld:CssParameter>
                            <sld:CssParameter name="fill-opacity">0.45</sld:CssParameter>
                        </sld:Fill>
                        <sld:Stroke>
                            <sld:CssParameter name="stroke-dashoffset">5</sld:CssParameter>
                            <sld:CssParameter name="stroke-dasharray">2.0</sld:CssParameter>
                        </sld:Stroke>
                    </sld:PolygonSymbolizer>
</sld:Rule>

Hypothetic SLD I'd like to create for my SQL Server table

Here I just want to change the part with  tags. With something posibly like  tags. Ogc tags and Geoserver wont recognize MY_GIS_TABLE field name LOCALNAME  so I don't know how to edit the existing SLD to create a filter upon that field.

<MSSQL:Filter>
    <MSSQL:PropertyIsEqualTo>
    <MSSQL:PropertyName>LOCALNAME</MSSQL:PropertyName>
    <MSSQL:Literal>MY_CITY_NAME</MSSQL:Literal>
    </MSSQL:PropertyIsEqualTo>
</MSSQL:Filter>


texts:
I'm new to the GIS world, and I'm trying to work out what tool I should start with in order to analyse some data and run custom spatial queries on it.
I have a data set representing the recording of a construction site vehicle and its activity over a period of time. Each record in the data contains such fields as lat/long, timestamp, and various parameters describing what activity the vehicle was performing at that particular point in time.
Without spending weeks developing my own tool, I would like to be able to display a map, select a start and end time (to filter out the relevant records), and render a trace of where the vehicle went and somehow indicate what it was doing where. Eg, it may be a grader machine, and I simply want to display whether the blade was engaged to the ground or not (on/off). Maybe just green points where it was engaged, and red points where it was not.
I'd also like to be able to run some spatial queries, eg within a certain bounding box, how much time was spent with the blade engaged to the ground. Or how far did the vehicle travel during a certain time period, and for how much of that distance was the blade engaged. Simple queries like that.
A bonus would be if I could achieve an animation of the activity, ie a play-back of the vehicle moving and plotting its route and activity, but this is not necessary.
Could I achieve this with something like QGIS, or what tool would you recommend, if any?


texts: Shapefile not displaying as a layer on leaflet map

texts:
I've been trying to learn how to display my shapefiles as a layer on leaflet map then i came across calvinmetcalf github . Using his example i was able to display the councils.zip file as a shapefile layer on the map. But, trying to the do the same with my own shapefile produce no result as the shapefile refuse to load on the map. I couldn't see an error in the console either, so i am stuck on what i could be doing wrong. Here is the code below.

    <script>
    var map = L.map('map',
    {
        center:[11.47,8.20],
        zoom: 10
    });
    L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png').addTo(map);
    var shpfile = new L.Shapefile('lg_polygon2.zip', {
        onEachFeature: function(feature, layer) {
            if (feature.properties) {
                layer.bindPopup(Object.keys(feature.properties).map(function(k) {
                    return k + ": " + feature.properties[k];
                }).join("<br />"), {
                    maxHeight: 200
                });
            }
        }
    });
    shpfile.addTo(map);
    shpfile.once("data:loaded", function() {
        console.log("finished loaded shapefile");
    });
</script>


texts: Getting Intra-Maps roam to interface with internal GPS

texts:
I've got Intramaps Roam installed on a Panasonic Toughpad fz-g1 with integrated u-blox M8 gps.

When I open up GPSviewer I can see my position, everything seems fine. When I enable GPS with intra-maps roam, it says 'GPS Active' but the lat lon reading is 0,0.

I'm connected on COM3 with 9600bps.

Any ideas what's going wrong?


texts: How can I know the specific information of the monthly average data?

texts:
I know that the "google earth engine" site provides NOAA's VIIRS-DNB VCMSLCFG information.
This is the monthly average data.
The satellite photographs the entire earth twice a day, and I want to know which data was calculated.
In other words, I want to know which images was calculated as monthly average data.
But I can't find out about this.
I ask for help with this.


texts: ArcGIS Silverlight - custom animated composite line symbol?

texts:
I want a 'marching ants' line symbol with a white background and gray foreground. To make custom symbols, see the API Reference, and for examples see the symbol gallery.

I got one in a single color easily enough like this:

    <esri:LineSymbol x:Key="ReferencePointLeaderLineSymbol">
        <esri:LineSymbol.ControlTemplate>
            <!-- see the API docs for the Symbol class. A LineSymbol is expected to have a Path element
            with the name Element at the root of the template. A VisualStateGroup called 'CommonStates'
            can use a state called 'Normal' for normal visual states -->
            <ControlTemplate xmlns:vsm="clr-namespace:System.Windows;assembly=System.Windows">
                <Grid>
                    <vsm:VisualStateManager.VisualStateGroups>
                        <vsm:VisualStateGroup x:Name="CommonStates">
                            <vsm:VisualState x:Name="Normal">
                                <Storyboard RepeatBehavior="Forever">
                                    <DoubleAnimation BeginTime="0:0:0" Storyboard.TargetName="Element" Storyboard.TargetProperty="StrokeDashOffset" To="1000" Duration="0:3:0" />
                                </Storyboard>
                            </vsm:VisualState>
                        </vsm:VisualStateGroup>
                    </vsm:VisualStateManager.VisualStateGroups>
                    <Path x:Name="Element" StrokeDashArray="2,1" StrokeDashOffset="0" Stroke="Gray" StrokeThickness="2" />
                </Grid>
            </ControlTemplate>
        </esri:LineSymbol.ControlTemplate>
    </esri:LineSymbol>

And I mocked up a static version of what I want like:

    <Grid HorizontalAlignment="Left" Margin="8,8,0,0" Width="100" Height="100" VerticalAlignment="Top" Background="ForestGreen">
        <Line StrokeDashArray="1.6,.4" StrokeDashOffset="0.15" Stroke="White" StrokeThickness="3" X1="10" Y1="10" X2="90" Y2="90" />
        <Line StrokeDashArray="2,1" StrokeDashOffset="0" Stroke="Gray"  StrokeThickness="2" X1="10" Y1="10" X2="90" Y2="90" />
    </Grid>

Which looks like this:

But the problem I am facing is that by the above docs, a line symbol needs a SINGLE Path element with the name 'Element' at its root. How can I do a 'white background, grey foreground' symbol without using two Paths? And how would I animate them to march as a unit?


texts: Incorporating Python script into ModelBuilder?

texts:
Within ModelBuilder in ArcGIS, I am trying to add a script that will randomly assign "model types" (which is text) into the new field I have created in the attribute table. I checked to make sure the whole model ran without any problems before adding the script. The script has ran without any errors as well, but from the python window within ArcGIS.

Does anyone know the proper way to add this script to the model?


texts: Comparison of wildfires areas using NBR index obtained by Landsat 8 for 2 dates?

texts:
My aim is to compare two images of NBR index obtained by Landsat 8 for two dates: April the 1st and October the 26th (the beginning of the wildfires period and the ending of the period). After that I want to compare the difference NBR image with MODIS vector products (as it's shown below)
The "blue" pic is a result of ISODATA process (for detecting of the areas with the biggest difference which are supposed to have been burned). The questions are:

Is it a good idea to use ISO DATA for this purpose (I obtained just 3 classes).
How I can compare my vector data (MODIS areas) and raster data (Landsat 8) to find the difference between burned area?


texts: Extra detailed bounding polygon from many geometric points

texts:
I have a database of around 100,000 postcodes each with lat/lon which I have converted to geom type in PostGIS. Each postcode also has one of 5 "regions". I wish to generate an accurate bounding polygon for each region, and have attempted to query it as so:

select region, count(*), ST_AsGeoJSON(ST_ConvexHull(ST_Collect(geom)))
from postcode
group by region;

Despite the counts all being above 10,000 each, the resulting polygons are disappointingly simple and broad (e.g. around 10 sides only), with lots of overlaps.

Is there any way in PostGIS or other where I can achieve this goal of extracting highly detailed polygons? I need to do this because unfortunately the government body in question does not publish open data for these respective regions they have defined.


texts: Opening project shows layers on layer panel but no map

texts: I created a project using QGIS 3.0.3 on Ubuntu 18.04. I can open the project and the layers show up in the Layers panel. However, I don't see any map. I can see the data if I go into Attribute Table but I can't go from there to seeing the map again.


texts: DEM terrain model tool: Band data type Byte not supported in GMT

texts: I am trying to create a hillshade. I am using the DEM terrain model tool, choosing the right input file and write output, but keep getting this error message "ERROR 1: Band data type Byte not supported in GMT, giving up.". What does it mean? Is there a plug-in that I have missed to install or am I writing something wrong?


texts: Calculating percentiles from timeseries raster data using ArcGIS?

texts:
My question is similar with Calculating Percentiles in ArcMap?
But I would like to calculate the percentile from timeseries raster data.

I have timeseries raster data from 1981 to 2016, its the max of daily rainfall in a year, so its just 1 raster data for 1 year.

The idea is to get the rainfall threshold and rainfall-triggered flood with return period, let say for 10 and 25 year. If my data is in excel format, I used =PERCENTILE(DATARANGE,0.9) for calculating the rainfall threshold for 10 year return period of rainfall-triggered flood.

For each cell in x,y location from 1981 to 2016, I would like to calculate the percentile of i.e 90%.

How to calculate the percentile for each cell from timeseries raster using ArcGIS?


texts: GeoExt: Using children node not expanded problem in check function

texts:
I'm using children in a Geoext tree definition.

If the node is expanded the function to see if children is checked or not (see the code) works, however when the node is not expanded it seems that the children are not defined.

Any suggestions as to why this is occurring?

    var childrenTotal = 0, childrenChecked = 0;
    var self = this;
    this.cascade(function (child) {
        // Don't check ourselves nor any other directory nodes that are not cascaders
        if (child == self || (child.hasChildNodes() && !child.attributes.xtype == 'hr_cascader')) {
            return;
        }

        childrenTotal++;
        if (child.ui.isChecked()) {
            childrenChecked++;
        }
    }, null, null);


texts: Drop Linestring geometry from a GeoDataFrame

texts:
I try to create a shapefile of the FAO/GAUL/2015/leve1 country regions of the LMIC's countries.
For that I decided to use Python as it gave me access to the earthengine API and the geopandas librairy.
To create the geodataframe I proceed as such :
import ee
import geemap
import geopandas as gpd

ee.Initialize()

for i, row in df.iterrows():
    # select the country
    ee_country = ee.FeatureCollection("FAO/GAUL/2015/level1").filter(ee.Filter.eq('ADM0_CODE', row.GAUL))
    country_json = geemap.ee_to_geojson(ee_country)
    tmp = gpd.GeoDataFrame.from_features(country_json)

    # add the country to the global gdf
    gdf = tmp if i == 0 else pd.concat([gdf, tmp])

gdf.crs = "EPSG:4326"
gdf

(the df is containing the FAO GAUL ADM0_CODE)
Then I try to transform this GDF into a shapefile using
gdf.to_file('FAO_GAUL_2015_level1.shp')

But I get an error because I have shapely.geometry.collection.GeometryCollection in my gdf (with linestring embed so not shapelyfiable).
Is there a way to only kee the Polygon for revery lin that contain a Geometry collection ?
or Is it possible to drop all the linestring from the geometry column ?


texts: Calculating average value into new field using layer file?

texts:
I have a layer file with a BE field. Some BE have more than one BU. Each BU has a scorecard number.  I need to calculate the average of those scores per BE.  I need this in the layer file to be able to symbolize the data.

I copy a  clip from the table


texts: leaflet change attribution prefix on baselayer change

texts: In leaflet is there a way to change the attribution prefix when a specific baselayer change happens using the default layer control?


texts:
I am using the QGIS Distance Matrix to calculate the distance of points along a line to that specific line. I used the tool before with 2 points layers and it worked well.

This time I created a points-chain layer of the line layer (same coord. system as the points layer, using the QChainage plug-in) and get the following traceback:

Log:

Processing algorithm…
Algorithm 'Distance matrix' starting…
Input parameters:
{ 'INPUT' : 'G:/SMALLWORLD-GTO/working HPa/L8/BH/v3 evaluate/BHwelds_v2.shp|layername=BHwelds_v2', 'INPUT_FIELD' : 'Item_nr', 'MATRIX_TYPE' : 0, 'NEAREST_POINTS' : 1, 'OUTPUT' : 'memory:', 'TARGET' : 'G:/SMALLWORLD-GTO/working HPa/L8/BH/v3 evaluate/20181001_GTO_Pipe Segment Route_dline_CHAIN.shp|layername=20181001_GTO_Pipe Segment Route_dline_CHAIN', 'TARGET_FIELD' : 'cngmeters' }

Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\qgis\PointDistance.py", line 141, in processAlgorithm
matType, nPoints, feedback)
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\qgis\PointDistance.py", line 183, in linearMatrix
index = QgsSpatialIndex(target_source.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(source.sourceCrs(), context.transformContext())), feedback)
Exception: unknown

Execution failed after 6.99 seconds

Loading resulting layers
Algorithm 'Distance matrix' finished

Log end.

Questions are:

1) What does the traceback mean?

2) How can I solve/avoid it?


texts: I remember the browser version coming with the install package in QGIS 2.1.8 (or so??) but cannot find the install for version 3 QGIS anywhere!  Am I missing something?  Please let me know where I can download the browser application..


texts: Highlighting and exporting raster layer based on vector polygon

texts:
I am a Quantum GIS beginner. I have a set of adjacent polygons in a vector layer that is overlaid on a raster layer of OS StreetView OpenData. I want to be able to select each of these polygons and separately export it as an image file (JPEG, PNG, whatever).

In the final image, I do not want any fill colour inside the polygon, but I would like the area immediately surrounding the polygon to appear filled (perhaps with diagonal lines), to clearly mark the boundary.

What is the easiest way I can achieve this with QGIS?

Thanks.

UPDATE: The mapbook plugins highlighted by Darren Cope provide a neat way for me to export the polygons in bulk, which is one major part of the solution. However, the other part of the problem still remains. I need to find a way to shade/fill the area outside each polygon.

I am wondering whether it would be possible to programatically create a polygon layer mask on-the-fly for each polygon, using something like Vector -> Difference.

Any guidance appreciated.


texts: Converting ArcGIS Online Basemap to new coordinate system for ArcMap?

texts:
I have a arcgis we application that depends on all coordinate systems of the layers and basemaps to be the same, which I would think is standard practice.

I like to just use the esri online basemaps such as this one, but it uses 102100 for a coordinate system instead of 4326, which is what I have all my other layers in.

Is there a way I can transform these esri basemaps into a new coordinate system and publish them to my server via ArcMap?

Am I doing something wrong here and I should be taking a whole different approach, because I would expect other people to run into the same issue as this all the time if they want to use ESRI basemaps.


texts: Merging shapefiles as Virtual Layer in QGIS

texts:
In QGIS 3.6.3 I'm trying to create a Virtual Layer obtained by merging 12 shapefiles (with the columns named and ordered in the same way) using the expression

SELECT *
FROM _January, _February, _March, _April, _May, _June, _July, _August, _September,
_November, _October, _December

but the software crashes.

If I wanted to create a Virtual Layer by taking data from a single column called "TEST", how can I modify the expression? I could not find the correct form.


texts:
I want to download a tiff image from this WMS server with the GetMap-Method:

http://services.lgrb-bw.de/index.phtml?SERVICE_NAME=lgrb_corg&REQUEST=GetCapabilities&VERSION=1.1.1&SERVICE=WMS

After searching for answer, especially here on Stackexchange I managed to write this request and could download the image: http://www4.lgrb.uni-freiburg.de/serverbase/services/index.phtml?SERVICE_NAME=lgrb_corg&amp&request=GetMap&service=WMS&version=1.1.1&layers=cg004&format=image%2Ftiff&STYLE=default&SRS=EPSG:4326&bbox=6.83803,47.2802,8.0206,48.9963&width=2048&height=2048.

But unfortunately the map is completely white. Any suggestions what could be the reason?


texts: Using If Then Else statement against selected records in ArcPy?

texts:
I am working on a Python script that will count how many selected features there are and if there are it will make a feature layer from them.  If the script runs and finds there are 0 selected records it should print 'NO SELECTED RECORDS'.  However when I run my script it is ignoring my If Then statement and creating a duplicate layer of "Locates", but still telling me I have no records selected.

Here is my code:

import arcpy

totalCount = arcpy.GetCount_management("Locates")
print 'Total Record Count = ' + `totalCount`
try:
    arcpy.SelectLayerByLocation_management("Locates","INTERSECT","Boxes",0,"NEW_SELECTION","INVERT")
    arcpy.SelectLayerByLocation_management("Locates","INTERSECT","Lines",0,"REMOVE_FROM_SELECTION")
    ##GET COUNT OF SELECTED RECORDS
    selectedCount = len(arcpy.Describe("Locates").FIDSet.split('; '))-1
    print 'Selected Records Count' + `selectedCount`
    if selectedCount != totalCount:
        ##IF SELECTED ROWS NOT EQUAL TO 0 THEN MAKE FEATURE LAYER ELSE PRINT 'NO LOCATES TO CLOSE'
        arcpy.MakeFeatureLayer_management("Locates","NeedToCLOSE")
        for row in arcpy.SearchCursor("NeedToCLOSE"):
            print 'Close Work Order #' + `row.WorkOrderID`
    else:
        print 'No Locates to close'
except:
    print(arcpy.GetMessages())

What do I need to change in my code for this to work?


texts: How do I change the display unit without changing the coordinates numbers?

texts:
(X= 716 562,800911  Y= 3 568 493,622591 "Decimal Degrees")
My problem is that the coordinates are correct but displayed with the wrong unit, it should be meters. When I try to use the classic way to change it (View-> Data Frame Properties... etc.) it changes the coordinates also.
Is there any way I can change only the display unit?


texts: How to get 25 points within 1 square kilometer in Arcmap

texts: I tried to generate 25 sampling points within  1sqkm using distance point to point is 250m but i got only 15 points not 25 points as i expected because starting point was not at the corner. The tools I used was  Creat fishnet.


texts: ST_ConcaveHull produces "No function matches given name" in PostGIS?

texts:
I'm using Postgis 2.2.2. This is my table:

county       | character varying       |
wkb_geometry | geometry(Polygon,27700) |

I want to get the area of the concave hull around all polygons in a particular county.

This is what I'm trying:

select st_area(st_concavehull(st_multi(st_collect(wkb_geometry))))
from mytable where county='XX';

But I'm getting:

HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

I think the problem is with st_concavehull. I've also tried just st_concavehull(wkb_geometry), but that produces the same error message.

What should I be doing instead?


texts: PostgreSQL PostGIS Query is taking long time

texts:
I have a table that included only 472 rows and geometry of table is Multipolygon. I am sending query to database via PgAdmin III

select * from mytable // this query takes 99413 ms

This is a very long time.

What do you think about long queries on PostgreSQL? Server resource, disc, CPU etc.


texts: Calling interpolation plugin from Python console of QGIS

texts:
I would like to call the QGIS interpolation plugin function (TIN method)  (Raster->Interpolate)  from  python console.

I can't find the corresponding function within QGIS API or within the processing algorithm list.  I found the SAGA Triangulation algorithm, which works fine but is 5-10 x slower and speed is important in my case.

Any idea how to execute it?


texts: Creating summary table of symbology of all datasets in multiple MXDs?

texts:
This is related to the following Retain ESRI symbology in Mapinfo TAB . As an initial step (in arcpy) I want to access the symbology of all the datasets and create a document which has the details of the symbology as per the graphic below.

Is this possible (should I invest the time to work it out) or if it isn't possible I will just have to get cracking on doing it manually on hundreds of MXD's (:-

To break down the steps -hopefully I can automate at least part of the process.

scrolling through mxd's and getting the datasets that are on is no problem

lyrFile = arcpy.mapping.Layer
if lyrFile.visible == True:
as it finds a dataset it should create a .lyr file (as most are just shp/gdb files) and then AttribXXX=lyrFile.xxx

where XXX and xxx are

BRIGHTNESS = lyrFile.brightness

(above format to following)

CONTRAST —contrast
DATASETNAME —datasetName
DATASOURCE —dataSource
DESCRIPTION —description
LABELCLASSES —labelClasses
LONGNAME —longName
NAME —name
SHOWLABELS —showLabels
TRANSPARENCY —transparency
WORKSPACEPATH —A layer's workspace or connection file path

The issue is I can't find anything that gets the RGB, line/point/polygon type, line weights, font etc.  If these can be found then the next step

write each of the values above to a CSV (easy)...for example (this has to be expanded)

                f = open(origPath+'\\LayerAtrributes_LOG.txt', 'a')
                f.write(str(DATASETNAME)+","+str(DATASOURCE)+","+str(DESCRIPTION)+","+str(MXDName)+"\n")
                f.close()

If possible a jpg of the symbology should be produced as in the example and the link to it added to the csv.

Am is asking to much of arcpy?


texts: Toggling layers with legends in CartoDB.js?

texts:
I created a map on the CartoDB website with multiple layers, each with its own legend. I now wish to use CartoDB.js to be able to toggle between the layers, so the legend will toggle on with the corresponding layer.

I used cartodb.createVis() to load the map and its layers. I'm able to show and hide layers using layer.show() and layer.hide(), but the legend(s) shown are those for the visible layer(s) on the CartoDB website version of the map.

Could you please tell me how to toggle the legend to correspond with the visible CartoDB.js layer.


texts: How to calculate area of a polygon in a shapefile in QGIS

texts: I want to calculate areas of different polygons in QGIS in tabular format.


texts: Arguments error when running LineDensity scripts in arcpy

texts:
I'm trying to use LineDensity tool on polyline features. The shapefile is small and includes some invalid geometries.

The code I use is an adaptation of an ArcGIS Help script.
Here is my code:

    import arcpy
    from arcpy import env
    from arcpy.sa import *

    env.workspace = r"Z:\...\...\fdc\Grand_Ouest\33"

    # Parameters
    inPolylineFeatures = "_polyline.shp"
    populationField = "NONE"
    cellSize = 200
    searchRadius = 500
    areaUnitScaleFactor = "SQUARE_KILOMETERS"

    # Check out the ArcGIS Spatial Analyst extension license
    arcpy.CheckOutExtension("Spatial")

    # Execute LineDensity
    outLineDensity = LineDensity(inPolylineFeatures,
                                 populationField,
                                 cellSize ,
                                 searchRadius ,
                                 areaUnitScaleFactor)

    # Save the output
    outLineDensity.save(r"Z:\...\...\test.gdb\linedensity")

Here is the error I get from the code:

  Traceback (most recent call last):
   File "Z:/.../.../python/in_progress/linedensity_test.py", line 22, in
     areaUnitScaleFactor)
   File "C:\Program Files (x86)\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 446, in LineDensity
     area_unit_scale_factor)
   File "C:\Program Files (x86)\Desktop10.3\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
     result = wrapper(*args, **kwargs)
    File "C:\Program Files (x86)\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 439, in Wrapper
      area_unit_scale_factor)
   File "C:\Program Files (x86)\Desktop10.3\ArcPy\arcpy\geoprocessing_base.py", line 504, in
     return lambda *args: val(*gp_fixargs(args, True))
  ExecuteError: ERROR 999999: Error executing function.
  Failed to open raster dataset
  Failed to execute (LineDensity).

I tried and modified my code more than once, primarily regarding the call of parameters and the use of environment module.

My last move was to discretize the arguments of LineDensity function by organizing them in separate lines. Now here we are: the only argument reported in the error is the last one, the area unit scale factor.

But that's the beginning of another question: cause the projection (RGF_1993_Lambert_93, epsg:2154), its measure unit (meter) and the scale factor (SQUARE_KILOMETERS) are perfectly consistent.

I suppose I've to try and find a new approach, with different data.

Little is to be found online on that matter except the link I delivered you at the beginning of the question.


texts:
How can I prepare data to generate a DTM using points and breaklines, but where every breakline has nodes at different heights. For example, a river valley, which regularly falls,  or a dorsal line, where each vertex is a different height.

Thanks for the answer. Tom


texts: How to batch generate thumbnails of raster layers with R?

texts:
I have been using R to maintain a rather large stack of raster layers (over 700 layers of data for sub-Saharan Africa). Each layer looks as follows http://hcapi.harvestchoice.org/v2.0/core/?varCode=SALUM&format=tif

I need to batch generate PNG images for print media, as well as thumbnails of the layers.

Using R raster package, I can easily plot each layer and apply custom color tables and legends. I run into problems when trying to generate small thumbnails with no title and axis (see an example below). I cannot find any option to have the raster plot occupy the full size of the graphic device, regardless of any graphic parameter I pass to par().

Here is reproducible code:

library(rgdal)
library(raster)

url <- "http://hcapi.harvestchoice.org/v2.0/core/?varCode=SALUM&format=tif"
download.file(url, "./r.zip", "curl")
unzip("r.zip")
r <- raster("./SALUMSSAtif.tif")
cc <- c("#AFF1E9", "#9ddfa3", "#87CD58", "#CDAF19", "#761200", "#94725F", "#D61F00")

png("./SALUMSSAtif.png", width=200, height=200)
par(mar=c(0,0,0,0), oma=c(0,0,0,0), xaxs="i", yaxs="i")
plot(r, col=cc, main=NULL, legend=F, axes=F)
dev.off()

This produces the attached plot that's badly off-center and too small. Any idea how to remove all the plot margins using raster::plot or another method?

Thanks, --Mel.


texts: Creating a custom "flat" projection definition in proj4

texts:
I want to create a simple, custom projection in proj4. I need it to be completly flat ( rectangular coordinates ). Starts in point (0,0), meters as unit, everything else - as simple as possible. What parameters my definition should have to achieve it?

Am I even able to do it? Ignoring the curvature of the Earth?


texts: "Couldn't load PyQGIS" Error running QGIS on Ubuntu 14.04

texts:
I have installed QGis on Ubuntu 14.04 from the QGis stable repository; it went fine with no errors reported. But when I launch QGis I get this huge window with the following message:

  Couldn't load PyQGIS. Python support will be disabled.

  Traceback (most recent call last):   File "", line 1, in    File
  "/usr/lib/python2.7/dist-packages/qgis/init.py", line 35, in
      from qgis.core import QgsFeature, QgsGeometry ImportError: /usr/lib/python2.7/dist-packages/qgis/core.so: undefined symbol:
  _ZN16QgsComposerShape12setSceneRectERK6QRectF

  Python version:
  2.7.6 (default, Mar 22 2014, 23:03:41)  [GCC 4.8.2]

  QGIS version:
  2.2.0-Valmiera 'Valmiera', exported

  Python path: ['/usr/share/qgis/python',
  u'/home/desousa/.qgis2/python',
  u'/home/desousa/.qgis2/python/plugins',
  '/usr/share/qgis/python/plugins',
  '/usr/local/lib/python2.7/dist-packages/MapScript-6.0.1-py2.7-linux-x86_64.egg',
  '/usr/local/lib/python2.7/dist-packages/virtualenv-1.8.4-py2.7.egg',
  '/usr/local/lib/python2.7/dist-packages/pip-1.2.1-py2.7.egg',
  '/usr/local/lib/python2.7/dist-packages/pytz-2012j-py2.7.egg',
  '/usr/local/lib/python2.7/dist-packages/python_dateutil-2.1-py2.7.egg',
  '/usr/local/lib/python2.7/dist-packages/six-1.3.0-py2.7.egg',
  '/usr/lib/grass64/etc/python/grass/script', '/usr/lib/grass64/etc',
  '/usr/lib/grass64/etc/python', '/usr/lib/grass64/lib',
  '/usr/lib/grass64/bin', '/usr/lib/grass64/scripts',
  '/usr/lib/python2.6', '/home/desousa/Eclipse/Workspace.Tudor/OWSLib',
  '/usr/lib/python2.7', '/usr/lib/python2.7/plat-x86_64-linux-gnu',
  '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old',
  '/usr/lib/python2.7/lib-dynload',
  '/usr/local/lib/python2.7/dist-packages',
  '/usr/lib/python2.7/dist-packages',
  '/usr/lib/python2.7/dist-packages/PILcompat',
  '/usr/lib/python2.7/dist-packages/gst-0.10',
  '/usr/lib/python2.7/dist-packages/gtk-2.0',
  '/usr/lib/pymodules/python2.7',
  '/usr/lib/python2.7/dist-packages/ubuntu-sso-client',
  '/usr/lib/python2.7/dist-packages/wx-2.8-gtk2-unicode']

The packages I have installed:

$ dpkg -l | grep qgis
ii  libqgis2.2.0                                          2.2.0-1~precise1                                    amd64        QGIS - shared libraries
ii  python-qgis                                           2.2.0-3+trusty1                                     amd64        Python bindings to QGIS
ii  python-qgis-common                                    2.2.0-3+trusty1                                     all          Python bindings to QGIS - architecture-independent files
ii  qgis                                                  2.2.0-3+trusty1                                     amd64        Geographic Information System (GIS)
ii  qgis-common                                           2.2.0-3+trusty1                                     all          QGIS - architecture-independent data
ii  qgis-plugin-globe                                     2.2.0-3+trusty1                                     amd64        OSG globe plugin for QGIS
ii  qgis-plugin-globe-common                              2.2.0-3+trusty1                                     all          OSG GLOBE plugin for QGIS - architecture-independent data
ii  qgis-plugin-grass                                     2.2.0-3+trusty1                                     amd64        GRASS plugin for QGIS
ii  qgis-plugin-grass-common                              2.2.0-3+trusty1                                     all          GRASS plugin for QGIS - architecture-independent data
ii  qgis-providers                                        2.2.0-3+trusty1                                     amd64        collection of data providers to QGIS
ii  qgis-providers-common                                 2.2.0-3+trusty1                                     all          collection of data providers to QGIS - architecture-independent files


texts: Using Google Earth KMZ file in Power BI, vector and Earth images combined

texts:
I have a Power BI project that contains a KML file from Google  Earth (kmz file). I just need the Power BI to show the shape in the same way Google Earth shows it. When I try to put it in Power BI I lost the satellite feature presented by the Google Earth software, I just get the vector/polygon map with no satellite image.

Is it possible to present the KML file in Power BI in the same way Google Earth presents it?

I uploaded a video to try to explain it better...

link with the question:
https://youtu.be/hp1zry38EZM


texts: Is there a dissolve / ST_Union function that will close gaps between features?

texts:
I've been reading a number of related threads about ST_Union, ST_SnapToGrid and ST_Dump but I haven't been able to solve my problem.

The problem is I am trying to merge Tiger/line census tracts that have been simplified but I get gaps where the tracts aren't flush. I simplified the tracts with ST_SimplifyPreserveTopology(geom, .0003) and I'm running ST_Union(ST_SnapToGrid(geom, .00001)) but I still get gaps. I've also tried using ST_SnapToGrid around of the initial simplification but I don't see a difference. I've also tried wrapping ST_Dump around the Union but I don't see a difference there either.

I thought ST_ExteriorRing was promising but my GROUP BY creates three concentric polygons (doughnuts) and ST_ExteriorRing seems to create three overlapping circles. Unless I'm doing something wrong.

I'm outputting my query as geojson and putting that on a leaflet map. Here is what those gaps look like at different zoom levels.

If anyone has any ideas on making a gapless polygon of census tracts under 80mb I'm all ears.


texts:
I have a WMS query (GeoServer) that returns info from my layer into a popup in Leaflet.JS.However I was wondering how would do this outside of leaflet, i.e. on a standard html/php page.

I assume I would use the ST_buffer function in POST GIS when I do my query, order by distance and then limit by 1. For example I have a table with the_geom and lets say house name, How would I do a simple buffer that selects the nearest house to that given point which happens to be from the same table using the same geom column. I think if I can just get that working I should be able to extend the query myself to provide the info I need.

Thanks


texts: Exporting ArcMap file to PDF with all layers

texts:
Below I posted pictures of my arcmap layers and PDF export.
My map is measuring sea level rise so all my layers are rasters.
How would I export my ArcMap with my layers so I can turn it on and off in PDF?
I'm using 10.6.1 for ArcMap and my Adobe Reader is 9.3 and seems to be updated. What can I do to incorporate all my layers?
I'm a beginner to GIS and I need my project to be exported.
As for as exporting, I am aware of the general and advanced tab while saving it. Even with all those options clicked, nothing is working. The only thing that gets exported are the vector layer and the sea level rise (10m) layer gets exported as an image.
If I have to export it one layer at a time as an image, how do I import all the images in one PDF file?


texts: remove whitespace Python code results in odd character - Field Calculator

texts:
I need to remove whitespace from the end of a string in a field in ArcMap.

I tried first the VBA function RTrim, but whitespace still remained.  I thought maybe the whitespaces were tabs, so I tried 'Replace([field], vbTab, ""), but no luck.  Then I tried the Python function rstrip, as so:

!FIELDNAME!.rstrip()

..but this resulted in the whitespace being replaced by this character:

Â

There are several of these characters in each of the fields cells trailing the string.
I have looked online for some meaning of this symbol, but I haven't found any.  Doe anyone know what this means? I have still not successfully removed the trailing whitespace.

here is a screenshot of the attribute column after I run the Python code given above:


texts:
I have what appears to be a simple question, but which I can't seem to determine an answer myself.  I have attributes in the polygons which represent multiple points over time (say January has a value of 1000, February has a value of 900, etc...)

I would like to know if there is an easy way to display/toggle only January time points, or February time points, etc.  I could theoretically create a layer for every month, but this is an impractical solution because the file would become too large.

The ultimate goal would be to export the map using qgis2web, and let the user toggle the month he/she wishes to see, with the corresponding attribute displayed.


texts: 10 GB Raster Aerial Image in Quantum and POSTGIS

texts:
I have a problem... I have been successful in uploading the raster in postgis but when I access it through Qgis. it does not open even after several hours.

I have 10GB aerial image raster and everytime i try to open it through QGIS it does not open even after many hours and QGIS crashes...

why is this happening...?? is it really due to the file size..??

Regards,
Omer


texts: Performing a function over all pair points of a grid, R

texts: I need to perform a mathematical formula over each pair points of a grid of 100*100 points, for example calculate distances of all two points and having them in a list. How I can do this in R?


texts:
I'm customizing the cmv viewer to works with my services.
I have service that have 3 layer (0,1,2) and is a dynamic map service layer.
It's possible to set only one layer to show?


texts: Reason of mismatch of projection system of rasterlayer and featurelayer?

texts:
What I would like to check is if the raster layer and feature layer have the same projection.
Code developed so far:
public bool IsSpatialReferenceSame(ISpatialReference sp1, ISpatialReference sp2)
{
  return sp1.Equals(sp2);
}

In SpatialReference Class:
public bool Equals(ISpatialReference sp2)
{
  var sp1Clone = this.Clone;
  var sp2Clone = sp2.Clone;
  return sp1Clone.IsEqual(sp2Clone);
}

When I checked the feature layer and raster layer spatial reference in layer properties, I saw this.
Feature Layer Properties:

Raster Layer Properties:

But when I call IsSpatialReferenceSame it returns false.


texts: ST_MakeLine based on a field or equivalent

texts:
I want to create a set of lines taking into account the value of an attribute. As soon as the value of this attribute changes, then I start a new line. Currently, I use the FME LineBuilder tool where I specify a field in the "Connection Break Attributes" settings.

I would like to do the same thing using Python or PostGIS.
I have explored several courses on the PostGIS side. I saw that I could use the ST_MakeLine function by using the value of a field. So, I tried to achieve the following result first. As long as the value of the "NumAntennas" field does not change then the "Shape" field takes the same value, otherwise I have an increment of 1, etc.

I tried to do this with functions that would look like row_number() or dense_rank() but I didn't succeed. I don't know if there is a function that would meet my need or if I need to make a loop in PL/pgSQL or Python.


texts: Making a grid with points that are staggered every 800 feet

texts: I need to make a grid pattern with points that stagger to make a diamond pattern instead of the regular pattern that a fishnet grid makes. It is sometimes referred to a "staggered start". My points need to be 800ft away from each other. I am creating my map in ArcGIS 10.2. Is there something else I can try?


texts: Extracting GPX / KML routes from embedded Google Map

texts:
I want to extract the route data that is shown below to a GPX or KML file. Being that this map is programmatically created and not a "Google My Maps" type of embedded map, is there a way to download these routes for offline use?
http://gts.groomertrackingsystems.com/map


texts:
I am having the following issue writing points to a SQLite file using ogr with SPATIALITE enabed:

  ERROR 1: Cannot insert feature with geometry of type POINTZ in column GEOMETRY. Type POINT expected

Below is my code, this code was previously used to create a GeoPackage with exactly the same fields. (But not with Spatialite enabled because it doesn't have it)

I have commented the code to indicate where the point is created and where the error message comes up

import ogr
import numpy as np
import os
DataSource= ogr.GetDriverByName("SQLite").CreateDataSource(os.path.join(OutputFolder,"DataFile.sqlite"), options = ["SPATIALITE=YES"])
DataLayer = DataSource.CreateLayer("Data_Points",proj,geom_type=ogr.wkbPoint)

DataLayer.CreateField(ogr.FieldDefn('ID',ogr.OFTInteger64))
DataLayer.CreateField(ogr.FieldDefn('buffer_dist',ogr.OFTInteger64))
DataLayer.CreateField(ogr.FieldDefn('arr_id',ogr.OFTInteger64))

Data_layer_defn = DataLayer.GetLayerDefn()
DataSource.StartTransaction()

data = np.load(os.path.join(OutputFolder,"Output0.npy"))
arr_id = 0
for row in range(len(data)):

    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(float(data[row,1]),float(data[row,2])) #x,y Setting the point

    Data_Feature = ogr.Feature(Data_layer_defn)
    Data_Feature.SetField("ID", int(data[row,0]))
    Data_Feature.SetField("buffer_dist", int(data[row,3]))
    Data_Feature.SetField("arr_id", arr_id)
    Data_Feature.SetGeometry(point)
    DataLayer.CreateFeature(Data_Feature) #Error message happens on this line

    point = None
    arr_id +=1

DataSource.CommitTransaction()
DataSource.FlushCache()


texts: Inaccurate raster values from "SAGA Extract Raster Values to Points"

texts:
I am perplexed by the results generated by the SAGA 2.3.1 Extract Values to Points tool in QGIS 2.18.

I am working with Vegetation Index raster layers, which means the values are between -1 and 1.  According the min-max of the raster values (prior to the SAGA calculation), all the minimum values are above 0. The Extract Values to Points tool generates negative values up to -3.  I have tried saving the No Data values to 0, but that doesn't change the output.

Has anyone else run into this problem? I understand that there was a known bug in the algorithm; however, I don't believe this is my issue.  The tool is generally working, except for the occasional new negative numbers (I say new because they were not present in the original input raster). Additionally, I know people also recommend using the Point Sampling plugin in QGIS; however, I am working with huge datasets and (according to this forum), the plugin cannot handle that type of volume.

This issue has me questioning my results.


texts: QGIS core plugins vs. normal plugin: what is the difference?

texts: In QGIS, besides the regular plugins, there are so-called core-plugins that can't be de-installed. So from a user point of view it's clear what is the difference: you don't have to install (just activate) them. However, less clear to me is what from a technical point of view is the difference: what is the advantage of implementing some of the plugins as core plugins?


texts: Filtering elements of FeatureCollection using fusion table that contains list?

texts:
I'm new to Google Earth Engine.

I am trying to filter all countries in the Americas from a Fusion Table that contains all countries in the World:

    // Load countries
    var countries = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw')

This fusion table has a column named ISO3 with the code for each country but doesn't have a variable indicating region.

In addition, in GitHub I found a table that contains the ISO3 codes together with a variable for each region (Americas=19 in column "region-code"):

ISO-3166-Countries-with-Regional-Codes/all/all.csv

I don't know how to get this table directly into Google Earth Engine, but I can manually add it as a fusion table name ISO3_Americas containing 1 column "ISO3" with the list of all countries in the Americas:

https://www.google.com/fusiontables/DataSource?docid=1_BoMJBCnMW72uy85jczMBHwH7lFe64M5_b_n6BSy

I would like to take the FeatureCollection "countries" and apply a filter: "keep if ISO3=any of the countries listed in ISO3_Americas"


texts: ArcGIS 10 MapControls in MS Access, can it be done?

texts:
In the past (9.3 ArcGIS) I was able to develop forms in MS Access where I could drop a MapControl on and have GIS functionality all within Access as long as ArcMap was installed on the same machine.

I now have ArcGIS 10 and Office 2003 and when I start Access go to  VB editor and create a form I am unable to add a MapContol, it says I do not have the appropriate license.  I've had a search on the ESRI website and Help and I have not found a clear answer.

Am I correct in understanding that they have changed something behind the scenes with regards to the MapControl, one cannot simply use it anymore? I've tried dropping an ESRI LicenseControl on the form and I get the same message rejecting me as I have no license to use the LicenseControl!

It would be good if someone clarifies this and explain how one can embed a MapControl into a MS Access form?


texts: Methods for burning depressions (not rivers) in a DEM

texts:
I am using a 2x2m DEM to simulate runoff on a 2D hydrodynamic model. After a baseline simulation, I am interested in evaluating the impact of adding depressions in specific areas that may work as small catching basins that accumulate surface runoff. See for example the photo below showing a flooded street near a park, and its DEM:

What I would like to do, is to carve a depression so that the park in the intersection acts as a catchment of the water that has accumulated on the roads. A quick and dirty approach would be to digitize the park's boundaries, convert it to a raster with a flat, constant value, and then use that value to substract it from the original DEM. However, this would mean that the park's topography gets lowered, but remains in shape (e.g. the raised tree pits on the Northeast that block water would be lower, but still exist).
In addition, just lowering the park's limits could mean that water does not accumulate in the new depression if the slope on the roads does not lead to it, so I infer that I should make small alterations on the roads' DEM area in order to redirect flow towards the central part of the park.
Acknowledging that DEM manipulation can be tricky, I am unsure how to approach this. regarding issue #1 (flattening the park before lowering it to create the depression) I think I could just create a raster and calculate the difference between each cell's value and the area's minimum. This way i would remove the park's topography.
Regarding issue #2 (how to make small, suddle alterations on the surrounding area in order to direct flow towards the new basin) I am completely clueless on how to operate in an efficient manner. Is there a method / tool that could be used to, for example, convert a step into a contiguous slope? (see diagrams)
I am working coordinately on ArcGIS and R, so solutions that use any of these options / combine them are welcome.


texts:
On a shared computer, I have a QGIS project which links up to a PostGIS database that contains my data. If a user wants to access layers from the db, they have to enter the db credentials, so that's OK.

However, I have previously logged on and imported layers from my db into QGIS, and manipulated them in QGIS (styling options based on meta-data, etc.)... Anybody who opens my QGIS project can view these layers in QGIS, and can therefore see the data I wish to be protected in my db.

Is there a way to password protect my QGIS project?


texts: Splitting polygon into multiple fixed-width polygons along edge in R?

texts:
I have a polygon and would like to spit into grid with fixed width along with one edge. The result will look like the figure below.

Are there algorithm to achieve my target? My language is R, but welcome to any other languages.


texts: Performing statistics or simple pie chart in QGIS?

texts: Does anyone know how to perform pie charts based on table values within QGIS and include them in the print outs?


texts:
I am using QGIS 2.18.9 and I have a problem. I dissolved a small map with about 1000 small regions into one big polygon, with Vector->geoprocessing tools->dissolve.

But after the dissolving there are little dots everywhere. After zooming in I found that these are all little holes in the layer.

Now I want the dissolved layer to not have any holes. Is there an easy way to do that?


texts: Trying to iterate through the nodes of a multipolygon

texts:
I've got a multipolygon which I want to iterate through the nodes for. I will be editing the nodes later on, but right now I just want some code which works.

My relevant code is just this bit:

for featureW in layerWorking.getFeatures():
if featureW['SW_MEMBER'] == My_SW_MEMBER:
    print('Found <<TARGET>> feature')

So I've found the feature I want to target, and want to iterate through the nodes in a loop to, say, get their coordinates.

It should be noted that, whilst the layer TYPE is multipolygon (as there are a handful there), the shapes I'll be looking at are just individual polygons. I'm also using QGIS 3.4.5


texts: Need Dominion Land Survey (DLS) grids loaded into QGIS for province of Saskatchewan in Canada

texts:
I found a source that allows the download of a shapefile for the grid system for Alberta, but am still in need of the same thing for the entire Dominion Land Survey Grid, or just Saskatchewan.

I'm working in the Agriculture industry on boundaries for fields. There are resources for locating specific locations, but I need to be able to create the boundaries in QGIS.


texts: Incorrect spatial extent for imported HDF image on ArcGIS Desktop?

texts:
I imported an hdf image onto ArcGIS but I cannot make correct projection of it. Also, it does not have correct spatial extent (please see the enclosed figure, the small one is projected under WGS84). I try to define the projection of the hdf raster but the tool just simply doesn't run. The hdf raster has 7200 columns and 3600 rows while the small one has 3600 columns and 1800 rows.

Do you have any ideas how I can fix this?


texts: When was Esri File Geodatabase API made available?

texts:
When was the Esri File Geodatabase API made available?

This question was originally asked to try and determine when the File Geodatabase API would be made available.


texts: Defining area not to show labels in using QGIS?

texts:
In a QGIS 3 layout is it possible to define an area (that's being used as an inset)  where labels should not be shown?
So I have

I'd assume that Position -- Above Map labels would work but the labels still show. The stacking layer is activated if we choose Above Map Layer but then the overview box disappears if I choose the labelled layer.
I know I can make it a solid block but I want the data to show underneath but not the labels.


texts: gdal2tiles is not recognized as an internal or external command?

texts:
I try to follow this guide ► Creating super-overlays with gdal2tiles, everything goes well until it reachs gdal2tiles part.
When I type the command gdal2tiles.py -p geodetic -k bluemarble2.vrt or gdal2tiles.py or gdal2tiles in OSGeo4W Shell it gives following result:

'gdal2tiles' is not recognized as an internal or external command, operable program or batch file.

I tried almost everything, searched Google for 2 days, added paths in Windows Environment Variables, installed-uninstalled separate OSGeo4W, experimented with different python versions. And still no any improvement.
Details: At the moment my OSGeoW4 (path: "C:\Program Files\QGIS 3.12\OSGeo4W.bat") is the one which comes with QGIS 3.12, and this installation has its own Python 3.7 (path: "C:\Program Files\QGIS 3.12\apps\Python37\python.exe") as well.
As google searching showed I am not the only one who has such problem, so does anyone know how to solve this issue?
► UPDATE [08:20 2020-06-18]
After @ian-turton's answer I've installed python bindings, but now when I run the command gdal2tiles.py -p geodetic -k bluemarble2.vrt it opens actual python file gdal2tiles.py in my VS Code editor.


texts: Routing Overlapping Transit Lines in ArcGIS Desktop?

texts:
I am trying to model a local bus network in ArcGIS 10.2 and using the Network Analyst Extension to find the shortest path. There are sections of a line in which a bus passes and then returns to form a kind of loop on the same road section before continuing on the journey. I tried to model this by copying the edges twice and giving a 0 Z-elevation on the way in and 1 Z-elevation on the way out.

Although this method has worked when i put the stops outside this "loop"  sometimes it is not working correctly (shortest time) when the stop is inside the loop, rather than stopping on the first pass, it goes all the way inside and then comes back and stops on the way out. The impedance value is set to Minutes.

Any ideas why this may be so?


texts: Group geometries of line data based on one common attribute and compute length

texts:
I have a vector file of river networks in QGIS. Some of the rivers are split up it multiple rows. What I'd like to do is to aggregate the rivers based on the attribute in one column (S_NAME) and afterwards calculate the length of each aggregated river.
I know there must be some PostGIS commands to do that. I'd be fine with any pointer to anything that relates to it. The pseudo-code I have in my head would look something like this:
output = layer -> groupyby(S_NAME) -> st_distance(.)

I know this is really bad, but maybe someone has any idea:.
The link to the data (18k rows...) is here:
https://geoportal.bafg.de/inspire/download/reporting_units/riverWaterBodyForWFD/datasetfeed.xml


texts: Problem with adding non spatial view from MSSQL to QGIS

texts:
I have QGIS 2.12.3 app installed on one server (servername1) and MSSQL server database installed on servername2. The database is password protected and I have made connection successfully using QGIS function to Add MSSQL spatial layer.  I want to add non-spatial views to QGIS and then make joins with spatial data layers. Problem is that some views I can add, and some I can not. The message I got from QGIS is:sql= is an invalid layer and cannot be loaded. I tried to use phyton script to add layers, but it doesn't help:

  uri = MSSQL:server=servername2;database='LCSPI_test_2016';tables=vw_LCZAJ_KIS_GIS_Anagrafi2;trusted_connection=yes"
  qgis.utils.iface.addVectorLayer(uri,'vw_LCZAJ_KIS_GIS_Anagrafi2','ogr')

I have also tried to insert in phyton script username and password of the database but nothing.

When I make all this on my computer (restore the same database, instal QGIS) I can add problematic views using Phyton script and trusted connection (so layers are ok, not invalid). Does anyone know what can be a problem?


texts: How to add OpenSeaMap tiles as an overlay on top of osm base map with Leaflet?

texts:
After generating the OpenSeaMap tiles, I tried to add them as an overlay upon the openstreet basemap. The idea is after generating the OpenSeaMap tiles to serve them as "static" tiles to the client, from an image folder under /var/www/html from apache.

I copied the tiles/ folder under /var/www/html/map/ (in the same location with index.html) and added

in the section script of index.html the following:

var map = L.map('map');
map.setView([53.45 , 14.55],5);

var basemap = L.tileLayer('http://localhost/osm_tiles/{z}/{x}/{y}.png', {
maxZoom: 18
)}.addTo(map);

var sea = L.tileLayer('tiles/{z}/{x}/{y}.png', {
minZoom: 9,
maxZoom: 18
)}.addTo(map);

but I see in browser only the openstreetmap layer not the seamarks layer. But no error (pink page). I may also open openseamap tiles and view as pictures. Is something wrong with the file structure (see below)? What is wrong?

I also noticed that the x , y (coordinates) of the OpenSeaMap tiles are large numbers like this:

under 18/ (zoom)
I have  219132/ (x)
and then 120318.png

Instead the structure of osm tiles is a little bit different:

18/ (zoom)

under 18/ there exist 33/72/75/80/128.meta


texts: OpenLayers not displaying map from WMS layer in GeoServer

texts:
I've adapted the working example from geo-solutions website but still it is not working, it seems the problem is with jetty. below is the previous question

I'm using Geoserver with jetty web server, I've published WMS layer in GeoServer and when i call it from QGIS it works perfectly, but when i call it with OpenLayers it is not loading, it says HTTP ERROR: 404 Not found, both my html file (with the OpenLayers script) and GeoServer are running on localhost:8083.

         <title>Counties of enschede</title>
         <link rel="stylesheet" href="localhost:8083/www/openlayers/theme/default/style.css" type="text/css">

Counties of enschede


texts:
What kind of visualizations/statistics can I get from PostGIS' GiST index?

If we look at Wikipedia's R-tree entry, there are different R-tree variations and visualizations of bounding boxes of each using ELKI. GiST index is similar to R-Tree index. So is there some method/ extension for this kind of visualization (perhaps through QGIS)?

Edit: Here it is also visualized - how is this achieved?


texts: Copying attributes from point layer to memory linestring layer PYQGIS

texts:
I've stacked a little bit with copying attributes.

I have a Point Layer created from .csv file. For each feature in this layer, I'm creating a new Polyline geometry in memory layer. Cannot add attributes from Point layer to newly created Polyline.

Here's my code:

#mcvLayer - Point Layer

nl = Create_layer("New_line", "LineString")
nl.show_layer
hh = iface.activeLayer()
a = QgsPoint(0,0) # Here will be another point

for point in mcvLayer.dataProvider().getFeatures():
        pattrs = point.attributes()
        pgeom = point.geometry().asPoint()
        nl.create_line(a, pgeom, pattrs)

class Create_layer(object):

def __init__(self,name,type):
    self.type = type
    self.name = name
    self.layer = QgsVectorLayer(self.type, self.name, "memory")
    self.pr = self.layer.dataProvider()
def create_line(self, startpoint, endpoint, attributelist):
    self.feat = QgsFeature()
    self.feat.setGeometry(QgsGeometry.fromPolyline([startpoint,endpoint]))
    self.pr.addFeatures([attributelist])
    self.layer.updateExtents()

But all I get is:

  TypeError: QgsVectorDataProvider.addFeatures(list-of-QgsFeature) ->
  (bool, list-of-QgsFeature): argument 1 has unexpected type 'list'

I'm using QGIS 2.8.1


texts: Opening and Viewing "Additional Files" in ArcGIS Layer Package

texts:
In ArcGIS I've created a layer package with several shapefiles (that include symbology and labels) so I can send it to a coworker. In addition, I've added a spreadsheet document under the "Additional Files" tab (see image below) so my co-worker can cross check it with the attribute table.

My question is, how do my co-worker open or view that spreadsheet from the layer package? When I opened the layer package in ArcMap there is no mentioning of the attachment in any of the menus.

I've researched ESRI's help pages on the subject and none of them specify how to open the additional files, only how to add them.


texts: Obtaining shape or width of road in GeoJSON format?

texts:
I'm new to openstreet maps and the geojson format.

I'm drawing buildings and roads using GeoJson. A Building is represented as a Polygon, thus i can create the proper geometry as needed, but a Road is represented as a LineString, thus i only get to draw them as lines, and i lack information about the width of it.

Drawing roads with different widths from GeoJSON using C#? seems to imply there's a way to get the width of a road using GeoJson (and also, the look of the roads as shapes instead of lines as the link's picture is what i'm looking forward in my project).

(image borrowed from above link's)

It's there in the GeoJson format data for road width/shape I'm failing to find in the GeoJson file ?

If not, it's there any sort of convention for road width i can use to draw my roads?, like estimated size

To be clear, i'm not trying to get the very exact shape of any given road, but i would like to draw roads as shapes just as the picture when possible, or at the very least understand how i can come close to achieve that.


texts: GeoServer property desciptor is null did you request a property that does not exist - problem with raster in group

texts:
In GeoServer 2.9.1 I have thematic groups configured. Each group consists of a series of vector layers and each of them has at least one raster. In version 2.9.1 cql filter deals with them perfectly and applies the filter only to vector layers, while in version 2.18 it turns out that it returns an error = property descriptor is null did you request a property that does not exist?

Is it a geoserver problem or a question of additional settings?


texts: GetFeatureInfo for more than one layer using OpenLayers 3?

texts:
I am trying to work out how to get the feature information onclick for more than one layer, apart from variable etc my code is very simmilar to the openlayers example (http://openlayers.org/en/latest/examples/getfeatureinfo-tile.html)

map.on('singleclick', function(evt) {
    document.getElementById('info').innerHTML = '';
    var viewResolution = /** @type {number} */ (view.getResolution());
    var url = wmsSource.getGetFeatureInfoUrl(
        evt.coordinate, viewResolution, 'EPSG:3857',
        {'INFO_FORMAT': 'text/html'});
    if (url) {
      document.getElementById('info').innerHTML =
          '<iframe seamless src="' + url + '"></iframe>';
    }
});

For example instead of just querying the wmsSource layer/variable I would like to query more than one?


texts: Tile Mismatch in QGis Server and WorldWind

texts:
I'm using a QGis Server to provide WMS layers to WorldWind 2.0.0.

I have a raster layer which appears correctly in the QGIS Desktop viewer, but as soon as I serve it to WorldWind and pick it up the tiles are around the wrong way and if I zoom in the land masses move about. I have checked that I'm serving the layer as WGS84 which WorldWind uses to map to the ellipsoid. I'm obviously doing something rather silly, does anyone have any suggestions what could be causing this. I attach a screen shot of the WorldWind rendered map for clarification.


texts: How to call gdaltransform with its input?

texts:
I want to convert between different coordinates. If I'm in the command line, it's as simple as:

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32616
-122 46

And I receive the output

-2193988.77788563 5724599.39928024 0

I'm trying to use this in a python program, and it seems like I need the input in the same command. Something like this (although it doesn't work)

gdaltransform -s_srs EPSG:4326 -122 46 -t_srs EPSG:32616

In my python program, I have tried to make separate calls:

command = "gdaltransform -s_srs EPSG:4326 -t_srs EPSG:" + data['crs']
os.system(command)
command = data['lon'] + ' ' + data['lat']
output = subprocess.check_output(command, shell=True)

but get an error:

CalledProcessError: Command '-122 46' returned non-zero exit status 2

I'm not sure of the correct way of doing this.


texts: I have a string of values (over 10,000) I want to label them but I don't want to label every single one. I need to label every fifth number.


texts: Where exactly are qgis.gui and qgis.core located?

texts:
I'm trying to create a standalone QGIS app using PyQGIS. I'll be deploying it by copying the install folder of QGIS and source folder of my app, and so for the testing purposes I copied QGIS install folder (1.8.0) to some other computer without QGIS installation to begin with; then I tried to test the qgis.gui and qgis.core imports, but no matter what PYTHONPATH (or PATH) I set, I keep getting ImportError  "no module named ...".

What's weird is, if I modify qgis.bat to this new location and run it, QGIS opens without a problem and all imports work through console inside QGIS. Then I simply duplicated qgis.bat and removed the last line which runs QGIS then added the line which runs my application instead. For some reason, imports are still unavailable, so I'm wondering where are actual Python binding files for qgis_gui.dll and qgis_core.dll because I can't find them anywhere inside QGIS install directory.

My .bat looks like this:

@echo off
SET OSGEO4W_ROOT=c:\IMS\Qgis
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat
call "%OSGEO4W_ROOT%"\apps\grass\grass-6.4.2\etc\env.bat
@echo on
PATH %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin;%OSGEO4W_ROOT%\apps\grass\grass-6.4.2\lib;%OSGEO4W_ROOT%\apps\qgis\python\qgis;
SET PYTHONPATH = c:\IMS\QGis\apps\qgis\bin\
python c:\IMS\IMS\ims.py


texts: Shapes shift over Bing images on QGIS when scolling the screen

texts:
Using QGIS 3.4.3, the shape files created don't fit with Bing aerial plugin only on the right side of the screen, on the other hand, for the same screen, the left side looks ok. When I move the non coincident point (for instance a cross road digitized shape and the same point on bing aerial on background)located at the right side of the screen to the left of the screen, it becomes coincident.
Can anyone explain me why this is happening and how to avoid that?

EDIT
As required by lambertj, here it goes the screenshot:


texts: QGIS Polygon Shapefile to PostgreSQL (PostGIS)

texts:
I am working on getting the polygon data of different areas loaded into PostgreSQL so that I can do searches within those areas based on address points.
I have downloaded postcode data from: https://datashare.is.ed.ac.uk/handle/10283/2597

It looks good and have loaded one of them on QGIS:

I have connected my PostgreSQL to QGIS through PostGIS connection.
I have PostGIS extension present in the database when looking at Extensions within database.
When trying to import to database get error:

Not really sure what could be wrong, never worked with QGIS before. What is the proper way to get this data including any underlying fields with polygons converted to PostgreSQL data set, so I can do a lookup based on a field such as postcode area CR5 7 and get the polygon area?

My final goal is to have a search mechanism, so you can specify search within e.g. Manchester and it will only show objects within the city. Or e.g. SW20 and get all within that post code district.


texts:
I am using QGIS 3.10.2-A Coruña. I have a UK County Boundaries shapefile layer and a custom Grid (gpkg) layer with a fixed distance of 1.5 miles between intersections in both planes (eastings & northings).
I need to:

select one county at a time;
select squares in the grid which partly or wholly overlay that county, to create a new layer.

Is this even possible?


texts:
I am trying to make a webpage which can post the possible requests and services from a spatial server, like geoserver can.

However I can only find owslib that can provide me the details.
An example is:

[operation.name for operation in wfs11.operations]

that prints out

  ['GetCapabilities', 'DescribeFeatureType', 'GetFeature', 'GetGmlObject']

from docs

another example I can think of is from the xml of a geoserver-instance:

...geoserver/wfs?request=GetCapabilities"/>

where the request=getcapabilities are telling what services is supported.
This can be compared to this example:

<ows:Operation name="GetCapabilities">

That states the available services in another way.

From the two latest examples I conclude that it would be a very large task to simply scrape the information from the xml, the task that owslib can do.
I am restrained from using owslib in a web application because of webhotel restrictions, so I'm looking for another way of acquiring spatial server metadata.


texts: How do I import a NMEA file into QGIS 2.8?

texts: I have some NMEA data recorded in a file - this is not a realtime application. How do I import that into QGIS 2.8?


texts: Converting an LAS object back to a LAS files (.las) using the lidR R package

texts:
I'm using the R package lidR to classify noise in a LiDAR dataset, but I would like to save the final LAS object as a LAS file to open it on another software.
How can I accomplish that?
Note that I'm not trying to convert the LAS object into a point shapefile, but back to a LAS file.
Moreover, my project area has more than 1,000 LAS tiles, should I use a loop to classify noise on all files or read them as LAScatalog and run the tool classify_noise() on the LAScatalog?
If the latter, how can I save a LAScatalog as multiple LAS tiles (with exactly the same number and dimensions of the input LAS files)?


texts: Writing a table (no geometry) to GeoPackage using R

texts:
I'm trying to do something so incredibly simple: write a table (with no geometry) to a GeoPackage. I can store tables in a GeoPackage using QGIS, so it stands to reason I'd be able to write one from R. So far no luck. Here's an example:

library(sf)

t1 <- data.frame(numbers=c(1:5), letters=c(LETTERS[1:5]))

st_write(t1, dsn="NewGeopackage.gpkg", layer="t1")

I've tried converting the table to other formats (using st_as_sf, for example) but I always get the same error:

Error in st_sf(x, ..., agr = agr, sf_column_name = sf_column_name) :
  no simple features geometry column present

Any idea how to do this?


texts: Why do different USGS Topo maps show different elevations for the same control point?

texts:
So I have been investigating the area around Jerome and Clarkdale Arizona using different USGS Topographic map sources of the same scale, but different time periods, and find that in some cases, the same control point has a difference in elevation in one source than in another.

In the case of two topographical maps published by USGS, one in the 1900s and the other a number of years later, every control point I have examined (though not by any means all) the difference in values for each control point varies by the same amount between the maps at every point I examined (which I admit is by no means all).

This suggests to me some type of systemic error, perhaps recalculation of mean sea level. Considering another pair of maps, the difference in elevation given on the two maps for the same control point varies by a different amount than the difference in elevation given on the two maps for another control point. In the latter case, the difference in the elevations ranges from 0 to a couple of meters.

Where can I find an explanation of what the differences between these maps stem from, and how do I decide which is the most accurate elevation?


texts:
I'm trying to use a custom coordinate transformation in FME 2016.

I have defined the following transformation  in the file "MyCoordsSysDefs.fme":

XFORM_DEF Test \
SRC_DTM ED50 \
TRG_DTM ETRS89 \
DESC_NM "7 parametros" \
GROUP "Custom" \
INVERSE Yes \
MAX_ITR 8 \
CNVRG_VAL 1e-09 \
ERROR_VAL 1e-06 \
METHOD 7PARAMETER \
DELTA_X -130.5090 \
DELTA_Y -89.9490 \
DELTA_Z -145.9030 \
BWSCALE 5.2690 \
ROT_X 3.1126 \
ROT_Y -0.4718 \
ROT_Z 3.3181

I can't find how to use this in the offical documentation . How can I tell FME to use this transformation instead of the default one?


texts:
Can I rearrange the elements within the Layout toolbar?

I want to put the Align Edge tool elsewhere, away from Shared features, like at the end of the Toolbar

I have tried customize, and I cannot seem to be able to reorder this Toolbar.


texts: Unable to install the PostGIS extension to PostgreSQL on Ubuntu 12.04 LTS

texts:
I just installed Linux Ubuntu 12.04 LTS (on a virtual machine). I then successfully installed the PostgreSQL 8.4 database management system. I am now trying to install PostGIS, which is an extension to PostgreSQL used for Geographical Information Systems capabilities. I followed the instructions for installing PostGIS on this webpage:
http://postgis.net/docs/manual-2.0/postgis_installation.html

Section 2.4.1 ("Configuration") on this webpage tells me to type "./configure".
I did this. The output (about 100 lines of mostly incomprehensible statements) included about one dozen instances of the following error message:

You need to install postgresql-server-dev-X.Y for building a server-side extension or libpq-dev for building a client-side application.

At the end of the output, the following lines were outputted:

checking libpq-fe.h usability... no
checking libpq-fe.h presence... no
checking for libpq-fe.h... no
configure: error: could not find libpq-fe.h

I have no idea what caused these errors or how to fix them. I would appreciate any assistance or pointers.

NOTE: the above-linked webpage includes a section 2.2 ("Requirements") which lists eight required software packages:

PostgreSQL 8.4
GNU C compiler (gcc)
GNU Make (gmake or make)
Proj4 reprojection library, version 4.6.0 or greater
GEOS geometry library, version 3.2.2 or greater
LibXML2, version 2.5.x or higher
JSON-C, version 0.9 or higher
GDAL, version 1.6 or higher

I just installed PostgreSQL 8.4, and I am pretty sure that the GNU C compiler is also installed on my computer. But I do not know whether I have the other six software packages installed or not (or if they are installed, whether they are properly configured, etc.).


texts: Filtering multiple values (more than 2) using XML filter for WFS?

texts:
I'm having a problem using WFS XML filter. I'm trying to read only 6 features from entire dataset and I get the following error:

I wrote query in SQL and need to convert it to XML:

"SAV_ID" IN (21,23,25,55,56,88)

Interestingly when I'm quering only 2 of thoses features, I get no errors. XML filter (working one) for 2 features is:

    <ogc:Filter>
    <ogc:Or>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>21</ogc:Literal>
        </ogc:PropertyIsEqualTo>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>23</ogc:Literal>
        </ogc:PropertyIsEqualTo>
    </ogc:Or>
</ogc:Filter>

I tried the following variants (for more than 2 features):

(1) (gives error)

    <ogc:Filter>
    <ogc:Or>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>21</ogc:Literal>
        </ogc:PropertyIsEqualTo>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>23</ogc:Literal>
        </ogc:PropertyIsEqualTo>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>25</ogc:Literal>
        </ogc:PropertyIsEqualTo>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>55</ogc:Literal>
        </ogc:PropertyIsEqualTo>
    </ogc:Or>
</ogc:Filter>

(2) (gives error)

    <ogc:Filter>
    <ogc:Or>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>21</ogc:Literal>
        </ogc:PropertyIsEqualTo>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>23</ogc:Literal>
        </ogc:PropertyIsEqualTo>
        <ogc:PropertyIsEqualTo>
    </ogc:Or>
    <ogc:Or>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>25</ogc:Literal>
        </ogc:PropertyIsEqualTo>
        <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>55</ogc:Literal>
        </ogc:PropertyIsEqualTo>
    </ogc:Or>
</ogc:Filter>

(3)

    <ogc:Filter>
        <ogc:PropertyIsEqualTo>
        <ogc:Function name "in">
            <ogc:PropertyName>SAV_ID</ogc:PropertyName>
            <ogc:Literal>21</ogc:Literal>
            <ogc:Literal>23</ogc:Literal>
            <ogc:Literal>25</ogc:Literal>
            <ogc:Literal>55</ogc:Literal>
        </ogc:Function name "in">
        </ogc:PropertyIsEqualTo>
</ogc:Filter>

Last one (3) downloads all the features. Any help on that?


texts: Create rectangle by entering width & height numerically in QGIS?

texts:
I need a quick way to create rectangles with set with & height in meters. Is there a straightforward way to do this? I am well aware that there are a number of ways to achieve this, but all that I have found (apart from the QAD plugin) are cumbersome workarounds. The keyword here is quick and easy.

Please note that the size of the rectangle is the priority. It will be moved around a bit until we decide on its final coordinates. So I would rather avoid solutions to my problem where the rectangle is drawn from coordinates.

Ideally I would like to click anywhere, enter height and with in meters and then move the polygon to its final resting place. I have used the QAD plugins rectangle tool to do this in the past but it seems to be currently broken.


texts: QGIS Print composer: Hiding locked Items in the Item Property items menu (not in in the layout itself) or putting another safeguard on it

texts: To make sure certain items cannot be worked on in layout (not being clicked on and edited) they are set to locked. While at first this would safeguard certain text and labels, for other people that have little bit of knowledge of QGIS they still can work around this and  access these in the Item Properties menu (for example the watermark label) and edit the text there. To make sure these items just cannot be removed or not interacted by I want to "hide" these items in the Item Properties menu, or at least put some kind of extra safeguard on that only the project author (aka me) can access it.


texts: ogr2ogr - Clipping postgis and creating geojson

texts:
I'm trying to create a geojson file data based on data from my local postgis database using ogr2ogr. I'm trying to clip to a specific area and when I start to run the command, everything starts fine, but I then get an ERROR 6: GeoJSON driver doesn't support creating more than one layer and I'm not quite sure how to solve it.

Do I need to select only one column in the SELECT statement?

Command:

ogr2ogr -f GeoJSON ./test-buildings.geojson PG:"host=localhost dbname=gis user=*USERNAME* port=5432" -clipdst -84.419932, 33.992761, -84.371342, 34.018029 -clipdstsql 'SELECT way, "building" as type FROM planet_osm_polygon ORDER BY ST_YMin(ST_Envelope(way)) DESC' -progress


texts: Calculating global horizontal irradiance (GHI) from surface solar radiation downwards from the ERA5 data

texts:
I have some hourly data of surface solar radiation downwards (SSRD) measured in J/m^2 from the ERA5 dataset (ECMWF - European Centre for Medium-Range Weather Forecasts).
How can I convert these to global horizontal irradiance (GHI) in kWh/m^2? Is it just simply SSRD/3600000?
I know the formula for GHI is:

where:
DNI = direct normal irradiance
DHI = diffuse horizontal irradiance
z = solar zenith angle
Is SSRD sufficient for calculating GHI or do I need other data?


texts: How to deploy custom .dll commands in ArcMap 9.3.2?

texts: In ArcMAp 9.3.2, I am trying to bring in an acreage calculator .dll buttom that I have used on previous versions.  The .dll doesn't show up on the list of available custom addons.  Is this a registration problem?  I'm at a loss of what could be causing the .dll to not be available on the list.


texts: QGIS 3.10: Evaluate default values (data source preference) not staying checked on?

texts:
Anyone having trouble with "evaluate default values" in QGIS 3.10.1 preferences, under "Data Sources"?

When I click the check box and then "OK", I go back to preferences it's unchecked! And default values are not being evaluated... I am a QGIS novice, so not qualified to call "bug"!

When I add new features, I want the field calculations in the table to fill automatically. Is this the way to go about it?

I don't want to use virtual fields as I want the calculated fields to be part of the geopackage/shapefile. I'd like the values to update as I modify features.


texts:
I have a gray scale tiff image. I can extract the values from this tiff using numpy and make a new image with blue scale. But how can I convert it to shades of blue instead of shades of gray directly from the Tiff? Is it even possible?

Edit:
Well I have two sets of tiff images(data), one with color-palette and another with gray scale ,I could generate tiles from color-palette one using generate_tiles.py,
 and show them on a map, but with gray scale one I can't do that, so I am wondering if I could change this gray scale tiff to color-palette one (with custom colours) so I can generate tiles and visualize the data.


texts: Filtering features using geometric attributes of polygons in ArcGIS Desktop?

texts:
Given a polygon featureclass containing many polygons of different sizes and shapes, which geometric attributes (area, perimeter, etc) could be used to filter the polygons? It's probably best to ask this in two parts:

Which geometric attributes do polygons possess?
How would you calculate these attributes? Geoprocessing tools/ArcObjects preferred, but I would be satisfied with conceptual calculations which could then be coded.

As an example, I am considering calculating the maximum length of a polygon by breaking the edge of the polygon into points at a small interval and then calculating the distance between all points - the longest distance would be the maximum length of that polygon.

I'm hoping to use these attributes to allow interactive user control of polygon filtering.


texts:
I am preparing a routing function with pgrouting with the travel model in public transport (bus).

I already have all the stops related to the table of ways (downloaded by osm2pgrouting) I have my table with the stops and the source and target of the segment of the line in which they are.

But my question is how to prepare a query where you can find the closest stop at the beginning of the route and at the end of the route and also start at that stop that will be the closest and end at the stop closest to the destination.

The only thing I have achieved is to represent the bus route, going through all its stops, but of course, when I launch my query it represents the entire bus route, and I don't want that. As I said before, what I want is that the beginning and end of the route be the segment of the line where the stops are, but not that the entire bus line represents me.


texts: Troubleshooting import arcpy gives ImportError: No module named numpy?

texts:
While I try to load the arcpy module in arcGIS 10.3.1 I get the following error message:

    Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 24, in <module>
    from arcpy.toolbox import *
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\toolbox.py", line 359, in <module>
    from .management import Graph, GraphTemplate
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 26, in <module>
    import _management
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\_management.py", line 14, in <module>
    import _graph
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\_graph.py", line 27, in <module>
    import numpy
ImportError: No module named numpy

How to fix this?


texts:
I need help choosing the appropriate Metadata standard for a GDB for a local government. what are your recommendations and why ? what are the main differences between ISO 19139, ISO 19115 and FGDC.

Thanks Geog


texts: How to rename field names in a shapefile from the commandline?

texts: I need to rename the field names (attributes) in a shapefile. Is there an easy way to do this from the command line?


texts: Has anyone successfully made a FeatureLayer from a view in DB2 and ArcSDE 10?

texts:
Lately it seems like every new ESRI feature fails for me - I must be cursed. I'm testing the use of a FeatureLayer to see if helps display and query features using the JavaScript API 2.0.

This layer was registered with ArcSDE from a view. The data is in a DB2 database. We have used the view for a long time in a MapService on our server.

I copied the MapService MXD to a new file and published it as a FeatureLayer. As a MapService it shows up in ArcCatalog and can be previewed. It also shows up in ArcGIS Service Manager. The MapService doesn't show up in the ArcGIS Services Directory page. If I enter the URL in a browser, it returns an error. If I try to add the layer to a JavaScript API map, it raises an error "missing } in XML expression" in Firebug for the script arcgis?v=2.0, line 52.

Is there anybody who's used a FeatureLayer successfully? Any pointers or hints?


texts: Does a SRTM DEM image show details such as fields?

texts:
I'm working through the QGIS Training Manual where I have a step to get ...

  ... raster images (SRTM DEM) which cover the region you have
  selected for your course.

  SRTM DEM can be downloaded from the
  CGIAR-CGI: http://srtm.csi.cgiar.org/

  You’ll need images which cover
  the entire region you have chosen to use.

It then shows an example image (Image 1) for Swellendam, Cape Town:

However, the only image I can find for my region (UK) on CGIAR-CGI is similar to the image shown below (Image 2), which only seems to show digital elevation and not any other detail such as houses, roads or fields like 'Image 1'.

Question:

Does a SRTM DEM raster image exist for the UK region on CGIAR-CGI that has the same level of detail as 'Image 1'?


texts: Importing/Exporting parameters from a text file into a batch process

texts:
I usually work my batches in Excel or TextPad, and then copy/paste it into ArcGIS batch UI.

Is it possible to do the same in QGIS? The tutorials I've seen show how to fill the interface line by line, only. I'm using QGIS 2.10.1.

BTW, right now I'm trying to merge all seven bands of Landsat imagery into one multiband image file. I have 152 scenes with 7 bands each.


texts:
Trying to extract MODIS Deep Blue Aerosols from a number of years.

I have managed to open the HDF files in QGIS and ArcMap; however, have failed to reproject out of the sinusoidal grid so that the UK matches up with reference imagery.

I have tried MRT, which gave me an error every time I loaded a 04_L2 file into as follows

  Error in Module: ReadHdrFile Opening Input Header File

I couldn't find an explanation except perhaps, that I needed to use MRTSwath instead.

When using MRTSwath, the file loads with no problems but only displays two bands: the Aerosol_Cldmask and the Cloud_Pixel_Distance. I imagine because these are the only two bands still in a swath rather than grid projection.

I have looked into doing the conversion through gdal, however, when getting the needed GCP's via a gdalinfo command on the subdirectory, the produced information fails to list anything to do with GCPs. I thought this might be due to the version of gdal I have, but I seem to have version 2.2.4 which is well past what I believe I need for GCPs to be produced.

Any alternatives for reprojecting?


texts: Deleting rows on dbf table with python for ModelBuilder?

texts:
Being not very familiar with python, and after carrying out unfruitful long research on internet, I submit to you my problem, hoping favorable reply.

I have a DBF table with geochemical data and I would like writing a script to be integrated in a wider model built with ModelBuilder, and that automatically delete some specific rows, specifically (please see attached example) the rows 1, 2, 4, 5 and 6, while the row 3 (Analyte symbol) is maintained. This is first.

Secondly, I would like to automatically fill in some empty cells in some coloumns (e.g., A, B, C) with the character 0 (zero). Could you please show me the way?

In fact, by deleting rows 1, 2, 4, 5 and 6, I want the row 3 to become the first row (headline) of the table (i.e., Au, Pt, pd, SiO2...). The latter contains 150 rows of data below the 6 first rows, that begins with #59101 in the provided example.

Thus, I do not want to keep just the row 3, but really deleting the rows 1, 2, 4, 5 and 6, so that the row 3 becomes the headline of 150 rows of data (as shown in the herein attached example).


texts: Spatial join of a single point to multiple polygons

texts:
I have a point layer which I want to spatially join to a polygon layer (bring in the attributes).

This works fine unless there are overlapping polygons, if so the QGIS spatial join tool will only bring back the features from the first one found (or some average). An average is not suitable for text attributes.

What I would prefer to happen is that if there is multiple polygons it can be joined to, it will, and result in a duplication of the point. Essentially the ArcGIS: join_operation = JOIN_ONE_TO_MANY

I am working in QGIS 2.8+.

Happy with a script/PyQGIS solution, as I don't think there are any ready made tools for it. But would prefer not to have a SpatiaLite solution as I want to keep it simple for the end user.

Ultimately if there is a solution to this problem that someone has already implemented I would love to see it. If not I will try to come up with the solution myself.


texts: Calculate accumulated distance with road network

texts:
I am working on a project in ArcMap and the task I am working on goes as follows:

Calculate the accumulated distance (in meters) of every pixel on the
map to Itaituba taking the road network into account (yellow vector lines).
The final output should cover only the state of Mato Grosso but as
the port of Itaituba is located in Pará, it will be necessary to use
other states area and road network for the modeling.
Use the Itaituba port as reference for destination
Resolution of the final output 100m

I have worked with the Network Analyst of ArcMap before but never done accumulated distance. Any suggestions on suitable tools to use?


texts: Geoserver 2.5.2 - How do I enable the GeoJSON bbox property?

texts:
Just started using Geoserver 2.5.2. On my old server which was Geoserver 2.1 when you ran a WFS request to return GeoJSON (or, for a quick example, clicked Layer Preview -> All Formats -> GeoJSON) the returned JSON contained a bbox property giving the bounding box for the returned results. I've written some applications that rely on this value being present but now they are all failing because the GeoJSON being returned does not conatain the bbox property.

e.g.
old:

{
"type":"FeatureCollection",
"features": [..features here...],
"crs":{"type":"EPSG",
    "properties":{"code":"4326"}
      },
    "bbox":[0,52,0,52]

}

New:

{
"type":"FeatureCollection",
"features": [..features here...],
"crs":{"type":"EPSG",
    "properties":{"code":"4326"}
      }

}

Is there an option I need to enable to encourage it to return this property or has the GeoJSON spec changed since GeoServer 2.1 to exclude it? That would be very inconvenient! Can I get it back?


texts: The buffer is covering unnecessary area of the polygon it was made from in QGIS

texts:
Consider a polygon that represents the area of New York State

I made a 125KM buffer from this layer

My problem is that the new buffer layer is covering the source polygon and I need only the outer borderlines of the 125KM radius. I need it as a new line layer for future work, so just changing the fill of the polygon into 100% transparent will not work .

How can I display only a line layer that is located in a distance of 125KM from the borderline of the New York State polygon?