texts: Specifying join on multiple datasets in a WFS URL request (CQL/Filter)?

texts:
I’m using a WFS and trying to join two/three datasets together in a URL query.
Two of the datasets are linked via an aspatial table in the middle like so:
| Table 1           |        | Table 2   |       | Table 3            |   | PFI | Item | Geom |    >   | PFI | AFI |   >   | AFI | Item2 | Geom | |  1  | House| foo  |        |  1  | 45  |       |  45 | Street| bar  | |  2  | House| foo  |        |  2  | 78  |       |  78 | Road  | bar  |
I need to join the tables in a WFS URL request and return all results that match the filter. I keep getting the error "Join query must specify a filter" when trying to specify how the join matches.

In SQL I would do something like:
SELECT * FROM TABLE1
LEFT JOIN TABLE2
ON TABLE1.PFI = TABLE2.PFI
WHERE TABLE1.PFI = 1

My main question is how do I achieve this in the WFS URL request?
WFS I'm using:
http://services.land.vic.gov.au/catalogue/publicproxy/guest/dv_geoserver/wfs?request=GetCapabilities&service=WFS&
I've attempted to use both CQL_Filter= and Filter=:
Example:
http://services.land.vic.gov.au/catalogue/publicproxy/guest/dv_geoserver/wfs?service=WFS&version=2.0.0&request=GetFeature&typenames=datavic:VMADD_ADDRESS,datavic:VMPROP_PARCEL_PROPERTY&CQL_FILTER=datavic:VMADD_ADDRESS/PROPERTY_PFI='41116831';datavic:VMADD_ADDRESS/PROPERTY_PFI=datavic:VMPROP_PARCEL_PROPERTYS/PROPERTY_PFI&

XML &Filter= variation:
<Filter>
  <And>
    <PropertyIsEqualTo>
      <PropertyName>/datavic:VMADD_ADDRESS</PropertyName>
      <Literal>41116831</Literal>
    </PropertyIsEqualTo>
    <PropertyIsEqualTo>
      <PropertyName>/datavic:VMADD_ADDRESS/PROPERTY_PFI</PropertyName>
      <PropertyName>/datavic:VMPROP_PARCEL_PROPERTY/PROPERTY_PFI</PropertyName>
    </PropertyIsEqualTo>
  </And>
</Filter>

EDIT:
The below query is the closest I've gotten to getting it working, but receives a different error:
http://services.land.vic.gov.au/catalogue/publicproxy/guest/dv_geoserver/wfs?service=wfs&version=2.0.0&request=getfeature&typenames=datavic:VMADD_ADDRESS,datavic:VMPROP_PARCEL_PROPERTY&CQL_FILTER=(datavic:VMADD_ADDRESS.PROPERTY_PFI=datavic:VMPROP_PARCEL_PROPERTY.PROPERTY_PFI)

ERROR:
<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="2.0.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://services.land.vic.gov.au/catalogue/publicproxy/guest/dv_geoserver/schemas/ows/1.1.0/owsAll.xsd">
<ows:Exception exceptionCode="OperationProcessingFailed" locator="GetFeature">
<ows:ExceptionText>
Error occurred getting features Cannot query this feature source with ADDRESS since it serves only VMADD_ADDRESS Cannot query this feature source with ADDRESS since it serves only VMADD_ADDRESS
</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>


texts:
I have a series of fields to check for a certain value, then calculate other fields with using UpdateCursor.

Fields are: SPECIES_CD_1, SPECIES_CD_2, etc.. up to SPECIES_CD_99

Here is the general idea (some pseudocode in here):

rows = arcpy.UpdateCursor(FC)
speciesLevel = "SPECIES_CD_"
for x in range (2,99):
    speciesX = speciesLevel + str(x) #this is a string
    speciesUnicode = speciesX.decode #Decodes the speciesX string to Unicode

    for row in rows:
    if row.speciesUnicode == 'P%':
        [do a whole bunch of things]
    elif row.speciesUnicode == None:
        [do other things]
    else:
        [do other things]

Of course the cursor cannot identify a row named speciesUnicode, so the cursor fails to run. I want speciesUnicode to act as a variable that subs in a field name. Is this possible somehow?


texts: Leaflet / Mapbox.js: incorrect Popup offset using divIcon with CSS class

texts:
I'm having an error which I can't seem to find elsewhere online with Leaflet / mapbox.js regarding the offset of a popup from a marker's custom icon which has been set with L.divIcon and using a CSS class to style:

Here's the invocation of the L.divIcon:

var icon = L.divIcon({
  'iconSize': [30, 40],
  'iconAnchor': [15, 40],
  'className': 'mapview-marker-icon'
})

The CSS class also has the width and height set:

.mapview-marker-icon {
  box-sizing: border-box;
  display: block;
  width: 30px;
  height: 40px;
  padding: 0;
  margin: -40px 0 0 -15px;
  background-color: transparent;
  border-width: 0;
  transform-origin: bottom center;
}

One thing I've found is that if I don't set the iconSize and the iconAnchor in the JS, then the marker's icon is overridden by mapbox.js / Leaflets's CSS styles applied directly to the element (mapbox.js documentation suggests that iconSize and iconAnchor can be set via CSS classes and don't need to be set in JS):

(note: I was also playing with display props to block and inline-block to see if that changed anything. It didn't.)

When I use L.mapbox.marker.icon the popup offset is correct:

var icon = L.mapbox.marker.icon({
    'marker-size': 'large',
    'marker-color': '#000'
})

Is there something I'm missing on the CSS or JS side of things? The popup is applied using L.Marker.bindPopup().


texts: Calculation of surface of external walls of buildings from a cadastral footprint shapefile

texts:
I need to calculate the surface (m2) of external walls of all the buildings in a town.
By external walls I mean all the walls of a building or parts of a wall which are not in contact with another building.
I have a 2D shapefile of polygons with the cadastral footprint of the buildings, which includes information of the building’s height in the attribute table.

I don’t really know how can I do it… My idea is to convert the 2D file to 3D and use the 3D shapefile to distinguish the external walls from the shared walls (parts of a building that are shared with other buildings) of all buildings, and from there, calculate the total area but I’m not sure whether this can be done with ArcScene.

Anybody knows how could it be done and if there is an automate way to do it?


texts:
I imported two layers to my Project (one is a shapefile from natural earth, the other is built out of lon/lat Points in a csv-data-file) QGIS fits the two graphics perfect together. Now my goal is to export this as an svg-file. This works fine, too.

Is there a way to store id-names for the paths and objects in the svg-file?

The ids should be based on the underlying attribute-tables of each layer. I have to do this automatically because of the amount of data-points. I read something about layer splitting. But I cannot even find the names of layers in my svg output file so I guess this will not help me. If there is no way to manipulate the svg directly, can it be done with a workaround? (I think about putting shapefiles with all informations together and transforming it to svg in another way?)

I am using QGIS 2.18.17 Las Palmas

I want to identify the polygons in the underlying svg data, like i have a polygon in my Image, and I open the svg-file with a text Editor, I want to have <path id="USA"><\path>. I want to have the Information stored in the metadata of the vector-graphic, but I don't want to see the data in the picture. Afterwards I want to adjust the file to fit into html and therefore I have to be able to identify single parts of the picture.


texts:
I am having an error in tilemill with postgresql and the Openstreetmap-carto style that I am having trouble figureing out the reason for. Everything is working perfectly up until I reach zoom level 16. When tilemill is zoomed to level 16 I am receiving this error.

I am not sure what the cause of it is, or why it is trying to find that attribute when it reaches zoom 16 and not at any other level. The reason this is an issue is I was looking to export the mbtiles at level 17.
I attached the mml code below.

http://pastebin.com/x7Sn0ks5


texts: Showing two sets of points from two different fields in QGIS?

texts:
I made a union of two adjacent shp (points) files with the same theme (Geochemistry) so the resulting attribute table gave duplicated fields, for example: ID, As_ppm, Ag_ppm, Pb_ppm, ... , ID_2, As_ppm_2, Ag_ppm_2, Pb_ppm_2,...
And the elements of the first shp give 'NULL' in the fields of the second shp and vice-versa, As in the screenshot.

So my goal is to show the combined points for a speciffic element, say Arsenic, so I have tried, without success, expressions like As_ppm || As_ppm_2 , As_ppm + As_ppm_2, combine(As_ppm ,As_ppm_2).


texts:
I have the below SLD and I I need to filter based on multiple (several dozen) that do not meet the criteria:

<PropertyIsNotEqualTo>
<PropertyIsLessThan>
<PropertyIsLessThanOrEqualTo>
<PropertyIsGreaterThan>
<PropertyIsGreaterThanOrEqualTo>

And I can not edit the DB table so I am left with <ogc:PropertyIsEqualTo>.  Is there a way to filter with multiple Literal statments within the function? if I add a second Literal call i get the error:

line 20: cvc-complex-type.2.4.d: Invalid content was found starting with element 'ogc:Literal'. No child element is expected at this point.

I can not duplicate the code because it creates a new legend category for each filter.  Using this post I have been researching functions but I still have been unable to filter with multiple values.

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd"
  xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

  <NamedLayer>
    <Name>test</Name>
    <UserStyle>
      <Title>Other</Title>
      <FeatureTypeStyle>
        <Rule>
          <Title>Other</Title>
          <ogc:Filter>
              <ogc:PropertyIsEqualTo>
              <ogc:Function name="in">
              <ogc:PropertyName>type</ogc:PropertyName>
              </ogc:Function>
              <ogc:Literal>0</ogc:Literal>
              </ogc:PropertyIsEqualTo>
              </ogc:Filter>
          <PointSymbolizer>
            <Graphic>
              <Mark>
                <WellKnownName>arrow</WellKnownName>
                <Fill>
                  <CssParameter name="fill">#0033cc</CssParameter>
                </Fill>
              </Mark>
              <Size>12</Size>
              <Rotation>
                <ogc:PropertyName>dir</ogc:PropertyName>
              </Rotation>
            </Graphic>
          </PointSymbolizer>
        </Rule>
     </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: Why is the dijit.popup not working correctly?

texts:
I take the following steps:

Click 2 overlapping features from 2 different featLayers.
Popup auto-selects and highlights map.infoWindow.features[first object]
The popup content is actually from map.infoWindow.features[second object]

After you navigate to the second graphic and back, the popup corrects it's self??

here's an example using esri data:
http://jsfiddle.net/92Q2G/1/

Any advice?
Thanks.


texts: Cutting features to do opposite of clipping in ArcGIS Desktop?

texts:
What´s the easiest/smartest way of getting the opposite output of the clipping tool?

I just need the area of feature A1 which is not covered by feature A2 in a new Layer/File.

Now, I get the result by union´ing both features and deleting the "wrong" parts by hand.

I'm using ArcGIS Desktop 10.


texts: Selecting service (WMS / WFS) through GeoServer REST API?

texts:
Is it possible to select a service (WMS / WFS) through the Geoserver REST-API?

I know how to add a workspace through the REST-API, like this

curl -v -u admin:geoserver -XPOST -H "Content-type: text/xml" -d "<workspace><name>workspacename</name></workspace>" %resturl%/workspaces

But I would like to add a service (WMS or WFS or both) through the REST-API. If this is possible I don't have to go the the webadmin-page the select the service manually.
Does someone know if this is possible?


texts:
I am trying to create a graphic demonstrating the difference between an idealized perfect sphere and the Earth's reference ellipsoid. (I know this is also an approximation, but nevertheless.) To this end I want to also exaggerate the eccentricity of the ellipsoid.

https://en.wikipedia.org/wiki/Latitude

Anyway, should I simply scale a sphere vertically such that the vertical axis is shorter than the horizontal axis? If I have a map of the Earth wrapped around the sphere, is squashing the sphere vertically sufficient? Or do I need to convert between geocentric and geodetic latitude first?

What's the proper way to handle this situation?


texts: WMS returns but WFS does not return same SRS and BBOX values

texts:
I am using WMS and WFS services on GeoServer. My WMS query returns the image result for features:

http://localhost:8080/geoserver/wms?service=WMS&request=GetMap&layers=Gis:Poi&version=1.1.1&useCors=false&height=256&width=256&srs=EPSG:3067&&BBOX=316600,6838986,316605,6838991&

But I am using same BBOX and SRS values for WFS service but does not return any value.

http://localhost:8080/geoserver/wfs?request=GetFeature&version=1.1.0&srsName=EPSG:3067&typeName=Gis:Poi&BBOX=316600,6838986,316605,6838991&service=WFS&


texts: What is the elevation unit of the Natural Earth's Shaded Relief raster?

texts:
I need to use this raster data to extract elevation for points and polygons, but I am really confused about the unit of the pixel values

http://www.naturalearthdata.com/downloads/10m-raster-data/10m-shaded-relief/

The pixel values are between 49 and 252, and oceans have higher pixel values than some land areas. How should I interpret the pixel values and convert them to meters?


texts: Figuring out map scale when you know the spatial resolution of a Satellite Image?

texts:
I'm interested in figuring out how map scale and satellite pixel size (spatial resolution) are related.

For example the TM or ETM+ has roughly 30m ground resolution, then, what could be the finest scale of the extracted map?

Is there any formula or rule of thumb for converting the satellite pixel size to map scale?


texts: Arc Hydro 2.0 for use with ArcGIS 10.1 (or later) for Desktop?

texts:
I have ArcGIS 10.1 for Desktop and the toolbox of Arc Hydro 1.1 beta.

I want to find and donwload the Arc Hydro 2.0 for use with ArcGIS 10.1 for Desktop but, after research on the net, I cannot find it.

I think that the new version is now released so where can I download it?


texts: Unable to create multi-user system using GeoServer

texts:
I am unable to secure my GeoServer installation and grant read and/or write access to different roles. On the most basic level, I would like to have a VIEWER and an EDITOR role that I can apply to different workspaces. Along with this, I would like to control how each role accesses WFS -- readers should be able to view while editors should be able to edit in QGIS.

After a thorough read of the GeoServer manual and this tutorial, I still cannot get multi-user access to work as expected. All users, regardless of data access rules, cannot access WFS layers once the following rules are in place.

Note that I am using a fresh install of GeoServer version 2.10.0 on an Ubuntu Server version 16.04 running Apache2 and Tomcat7.

Steps taken:

Create workspace called "test"
Create a Postgis store using admin credentials (just for testing, not production)
Publish a single layer.
Create two roles: EDITORS and VIEWER
Create two users: test_editor and test_viewer. Give test_editor the EDITORS role and test_viewer the VIEWER role.
Set up data security as follows:

*.*.r  *
*.*.w   *
test.*.r    VIEWER,EDITORS
test.*.w    EDITORS
Create the following service rules:

*.*    *
wfs.*   EDITORS
Launch QGIS 2.16.3 and add the WFS 1.0.0 service capabilities (e.g. http://myserverip:8080/geoserver/ows?service=WFS&version=1.0.0&request=GetCapabilities&) AND the test_viewer or test_editor login credentials.
With the test_editor credentials, the expected result would be the ability to edit and save the layer, and with the test_viewer credentials the user should be able to simply view the WFS layers in QGIS without the ability to edit. However, trying to get capabilities causes a popup error in QGIS: "Unexpected end of file". QGIS error log reads:

Download of capabilities failed: Error downloading http://myserverip:8080/geoserver/ows?version=1.0.0&&SERVICE=WFS&REQUEST=GetCapabilities&VERSION=1.0.0 - server replied: Forbidden

GeoServer log reads:

2016-11-17 13:59:22,480 WARN [wicket.Localizer] - Tried to retrieve a localized string for a component that has not yet been added to the page. This can sometimes lead to an invalid or no localized resource returned. Make sure you are not calling Component#getString() inside your Component's constructor. Offending component: [XMLUserGroupServicePanel [Component id = dummy]]
2016-11-17 14:04:10,769 INFO [geoserver.wfs] -
Request: getServiceInfo
2016-11-17 14:07:25,217 INFO [geoserver.wfs] -
Request: getServiceInfo
2016-11-17 14:08:39,836 INFO [geoserver.wfs] -
Request: getServiceInfo

BUT... and this is where it gets weird, the credentials work perfectly for WMS. No creds, no access to WMS layers in "test" workspace. Apply creds, and authorized users have access to layers. So, there is something wrong with WFS(-T).

Some further things to note:

All default GeoServer settings are in place (except for those changes in the above steps)
WFS Service level is set to Transactional
Local settings are not being used for the "test" workspace. Global settings are used, which are set to default.
Apache, Tomcat7 and Ubuntu Server have not been reconfigured. All default settings are in place.

Why are the security settings not working as expected for WFS layers but work fine for WMS? Further, it is this a GeoServer-, QGIS-, Apache-, Tomcat 7-, or Linux-related issue? I suspect that there are some Linux permission issues or something else under the hood preventing GeoServer from working as it should.


texts: Clip heat map so it does not extend past boundaries of another layer

texts: I have the following heat map that was created from clusters of vector points. I have a second layer underneath this heat map that shows the U.S. How do we prevent the heat map from extending past the boundaries of the U.S.?


texts: Convert SpatialPoints to SpatialPointsDataFrame

texts:
Give an object x of class SpatialPoints how do I convert it to a SpatialPointsDataFrame ?

class       : SpatialPoints
features    : 50
extent      : 3187836, 3827836, -489229.3, 950770.7  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs


texts: Python Error while Publishing mxd ArcGIS 10.1


texts: QGIS GRASS: Error when converting vector to raster

texts:
I am trying to convert a vectorial layer (a grid of regular points with values, representing probabilities, associated to each point).
I tried to use QGIS, but I got this error:

ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

I checked and I don't have any probability for any point as zero.

I then tried to use GRASS (in "source of raster values" I selected "attr" and then in "Attributes" I selected the column in the file which has the probabilitites), but it does not work again:

ERROR: Column type (CHARACTER) not supported (did you mean 'labelcolumn'?)

I used "v.info" to see the information about the vectorial layer and this is the result (I don't know if it helps):

Do you have any idea why this is happenning?

Here is a preview of the dataset:

ID  PROBUSE1
499 0.032075
500 0.020963
501 0.027376
502 0.0507
503 0.058432
504 0.043781
505 0.029973
506 0.028873
507 0.026719
508 0.025188
509 0.037346
510 0.028899
511 0.021097
512 0.017613
513 0.014795
514 0.017285
515 0.019617
763 0.035415
764 0.049729
765 0.050707


texts: SNAP error encountered with Sentinel-1 SAR from Google Earth Engine

texts:
I have downloaded Sentinel-1 SAR images from GEE in VV and VH polarization.
I have two groups of images: 17 images in VV polarization and 17 images with VH polarization.

I have a model that is suppose to correct those images in SNAP:

As you can see, whenever I try to run this I get this error message:

I thought maybe the error is because I have for each image only one polarization mode because my code extract it like this. Or maybe the problem is that I haven't downloaded it directly from Sentinel, however, I have seen many researches that have used GEE to download SAR images.

This is the code I have been using:

var img = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation','VH'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .select('VH')
        .filterDate('2018-01-01','2018-05-01')
        .filterBounds(geometry)
        .map(function(image) {
          var edge = image.lt(-30.0);
          var maskedImage = image.mask().and(edge.not());
          return image.updateMask(maskedImage);
        });

var desc = img.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));

var dataset=img.size();
print(dataset)

var clippedCol=img.map(function(im){
   return im.clip(geometry);
});

var listOfImages =(clippedCol.toList(clippedCol.size()));

var listOfNumbers = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16];

for (var i in listOfNumbers) {
  var image = ee.Image(listOfImages.get(listOfNumbers[i]));
  var toexport=image.visualize({min: -30, max: 5}).addBands(image);

  // do what ever you need with image
  Map.addLayer(image,{min: -30, max: 5}, i);
  Export.image.toDrive({
  image: toexport.toFloat(),
  description: i,
  scale:10,
  crs:'EPSG:4326',
  maxPixels:1310361348,
  region:geometry.geometry().bounds()});
}

Map.centerObject(geometry)

My end goal is to be able to process my SAR images that were download from GEE.

Additional information:
I am able to display those images from GEE, but can't run any analysis, and each image has 2 bands as you can see here:

Edit2 - I thought maybe it doesn't work because the images are corrected already. I have tried to create stack but got the next error message:


texts:
I'm splitting a national soils raster file into regional (county) rasters, and organizing the soils classifications into a database keyed by region, eventually plot- or farm-level detail.

At the moment I'm doing it the hard way: using QGIS 3.4.4 to clip big rasters to smaller polygons.

As my result will be hundreds of smaller rasters with maybe 40 - 100 soil classes each (varies by polygon), I'll be 1) making lots of errors and 2) losing precision.

Occurs to me I don't need to split the raster at all, I just need to query the classes by each polygon, and write that to my database. Something like:

FOR each polygon IN vector.layer:
   CREATE.MASK FROM polygon1
   SELECT COUNT(raster.classes) as soil.area
   WRITE_DB polygon1.name, soil.area
   NEXT polygon

Maybe PostGIS is the better tool? Python? Dare I dig into R's formidable spatial libraries?

EDITS

My original framing was vague as I didn't have a lot invested in one tool over the others. Since then I've migrated over to a PostGIS database, read up on spatial indexes, and replaced rasters with vector data.

Decision-making went like this: our datasets are big, but unlikely to get bigger than the 28 countries of Europe; if I understood correctly, R likes to load spatial data into memory - which could cause trouble in some of these big joins.

We are likely to add layers of data - perhaps water features within some countries, soil moisture data in others - and I was impressed by the speed of st_transform() function in postgis when adding layers. Hope this is a more useful thread for those researching the similar questions.


texts: String format for geodatabases in ArcPy scripts

texts:
Several ArcPy functions take strings as geodatabase identifiers.  I've looked through the ESRI docs, but can't find a good reference for the various acceptable formats.  I'm looking for examples of each of the different kinds of geodatabases (personal geodatabase, file geodatabase, remote geodatabase).

For example,

gdb = "C:/MyData/PersonalGeodatabase.mdb"    # personal geodatabase
gdb = "C:/MyData/FileGeodatabase.gdb"        # file geodatabase
gdb = ""                                     # named instance SQL Server Express
gdb = ""                                     # sde database in SQL Server

If possible I would prefer that there not be a corresponding entry in ArcCatalog to make this work.  The current working examples I have seen here all seem to rely on specifically named connections in ArcCatalog.


texts: Computing driving directions in pgRouting/PostGIS?

texts:
I have calculated shortest route/path by using pg_Dijkstra algorithm in PostGIS and it is working fine. Now how to calculate driving directions from starting edge to end edge. I have gone through the link. But unable to go forward.

Is there any algorithm or procedure to calculate driving directions in PostGIS?


texts: Automating procedure to calculate least-cost path to destination locations from multiple source locations using ArcGIS Desktop?

texts:
I have 2 sets of locations represented by two point shapefiles in ArcGIS 10.1. My aim is to calculate the least cost path from each source location (let's call it A) to the destination locations (let's call them B) implementing the Tobler's hiking function and using ArcGIS path distance and cost path tools, as described for instance in this site:

http://www.mapaspects.org/node/3744

When it comes to calculating the least-cost path between just one A and all the Bs, I have no problem. But, since the procedure has to be repeated n times, being n equal to the number of As location, I am badly wondering if there is any way in ArcGIS to automate the procedure.

I tried to follow @Sam3000's very useful suggestion. I built a main model containing the Path Distance and the Cost Path tools (see attached image), and two submodels (nested into the main model) each iterating feature selection on the origin and destination locations respectively. But arcgis returns "ERROR 000581: Invalid parameters. Failed to execute (Path Distance).. Now, from what I can hypothesize looking at what happen when I run the whole model, is that arcGIS is first running the first submodel (Submodel_origin; bottom-left in the attached image) on its own, that is without passing its output (row-by-row) to the Path Distance tool (which is part of the main model). Maybe, there is setting that I am missing as to the way in which the operations performed by the different chunks must related to each other.

What this model should do in principle, as I have conceived it, is:

selecting an origin location at time (submodel_origin);
passing it to the Path Distance tools (which also takes as inputs a DEM and a SLOPE raster);
step 2 should provide  two rasters (namely, an anisotropic raster and a backlink raster) that with be feed into the Cost Path tool;
along with the 2 rasters in the previous step, another submodel (Submodel_destination) should select one destination location at time to be feed into step (3);
the output of the Cost Path is then converted into a polyline feature.


texts: Differences between Raster Dataset, Raster Catalog, and Raster Mosaic

texts:
I have about 40 scanned, colorful Geology maps that are GeoTiffs.  I want to bring them into an Esri FileGeodatabase to lower their size, combine them as one and display them with a mostly transparent overlay that is also scanned.

Right now they are individual Tif Images and I'm not sure which is the best way to store them in a FileGeodatabase.  Does anyone have experience with these and could you tell me when I would want to use one and not the others?

Eventually we will have an ArcServer, but for now they will be served from a Geosoft DAP server.  They will mostly be displayed in Arc 10.0, or 10.1, but will any of these options display well in Geosoft Target?

Edit: So after some more research I've found that none of these options are useable by Geosoft.  At least not in a File Geodatabase.  So what we will end up doing is making two mosaic rasters and save them as GeoTiffs with Packbits compression (which Geosoft supports, LZW on the other hand is not supported by Geosoft).  One for the Geology and one for the Overlays.  The compression will work great for the overlays which are mostly just one color with some black features, so the file size will be a fraction of the uncompressed size.

As to the difference between the 3 options I was considering, I'm still fuzzy on it.  So if anyone has a good explanation, I would still love to know...


texts: Populate field with result from a function?

texts:
I want to make a new column in every shapefile using for loop and populate that column using a function.

#the function

def find_folder(fn, k=-2):
    return os.path.split(fn)[0].split(os.sep)[k]

I tried:

fi = path
d = defaultdict(list)
for root, dirs, files in os.walk(fi):
    for shapefile in files:
        if shapefile.endswith('.shp'):
            arcpy.AddField_management(os.path.join(root,shapefile), "field_Tmp", "LONG")
            #here fill new field with 'find_folder(shapefile)' <---- How to do this?
           d[os.path.basename(shapefile)].append(os.path.join(root,shapefile))

Useful notes:

When the function is used, ex.

find_folder(shapefile)

it will give always the part of path that has a code. 1873, 8383, etc

Therefore, the field_Tmp will be populated with a number like that in every iteration for every shapefile found in the loop.

For people who know pandas it would be like this:

df = geopandas.read_file(shpfile).assign(field_Tmp = find_folder(shapefile))

Edit

Finally I am using the dictionary for the merge of all the shapefiles.

for cat, lst in d.items():
        if cat.endswith('.shp'):
            pd.concat(lst,sort = True).to_file(folder/cat)


texts: Alpha Shapes in pgRouting - missing FROM-clause

texts:
I am attempting to follow Anita Graser's post on creating alpha shapes in pgRouting.

The article was written several years ago, so some of the functions have been updated and names have changed. I have updated the code and it all appears to work fine apart from the last step. I am getting the following error for the last three lines of code:

"ERROR:  missing FROM-clause entry for table "home_catchment10km""

    CREATE TABLE "os_m2".home_catchment10km AS
SELECT *
   FROM "os_m2"."m2_roads_vertices_pgr"
   JOIN
   (SELECT * FROM pgr_drivingdistance('
      SELECT gid AS id,
          source,
          target,
          cost_time AS cost
      FROM "os_m2"."m2_roads"',
      674062,
      600,
      false,
      false)) AS route
   ON
   "m2_roads_vertices_pgr".id = route.id1

   CREATE TABLE "os_m2".home_isodist (id serial, max_cost double precision);
    SELECT addgeometrycolumn('os_m2','home_isodist','the_geom2',27700,'POLYGON',2);

INSERT INTO "os_m2".home_isodist (max_cost, the_geom2) (
SELECT 600, ST_SetSRID("home_catchment10km"."the_geom",27700) <**ERROR HAPPENS HERE**
FROM pgr_alphashape('SELECT id, ST_X(home_catchment10km.the_geom) AS x, ST_Y(home_catchment10km.the_geom) AS y FROM "os_m2".home_catchment10km where cost < 600'));


texts: Adding spatial reference to raster data using ArcGIS Pro?

texts:
New to ArcGIS Pro - I want to add a spatial reference to a raster file I have downloaded (in this case from Nearmap.com).

When I click properties of the raster file through the new project TAB in ArcGIS Pro I can not edit/change the spatial reference.

Is ArcCatalog 10 still the only way to do this?


texts: How to join all of the polygons that intersects and leave others without change in postgis?

texts:
What I need is to unite intersecting polygons in a single geometry and leave others without change.

I am able to do this using the following query, which unites all of the polygons into a single geometry and then dumps it to the separate ones.

INSERT INTO unioned(geom)
SELECT ST_MULTI((ST_Dump(ST_Union(geom))).geom) from footprints;

The problem is that I have more than a million rows in a database and ST_Union operation never ends. Is there any way to achieve the same without joining all of the geometries?


texts: How to add legend heading to a GeoServer WMS layer in QGIS?

texts:
I am trying to publish some layers in QGIS (version 2.14.0 Essen) to create WMS in geoserver. In ArcMap I can define legend heading

. I have not encounter any option in QGIS.  Since I am learning QGIS and Geoserver, wanted to make sure if I am missing any option (or any plugin can help)?


texts: How to convert NAD 27 to WebMercator in using the ESRI Silverlight API? (C#)

texts:
I am working with the ESRI Silverlight API. I found THIS post in the ESRI forums and THIS OTHER which show how to convert between WGS84 geographic coordinates and WebMercator coordinates. I want to do the same thing but from NAD27 Geographic coordinates to WebMercator. I found the spatial reference for NAD27 at the spatialreference website with the WKID 4267 and changed the code accordingly with the radius of 6378206.4 and degree 0.01745329251994328 respectively. I ended up with this function (in c#):

private MapPoint GeographicToWebMercator(MapPoint point)
{
        double a = point.Y * 0.01745329251994328;
        double y = (6378206.4/2) * Math.Log((1.0 + Math.Sin(a)) / (1.0 - Math.Sin(a)));

        double x = 6378206.4 * (point.X * 0.01745329251994328);
        return new MapPoint(x, y);
}

Unfortunately, my dots are displaced from where they should be. So, this is not working as intended. I was wondering if someone can take a look at it and tell me what is wrong... or if I can do this projection this way to begin with. The function, as is, is not checking for null "point" or any other data just yet. But, please, don't focus on that. I am just interested in the calculation per se.
Thanks!


texts:
I have a line-layer and want to style the line by using a data-table with no geography - I am using GeoPackage.

With other words, how can i get a field-value (i have two fields for filtering the unique entry) from the 2nd table in the field calculator to set the stroke / line width?


texts: QGIS: How do you overlap shapefile and raster with different coordinate systems?

texts:
I have a shapefile created from a GPS log and containing long/lat data in the WGS 1984 format.  I would like to overlay it on a raster (.tif) that has the following metadata:

Spatial data:

+proj=utm +zone=33 +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +units=m +no_defs

Layer extent:

5516499.9999999925494194,4837799.9999999403953552 : 5542600.0000000819563866,4860000.0000000167638063

Google satellite aligns perfectly with the shapefile, but the shapefile and Google satellite are about 354 meters due west of the raster.


texts: Parameters not valid in batch Alter field Arcpy

texts:
I am trying to perform a batch of Altering field names through a stand alone python script using Arcpy but it seems I am making a basic mistake in my python code. Using a list I refer to the original field names and with another list I refer to the new field names. When I alter the fields through arcpy.AlterField_management I refer to elements in both lists.

When I run the script, I get the error that an 'Invalid value type for parameter field' is used.

What mistake am I making here?

This is the script:

import arcpy
orig_fieldName_List = ['C1S1_N_VL_SFS_cwIS', 'C1S1_R_VL_HPR_cwIS', 'C1S1_N_VL_APT_cwIS']
new_fieldName_List = ['C_N_SFS', 'C_R_HPR', 'C_N_APT']

fc = r'C:\Users\Jelle\My Documents\Test.gdb\test'
fieldList = arcpy.ListFields(fc)
for field in fieldList:
    if field.name == orig_fieldName_List[0]:
        arcpy.AlterField_management(fc, field, new_fieldName_List[0], new_fieldName_List[0])
    if field.name == orig_fieldName_List[1]:
        arcpy.AlterField_management(fc, field, new_fieldName_List[1], new_fieldName_List[1])
    if field.name == orig_fieldName_List[2]:
        arcpy.AlterField_management(fc, field, new_fieldName_List[2], new_fieldName_List[2])


texts: differentiating stone houses from barren rock in ArcGIS 9.3

texts:
I'm working in ArcGIS 9.3 trying to differentiate stone houses (often red brick roofs) from barren rock in a typical Mediterranean karst region via remote sensing. I have been using Landsat imagery. I have been trying out various band combinations, multiplying and filtering bands in unsupervised classification. So far no good.

Does anyone have some helpful suggestion?


texts: ArcGIS 10.3.x - Publishing GP tools as GP services using Python 32- or 64-bit dependencies

texts:
Desktop environment:

Windows 10 64-bit,
ArcMap 10.3.1,
Python 2.7 32-bit (installed by ArcMap)

Server environment:

Windows Server 2012 64-bit,
ArcMap 10.3.1,
ArcGIS for Server 10.3.1,
Python 2.7 32-bit (installed by ArcMap),
Python 2.7 64-bit (installed by ArcServer)

My custom GP tool script is calling a 3rd party set of open source python libraries that have been developed to work within Python 2.7.x 32-bit framework. The GP tool works as expected on my desktop PC. I would like to publish my custom GP tool as a service on my remote server with both ArcGIS for Server 10.3.1 and ArcMap 10.3.1 installed on it. I am well aware that ArcGIS for Server demands a 64-bit framework for Python, so my question is: what would be the best way to publish my custom script tool (working on my desktop PC) as a service to be served within a web application (i.e. GP widget) without encountering conflicts and errors?


texts: PyQGIS script returning blank layer when applying singleband pseudocolor

texts:
So I'm trying to automate applying singleband pseudocolor style using a python script on a height map (.tif file).

Here's the code :

rlayer = QgsRasterLayer(inp, "pseudocolor_layer")
fcn = QgsColorRampShader()
fcn.setColorRampType(QgsColorRampShader.Interpolated)
lst = [QgsColorRampShader.ColorRampItem(0, QColor(0,255,0)), QgsColorRampShader.ColorRampItem(255, QColor(255,255,0)) ]
fcn.setColorRampItemList(lst)

# to make 5 classes
fcn.classifyColorRamp(5)
shader = QgsRasterShader()
shader.setRasterShaderFunction(fcn)
renderer = QgsSingleBandPseudoColorRenderer(rlayer.dataProvider(), 1, shader)
rlayer.setRenderer(renderer)
rlayer.triggerRepaint()

But I'm getting a blank layer. I ran these commands separately in the QGIS python console as well with no luck.


texts:
See photos. Both are screenshots of the same geotiff (orthophoto) image. Standard Windows (and other viewers) display the image correctly. QGIS shows the same image distorted.

The image is from Agisoft Photoscan. Orthophoto and QGIS Projection is "WGS 84 (EPSG::4326)".

When measuring an rectangual object in the orthophoto I got 105 degrees instead 90 degrees in one corner.

What's wrong/missing?


texts: How can I show "loading" while doing WMSGetFeatureInfo?

texts:
I have the following control in my page :

new OpenLayers.Control.WMSGetFeatureInfo({
        url: xxxxxx,
        name: 'infoControl',
        queryVisible: true,
        eventListeners: {
           getfeatureinfo: function (e) {
.
.
.}})

I want to show a "loading" message while doing the ajax call for WMSGetFeatureInfo.
How do I call my "loading" function the moment I click on the feature?
the event listener in the code above happens once the data returns...


texts: Transforming GCS_OSGB_1936 to WGS_1984_UTM

texts:
I am trying to add a DEM data layer to my map in ArcGIS 10.3.

The options available to transform the coordinate system to are:  WGS_1984_1,  WGS_1984_2,  WGS_1984_3,  WGS_1984_4,  WGS_1984_8.

Which one should I select?


texts: Using ArcPy to zoom to Y,X (string) coordinate pair copy/pasted from Google Maps?

texts:
I frequently use google maps to find locations, and paste the coordinates into the "go to XY" toolbar in decimal degrees. It's a bit repetitive to paste the latitude into one box, and the longitude into another. Is there a quicker way of zooming to a location with a string in the "google maps" coordinate format such as:

-29.948153, 146.864276

I was thinking I could make a python tool, but I can only find code for "zoom to selected features" (http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00s300000003000000)

I'm using ArcGIS Desktop 10.0.


texts: Sentinel-1 noise rebuilding using XML auxiliary files

texts:
I want to use XML auxiliary files in Sentinel-1 files to remove noise of HV-pol images, I have extracted these noise GCP in noise XML files into a point shapefile using python( in below ). However,
I don’t know how to rebuild the noise raster files?

I just know linear interpolation means, and I have tried as enter link description here.

However, the S-1 image size is 8992 * 10376 , script running error shows, MemoryError . Why popup the error message?

I posted my script includes to set noise field to those points and from point to raster in the end.

def Create_Points( root_dir, noise_dir, noise_value, ROWS, COLS, num_rows ,num_cols, rows, c ):

    env.overwriteOutput = True
    outpath = root_dir + noise_dir + r'measurement'
    newfc = 'newpoint.shp'
    ## Create points shapefiles
    arcpy.CreateFeatureclass_management( outpath, newfc, "Point" )
    newfc = os.path.join( outpath , newfc )
    cursor = arcpy.da.InsertCursor( newfc , [ "SHAPE@" ] )
    for i in range( len( rows ) ):
        for j in range( len( c ) ):
            cursor.insertRow( [ arcpy.Point( float( rows[ i ] ), float( c[ j ] ) ) ] )
    del cursor

    ## Add a new field to store noise values.
    source = ogr.Open( newfc, update=True )
    layer = source.GetLayer()
    noise_field = ogr.FieldDefn( 'noise', ogr.OFTReal )
    layer.CreateField( noise_field )
    #To set value to noise field.
    fields = [ 'FID', 'noise' ]
    with arcpy.da.UpdateCursor( newfc, fields ) as cursor:
        for row in cursor:
            row [ 1 ] = float( noise_value[ row[ 0 ] ] )
            cursor.updateRow( row )
    # Close the Shapefile
    source = None

    ##Create noise tiff file based on point shapefile.
    # Reshape noise_value to a matrix of [ num_rows, num_cols ] size.
    noise_value = np.asarray( noise_value ).reshape( num_rows ,num_cols )
    position = np.chararray( ( num_rows , num_cols ), 10 )
    for i in range( num_rows ):
        for j in range( num_cols ):
            position[ i ] [ j ] = str( rows[ i ] ) +',' + str( c[ j ] )
    position = np.asarray( position ).ravel()
    location = np.zeros( ( num_rows * num_cols, 2 ) )
    for id in range( len( position ) ):
        location[ id ] [ 0 ] = int( ( position[ id ] ).split(",")[ 0 ] )
        location[ id ] [ 1 ] = int( ( position[ id ] ).split(",")[ 1 ] )
    #Convert noise_value to one dim array.
    rows_min, rows_max, c_min, c_max = [ min( rows ), max( rows ), min( c ), max( c ) ]
    n_rows = ( int( rows_max ) - int( rows_min ) + 1 )
    n_cols = ( int( c_max ) - int( c_min ) + 1 )
    xi = np.linspace( int( rows_min ), int( rows_max ), n_rows/20 )
    yi = np.linspace( int( c_min ), int( c_max ), n_cols/20 )

    grid_x, grid_y = np.meshgrid( xi, yi )
    noise_value = noise_value.reshape( num_rows * num_cols )
    noise_inter = griddata( location, noise_value, ( grid_x, grid_y ), method = 'linear' )
    noise_inter = np.kron( noise_inter, np.ones(( 20, 20 )))
    ##Create the output image.
    driver = gdal.GetDriverByName( 'GTiff' )
    outputRaster = os.path.join( outpath, 'noise.tif' )
    dst_ds = driver.Create( outputRaster, COLS, ROWS, 1, gdal.GDT_Int16 )
    band = dst_ds.GetRasterBand( 1 )
    band.WriteArray( noise_inter )


texts: Efficient way to multiply specific raster with specific variable python

texts:
I have 4 separate raster bands of the same image and 4 different coefficients (one for each band) stored as variables. I would like to iterate over the raster folder, and multiply each band with its corresponded coefficient. I can recognize each band by its file name. Let's say I have in folder c:\bands these bands:

sn2_b1.tif
sn2_b2.tif
sn2_b3.tif
sn2_b3.tif

And let's say these are my coefficients:

coeff_b1 = 10
coeff_b2 = 20
Coeff_b3 = 30
Coeff_b4 = 40

My goal is iterate over the bands, multiply each one with its corresponded coefficient and save the result as a new raster. My solution so far is quite cumbersome:

import gdal
import os
src = r"c:\bands"
dst = r"c:\cal_bands"
for i in os.listdir(src):
    if i.endswith(".tif"):
        orig_tif = os.path.join(src,i)
        band_ds = gdal.Open(os.path.join(src,i))
        band_DN = band_ds.GetRasterBand(1).ReadAsArray()
        geotransform = band_ds.GetGeoTransform()
        geoproj = band_ds.GetProjection()
        xsize = band_ds.RasterXSize
        ysize = band_ds.RasterYSize
        band_DN = band_DN.astype('f')
        format1 = "GTiff"
        driver = gdal.GetDriverByName(format1)
        dst_datatype = gdal.GDT_Float32
        if i.split("_")[1] == 'b1.tif':
            cal_band =  band_DN  * coeff_b1
        elif i.split("_")[1] == 'b2.tif':
            cal_band =  band_DN  * coeff_b2
        elif i.split("_")[1] == 'b3.tif':
            cal_band =  band_DN  * coeff_b3
        elif i.split("_")[1] == 'b4.tif':
            cal_band =  band_DN  * coeff_b4
        new_name = "cal_" + i
        output_path = os.path.join(dst,new_name)
        dst_ds = driver.Create(output_path ,xsize,ysize,1,dst_datatype)
        dst_ds.GetRasterBand(1).WriteArray(cal_band)
        dst_ds.SetGeoTransform(geotransform)
        dst_ds.SetProjection(geoproj)
        dst_ds = None

This script works great but is there a way to make it less cumbersome? In this example I have only 4 bands, which is easy to write all the ifs, but I might have 20 bands or more.


texts: Drawing series of lines from hub in QGIS?

texts:
I'm trying to create a route map and was wondering how best to do this. I have a series of 'routes' emanating from hubs to single destinations and wish to link them creating a spoke-like effect like an airline route map.

Is there a way in QGIS 3 to do this?


texts:
I want to upload waypoints recorded in the Maps.me App onto openstreetmap.org for further editing.

Already tried:

Export from Maps.me to .kml then conversion to .gpx with QGIS (second
try with geotracker app)
Then I added a fake trackpoint (because I had only recorded waypoints in maps.me) according to https://wiki.openstreetmap.org/wiki/Upload_Waypoints.
Then the .gpx upload to OSM seemed to work but only the fake trackpoint was loaded.

What else can I do?


texts: Cesium, how to persist baselayer selected by user?

texts:
For Cesium, I would like to know how to persist the baselayer selected by user?

For example, if a user selected "Mapbox Streets", how can we detect it was selected and persist it, so that the same "Mapbox Streets" will be used when it is reloaded.

Is there a callback method, so that we know which baselayer the user selected?


texts: I'm using Terrasync on a Trimble GeoXH 6000 - set for WGS 84 and UTM. It acquires satellites fine, but instead of displaying my current location in UTMS, the skyplot screen shows "12RVU" and then a single 10 digit number. For the life of me I can't figure out what 12RVU is, or how to make it display actual coordinates. I've been using Terrasync for years and have never had a problem like this before.


texts:
I am trying to get the network graph of a specific reason of USA. I have got a couple of options. One of them is using osm2po tool to extract the nodes and edge table from the osm map of a specific reason.

I tried with this one and I could get the graph in the database table. But I think that there is property called geom_way, that defines the intermediate nodes of each edge. I mean suppose I have an edge like this

edgeid = 1;
source_node = 1;
target_node = 152
source_lat_lon = -93.27368,44.98204
dest_lat_lon = -93.27254,44.98124
geom_way = 0105000020E6100000010000000102000000040000006F641

Then if I decode the geom_way it will be the coordinates of the intermediate nodes between source_lat_lon and dest_lat_lon.

Well I need to run the shortest path algorithm on this graph. So I was wondering if it could create some issues.

Also I was wondering if there is any alternative tool or way to get the graph of a specific area of USA except osm2po.


texts: Applying 3D-Helmert Transformation Parameters in QGIS?

texts:
I have received a set of local transformation parameters from my surveyor to go from WGS84/UTM36S to ARC1950/UTM36S (See picture below)

The way QGIS does the transformation is not the best. I am working in Malawi and the default transformation parameters give misplacements of 20 to 40 m. This is why I ask a very reliable local surveyor to provide the local parameters for the area of interest.

What is the proper way to write a proj line to introduce these parameters into QGIS?

I was trying: +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +x=249.7578 +y=2.8706 +z=203.0842 +rx=-14.90782 +ry=-6.25396 +rz=0.20011 +s=-7.3003

But it does not seem to work.


texts: Adding GPS coordinate points to point feature class using ArcGIS Pro

texts:
I have a raster, and in ArcGIS Pro 2.4, I created a point shapefile with 830 points. I forgot to specify that I wanted to show the GPS points in the attributes table of each point. However, I need to add the GPS coordinates as another field.

How can I do this without creating another point layer and creating all those points again?


texts:
I have been attempting to find a solution to the problem of limited access to GRIB_API(on Windows) through using GDAL in Python.

The main thing I want to achieve is to be able to set some conditions, then iterate through folder(s) for a given time period and produce a raster that would represent the amount of times the conditions were met in a format that is easily handeled by most GIS software, such as GeoTIFF.

I have written up two blocks of code either of wich should do this for me, but I at best get a blank raster. I assume it is either a flush to disk problem or I am missing something to translate integers into GeoTIFF.

  This is one version:

import gdal
import numpy

src_ds = gdal.Open('inca_20140101-0000.grb')
format = "GTiff"
driver = gdal.GetDriverByName(format)
dst_ds = driver.CreateCopy('01012014.tif', src_ds, 0)

def krog(baza,novi):
  x = 0
  while x < 25:
      baza = gdal.Open('01012014.tif')
      if x < 10:
          novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
          if novi > 1 and novi < 5:
              baza = baza + 1
          else:
              baza = baza + 0

      elif x > 10:
          novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
          if novi > 1 and novi < 5:
              baza = baza + 1
          else:
              baza = baza + 0
      x += 1
dst_ds = None
src_ds = None

  This is the other:

import gdal
import numpy

dst_ds = driver.Create('01012014.tif', 400, 300, 1, gdal.GDT_CInt32)
dst_ds.SetGeoTransform([4468800, 30, 0, 1223500, 0, -30])

srs = osr.SpatialReference()
srs.SetUTM(11,1)
srs.SetWellKnownGeogCS('LCC')
dst_ds.SetProjection(srs.ExportToWkt())

raster = numpy.zeros((400, 300), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)

def krog(baza,novi):
  x = 0
  while x < 25:
      baza = gdal.Open('01012014.tif')
      if x < 10:
          novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
          if novi > 15 and novi < 25:
              baza = baza + 1
          else:
              baza = baza + 0

      elif x > 10:
          novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
          if novi > 15 and novi < 25:
              baza = baza + 1
          else:
              baza = baza + 0
      x += 1
dst_ds = None

GRIB files have 4 bands (Temperature, Humidity, u and v values) and are provided for every full hour in a day (most of the time).

EDIT: Following advice I adopted the second version like this:

import gdal
import numpy
import osr

ime_rez = str(input("Poimenuj datoteko z rezultati:") + ".tif")
#User can provide a name for the output picture
format = "GTiff"
driver = gdal.GetDriverByName(format)
dst_ds = driver.Create(ime_rez, 400, 300, 1, gdal.GDT_CInt32)
dst_ds.SetGeoTransform([4468800, 30, 0, 1223500, 0, -30])

srs = osr.SpatialReference()
srs.SetUTM(11,1)
srs.SetWellKnownGeogCS('LCC')
dst_ds.SetProjection(srs.ExportToWkt())

raster = numpy.zeros((400, 300), dtype=numpy.uint8)

def krog(baza,novi):
    x = 0
    while x < 25:
        baza = gdal.Open(ime_rez)
        if x < 10:
            novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
            #This should open the next file in line
            if novi > 15 and novi < 25:
                baza = baza + 1
            else:
                baza = baza + 0

        elif x > 10:
            novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
            if novi > 15 and novi < 25:
                baza = baza + 1
            else:
                baza = baza + 0
        x += 1
dst_ds.GetRasterBand(1).WriteArray(raster)
dst_ds = None

This now results in a "ValueError: array larger than output file, or offset off edge". I believe this is a problem because the cell size in .grb files is 1000 m, but I have yet to find a way to manually set it in the newly created .tif.


texts: How to store the CHANGES to a map database

texts:
I'm looking into providing a map amendment service for a particular map series.
The service will help map owners look into the changes made to the operational features of the map, after the map was published.

There are three basic types of changes:
1) An object has been added.
2) An object has been deleted.
3) An object has changed, either an attribute or its geometry.

For object changes, the user should get both the old and the new value of the attribute or geometry.

Now, I'm looking into a clever way to implement this.

The basic idea is that the database administrator edits the map data in an interface, just as he would edit a normal map, and the rest is taken care of behind the scenes.

As I see it, there are two basic approaches:
1) Store the changes in a separate table, and keep the main database current.
2) Timestamp the records in the main database.  In case of changes, I cannot change a record, but timestamp the old one as expired and add a new urrent one.

My main concern is not performance (not huge amounts of data), but ease of implementation.  Is there a way I can make the DBMS (say, PostGIS) handle this automatically (for example by using triggers), Or can I build my own interface on top af an exisiting GIS (say, QGIS or ArcGIS) to handle the storage of records in the various tables as the database administrator makes changes to the map. Preferrably, I'd like a Web-interface.

Any Ideas anyone?


texts: Random assignment of attribute table values to features

texts:
I have the current SQL expression, which I have already had help with:
SELECT j_pipeline_midpts.*, j_pipeline_inc.*
FROM   j_pipeline_midpts
INNER JOIN
       j_pipeline_inc
   ON  j_pipeline_midpts.LSD = j_pipeline_inc.LOCATION
  AND  j_pipeline_midpts.Licence_No=j_pipeline_inc.LICENCE_NO

I am new to SQL expressions and QGIS tools and would like to get some help on how I could possibly solve the issue detailed in the attached picture.
However, for cases like I have pictured, the above expression does not work. There are multiple features with the same field key values.
All lines in the grid cell are attached to a grid cell number (LSD), and each line has a licence number (LICENCE_NO), but sometimes, most times, the licence number is the same for each line in that grid cell – but each line has a unique OBJECTID. I have placed midpoints along each line with the associated line attribute data, which I have called j_pipeline_midpoints.
The problem is, as expected, that QGIS attaches all associated attributes from the join layer (j_pipeline_inc [red diamond]) with the same licence number to each associated point target layer (j_pipeline_midpoints [green triangles]) with the same licence number; unfortunately, though, each row represents a different incident in time. Attaching all 9341 licence number entries from the join layer to all associated target layer midpoints with the licence number 9341 doesn’t make sense.
I would like to randomly place only one entry from the join table (pipeline_inc [red diamond]) to a separate point on the line [green triangle] with the same licence number.  It doesn’t matter what points on the line in the grid cell that the entry joins to, if the number of features is > 1, only that it is associated with an equal number of points as entries within the join table.


texts: Use of OpenStreetMap France in another country

texts:
I've just stumbled across OpenStreetMap France and was amazed to find that it has 2 extra zoom levels, 19 and 20, whereas OpenStreetMap only goes up to 18.

In our application we display potentially masses of geographic data over the top of basemaps. We use Ordnance Survey and OpenStreetMap. The Ordnance Survey mapping we use we can render off at higher zoom levels, which will (once we've upgraded our GeoServer) allow our application to zoom closer and see greater details in the overlaid data. But by default we use OpenStreetMap as it is quicker, looks better (especially far zoomed out) and we like trying to get people to use OSM. Sadly this means we can only go to zoom level 18 (unless I do the OpenLayers thing to 'force' it to zoom in which drops the resolution).

However, if we were to use OSM France we could get those extra zoom levels, and as far as I can see it is (almost) as quick and doesn't appear to have any language diffrences (apart from London being displayed as Londres, which I am not concerned about as we are based in Dorset).

I guess my question is twofold.

1 - Can anyone think of a good reason for me not to use OSM France

2 - Why doesn't standard OSM have these two extra zooms (capacity?) and will they ever?


texts:
I have a table with 2 fields, field1 contains site specific data, field2 contains a number which groups the field1:

Field1, Field2
180, 1
203, 1
220, 1
205, 1
160, 2
140, 2
135, 2

I want to sort low to high of field1 grouped by field2.

I then want to take the "rank# of field1 grouped by field 2" SUBTRACT the "median of field1 grouped by field2" and divide the answer by "the total # of points of the field2 group"

I'm looking for some direction, a few ideas or a formula I can use in QGIS field calculator to help solve this problem. Or some idea's on a python script for the advanced python calculator.


texts: Creating and Joining Points in Polygons Based on Counts

texts:
I have a polygon shapefile of aggregated commuting journeys that has >10,000 features in it. Each feature has a unique ID. I have a csv file that contains thousands of rows. Each row represents one feature from the shapefile. Each polygon can be represented by numerous rows, where each row is the count of journeys from the original polygon to the destination polygon. I have attached a picture that may better explain what I'm trying to articulate.

I want to create points in both the origin and destination polygons (joined by lines) to represent individual journeys. Does anybody know how using QGIS/Python etc. I could accomplish this?


texts:
I am doing a map which refresh dinamically all markers when they get modified.
I want to add popups to show some info.

I create popups this way:

        var feature = new OpenLayers.Feature(markers, location); //location = lonlat
        feature.closeBox = true;
        feature.popupClass = OpenLayers.Class(OpenLayers.Popup.FramedCloud, {
        'autoSize': true
    });;
        feature.data.popupContentHTML = "info";
        feature.data.overflow = (false) ? "auto" : "hidden";
            var markerClick = function (evt) {
            if (this.popup == null) {
                this.popup = this.createPopup(this.closeBox);
                map.addPopup(this.popup);
                this.popup.show();
            } else {
                this.popup.toggle();
            }
            currentPopup = this.popup;
            OpenLayers.Event.stop(evt);
        };
        mark.events.register("mousedown", feature, markerClick);

Now, I have the problem about what to do if one of my markers update its info. How can I destroy the popup? (I will create other one) or how can I modifiy it? Can I access them using lonlat?

Sometimes, I want to delete all of them, is there a way to do that? like deleting a layer or something?


texts: How to georeference PNG images using QGIS?

texts:
How can I georeference png images using QGIS, running on OSX?

I've looked at ExifTool, and the plugins GeoTag and Import Photos + Photo2Shape. These tools (and other answers) as far as I can tell all use georeferenced images or images with EXIF data. I would like to load the png somewhere, and drag it around on a basemap from OpenStreeMap to position it where it should be. I know the general location, and just need to position it.

(I think I know the projection, but will find out if that is the case. The images are of an area ~500m x 200m so the projection variation will hopefully not matter hugely if I am incorrect).


texts: Stop ArcGIS model until all values are collected and then merge?

texts:
This is similar to Using If precondition is false, continue iterating next input in ModelBuilder? but I want the iterate tool to run to collect all the values and only then continue to the merge. Any ideas? Currently it just runs the whole process for each file (which works but takes much longer than if the merge is done only once) or I have to manually break the collect-merge link, run the iterate, reconnect the collect-merge link and then continue running the model which is a pain as it is a lot of towers that need processing and this collect values-finish getting all values and then continue will be required in other models as well.

I tried the "Select Data" tool as well but this doesn't allow for wildcard searches like the Iterate does -my Iterate has the "_Graded*" as the search param as the gdb included many other files related to each tower chunk.

My STOP is set to "True" (tried false as well). I also tried to connect the Stop to the Iterate instead of Collect values but same result.

=== SCRIPT ADDED

I just created a python tool for this and it collects all the values but it passes them back to the model as a long list -how can I get it to output in a proper format?

    import sys, os, arcpy

GDB=arcpy.GetParameterAsText(0)
Search=arcpy.GetParameterAsText(1)
Type=arcpy.GetParameterAsText(2)

##GDB=r'P:\2012\044_DM_Radio_Propogation\Working\Propogation\Propogation_28355.gdb'
##Search='_m'
##Type='Polygon'

path=os.path.dirname(GDB)
arcpy.env.workspace = GDB
List=[]
for item in arcpy.ListFeatureClasses('*'+Search+'*',Type):
    List.append(path+"\\"+item)
    #arcpy.SetParameterAsText(3, item)

arcpy.SetParameterAsText(3, List)
arcpy.AddMessage("List: "+str(List))


texts:
I am trying to follow the tutorial posted here.  When I get to this part "Extraction of seasonal mean temperature" and run this code:

t.rast.extract \
  input=temperature_mean_1950_2011_seasonal \
  output=temperature_mean_1950_2011_spring \
  where="start_time = date_trunc('year', start_time \
         + interval '2 months')" \
expression="temperature_mean_1950_2011_seasonal/100.0" \
  base=temp_mean_spring

I get the following error:

ERROR: Unable to get map ids from register table
   <temperature_mean_1950_2011_seasonal_PERMANENT_raster_register>
Traceback (most recent call last):
File "/usr/local/grass-7.0.svn/scripts/t.rast.extract", line 91, in <module>
main()
File "/usr/local/grass-7.0.svn/scripts/t.rast.extract", line 85, in main
base, nprocs, register_null)
File "/usr/local/grass-7.0.svn/etc/python/grass/temporal/extract.py", line 95, in extract_dataset
rows = sp.get_registered_maps("id", where, "start_time", dbif)
File "/usr/local/grass-7.0.svn/etc/python/grass/temporal/abstract_space_time_dataset.py", line 905, in get_registered_maps
dbif.cursor.execute(sql)
sqlite3.OperationalError: near "'2 months'": syntax error

I am using the data specified in the presentation (mean temperatures dataset downloaded from http://eca.knmi.nl/download/ensembles/data/Grid_0.25deg_reg/).

Please help.

Thank you,
Bogdan Rosca


texts: OSM Display a particular changeset-version of a way

texts:
Is it possible to visualize an particular changeset of an way?

It's about this way:

http://www.openstreetmap.org/way/24837452

At the moment it seems to be a part of an MultiPolygon, but in the past it wasn't. I have an OSM file which is a bit older and contains this version:

<way id="24837452" version="64" changeset="13731788" uid="362915" user="micha555" timestamp="2012-11-03T09:11:23Z">

Even when the first and last node in this version of the way is the same, mongodb tells me this polygon is invalid. So I would like to visualize it.

If I try to open the history

http://www.openstreetmap.org/way/24837452/history

I always get errors


texts: Enable/Remove popup on mousedown/mouseup?

texts:
I'd like to do as the title suggests and can get the popup to appear on mousedown with the following code, but mouseup doesn't close the popup:

// Quick look at info on mouse hold
marker.on("mousedown", function(e){
    this.openPopup();
});
// TODO: can we get mouseup working?
$(marker).mouseup( function(e){
    this.closePopup();
});

Am I missing something?  If I log within the functions I see the mousedown hitting but mouseup never happens.


texts: Replicate the Google Earth Engine satellite basemap for analysis in code editor

texts: I need high resolution cloud-free raster imagery (<1m resolution) for Montana, USA. Visually inspecting the Google Earth Engine satellite basemap looks to be higher resolution than NAIP imagery (which is the highest resolution I am currently aware of). But I understand that raster analyses can't be conducted on the Earth Engine basemap.  How do I replicate what I'm seeing on the Earth Engine satellite basemap with imagery from the archive to get the same level of resolution?


texts: Change MultiLineString direction in QGIS/PyQGIS

texts:
I am looking for a method in PyQGIS (QGIS3.x) of reversing the direction of a MultiLineString feature. If there is more than one part, the order of the parts should be reversed.

I have found a similar question -
How can I switch line direction in QGIS? - for LineString/Polyline.  - which only works for QGIS2 - and modified it slightly for use in QGIS3. However, I can't work out how to adapt the approach for MultiLineStrings.

layer = qgis.utils.iface.mapCanvas().currentLayer()
fid = 1 # feature id
feature = layer.getFeature(1)
geom = feature.geometry()
nodes = geom.asPolyline()
nodes.reverse()
newgeom = QgsGeometry.fromPolylineXY(nodes)
layer.startEditing()
layer.changeGeometry(feature.id(),newgeom)
layer.commitChanges()

In theory the above approach could be used with
    newgeom = QgsGeometry.fromMultiPolylineXY()
but I can't see how to generate the necessary QgsMultiPolylineXY object.

I have also noted the reversed() function on the QGSMultiLineCurve class, but I couldn't attach the generated QGSMultiLineCurve to the original feature.

I'm not really interested in UI-based approaches, as the features to be reversed are determined programatically.

Any approaches are welcome - if there's a better general approach than above, I'm happy to hear.


texts: ArcPy x64 - Create Enterprise Geodatabase in SQL Server Fails

texts:
I'm using ArcGIS Desktop Standard 10.1 SP1 on Windows 7 SP1 x64. SQL Server Native Client 11.0 is installed for both 32 bit and 64 bit. I've also got SQL Server 2012 x64 installed (client tools and database engine).

Running an ArcPy based script in Python 2.7 x64 crashes Python.exe when it tries to create an Enterprise Geodatabase in SQL Server 2012 with this line:

arcpy.CreateEnterpriseGeodatabase_management("SQLSERVER", sDbServerName, sDbName, "OPERATING_SYSTEM_AUTH", "", "", "DBO_SCHEMA", "", "", "", sLicenseKey)

The error is:

  Unable to load server library (c:\program files
  (x86)\arcgis\desktop10.1\bin\sdesqlsrvr101.dll), Error = -324

Using SDE_SCHEMA and/or database authentication don't work either.

The workaround is to use 32 bit Python. However, our standard development environment is x64 (for performance reasons).

Any ideas how to get this working in x64 Python?


texts: Getting elevation data or Z value of point using QGIS

texts:
I know that there is 3d drawing in QGIS. However, I want to learn that how it converts 2d map to 3d.

Is there a way to get height value of points or does it convert longtitude, latitude to xyz format?


texts: GeoServer Limit SLD editing functionality to certain users in GeoExplorer

texts:
How would I go about preventing user X from editing layer styles? I have created a user in GeoServer that if they login they only see the layer preview screen. However as I am using the Boundless Suite (OpenGeo) if they add a new site in the GeoExplorer,they login and add a new point/ However the same use can also edit the style for that layer. I have hidden all the other layers, however this layer must be visible.

Is there a way to either remove the style editor from GeoExplorer? or prevent users from editing styles from the GeoServer side?

I did try creating a edit feature map using the OpenGeo SDK however found that adding the user login functionality pretty complex with little easy to follow documentation. I had all the functionality I required, as I could turn off the edit style feature, but then anyone who managed to access the layer would be able to add a new point.


texts: QGIS 3.x Python Plugin: Fetch Path to QgsMarkerLineSymbolLayer Marker

texts:
I am new and currently working on an experimental plug-in for QGIS 3.x. I want to fetch a custom SVG marker used in a MarkerLine (QgsMarkerLineSymbolLayer as in https://qgis.org/api/classQgsMarkerLineSymbolLayer.html) within a QGIS 3.x project from the local file system and copy it to a ZIP archive. Locating the marker in a point layer as in QgsSvgMarkerSymbolLayer is trivial, because QgsSvgMarkerSymbolLayer supports a path() function (https://qgis.org/api/classQgsSvgMarkerSymbolLayer.html). I've searched for a similar function in QgsMarkerLineSymbolLayer, but couldn't find it.

I seek any hints or code examples.


texts: Retrieve all feature IDs displayed on the viewport into JSON using Leaflet

texts:
I have build a map application using leaflet, where I display a WMS layer from Geoserver. What I am trying to do now is to get all the elements which are currently on the Bounding Box and store them in a JSON.
I have read here, that leaflet applies a BB on the fly according to the position of the viewport. But what I don't know is how I can export the fid's of the elements into a JSON or array etc.
Ideally I want to be able onbutton click to get the JSON with all the features displayed in the viewport.


texts:
I would like to calculate approximately, how many tiles of given size would need to be rendered for a given bounding box and zoom level range in a given projection (always 3857).

For example, for:

Tilesize: 256x256 pixels
Longitude: 15-25°
Latitude 10-20°
Zoom levels 5-15

I know that number of tiles at zoom level nx = 4^x.

I think the procedure goes like this: add 1 tile for every level in which one tile is still enough to cover the entire bounding box. For every level further down, assume tiles scale based on the equation above.

Therefore my question is: How do I know the last zoom level on which I can cover the bounding box with 1 tile?


texts: Python Stand Alone Script Parsing error SyntaxError: invalid syntax (line 1)

texts:
I'm pretty new to python, but I'm trying to transpose a time field
 for a set of data I'm working with in ArcMap. I changed the file locations and fields to my needs. However, when I try to run the script in the python window I get the error

  Parsing error SyntaxError: invalid syntax (line 1)

It seems like this is an error with importing arcpy, but I believe ArcMap contains ArcPy. Additionally, running just import arcpy in the scripting window works fine.

Here is the Code:

# Name: TransposeFields_Ex_02.py
# Description: Tranpose field names from column headers to values in one column
# Requirements: None

# Import system modules
import arcpy
from arcpy import env

# set workspace
arcpy.env.workspace = "H:\Users\bdf006\Desktop\GIS Review\Pittsburgh Roads\RoadsWithTime.gdb"

# Set local variables
inTable = "Pittsburgh_Roads_Counts"
# Specify fields to transpose
fieldsToTranspose = "1a;2a;3a;4a;5a;6a;7a;8a;9a;10a;11a;12p;1p;2p;3p;4p;5p;6p;7p;8p;9p;10p;11p;12a"
# Set a variable to store output feature class or table
outTable = "Pittsburgh_Roads_Counts_With_Timw"
# Set a variable to store time field name
transposedFieldName = "Transposed_Time"
# Set a variable to store value field name
valueFieldName = "Value"
# Specify attribute fields to be included in the output
attrFields = "Shape;Type"

# Execute TransposeTimeFields
arcpy.TransposeFields_management(inTable, fieldsToTranspose, outTable, transposedFieldName, valueFieldName, attrFields)


texts: Filter FeatureCollection by the number of vertices

texts: I have a FeatureCollection of polygons, I was wondering how I filter (remove) polygons that has the number of vertices smaller than a scalar.


texts:
I am looking to host lidar data on my own server and give clients the ability to view that data through a web browser as a 3D point cloud (file would be in .las format or equivalent), possibly with the ability to measure features.

Does anyone know of an easy to implement piece of software that can do this?

I am not nessecarily looking for open source, and I would prefer something plug and play because I am not a programmer.


texts:
I'm currently working on reachabilities (isochrones) to certain locations in QGIS. So far I've tested the ORS plugin and the Hqgis plugin. Now I'm looking for a way to calculate the boundary of between two or more reachabilities.
E.g. I drive from place A to place B with a car. And a second car starts from place B in direction of A.
The calculation should stop at the point the cars meet. In the result I want to have to areas.
Does anyone have any experience with the calculation?


texts: How to create a 4096x4096 png from a 1440x720 data set using GDAL?

texts:
I'm having a bear of a time trying to do something that seems fairly straightforward. I'm trying to take a netcdf file from ftp://eclipse.ncdc.noaa.gov/pub/OI-daily-v2/NetCDF/2011/AVHRR-AMSR/ and turn it into a 4096x4096 png so that I can tile it up for use w/google maps.  I've been using the following commands but it always seems to create blank areas on the east and west borders:

/usr/local/bin/gdal_translate -a_srs EPSG:4326 -outsize 4096 4096 -a_ullr 0 90 360 -90 NETCDF:"input.nc":sst -of netCDF output.nc
/usr/local/bin/gdaldem color-relief output.nc -b 1 -alpha colorscale.txt -of PNG output.png
/usr/bin/convert output.png -transparent #FFFFFF output.png

Am I just approaching this the wrong way?  Am I using the wrong tools?

FWIW - gdal_translate is v1.8.0

Thanks in advance.

UPDATE: "Answer" below from Dija - colleague of mine - is actually a response to the comments on the question itself.

UPDATE  Accidentally added an answer instead of a comment.  Moving it here.

I thought I'd weigh in here. I've been working on this issue for a little while now. I've tried many iterations of the process. I scaled these images to the 4096x4096 as it actually creates a better image than if not. Meaning, if I leave the image in it's native size, 1440x720 I get a larger gap at the prime meridian (for this data set) and at the 180th meridian for other data sets. I scaled the output image to be a multiple of 256 (to fit better with google) and get various results.

If I scale it to 1536x768 the gap shows up to the right of the 0 meridian. If I scale to 1280x640 it's to the left. If I leave it alone, it's pretty much splits it.

I have taken the base image that gdaldem spits out, cut it in half and matched the edges up...they line up perfectly. This leads me to think that there's something weird going on with gdal2tiles.py

I wrote my own tiler that I'm getting ready to turn loose on our images to see if the gap remains.

I'll put together some screen shots to show what I am talking about with how the gap as well as report on what my tiler does.

Thanks for any input you might have. It's appreciated.

UPDATE 2 I ran the maptiler I wrote against the original raster and it works fine with out a gap.  gdal2tiles.py seems to be the problem.  It introduces a transparent edge to the last tile that touches the 180th meridian.


texts: Building "floor switcher" toggle button in Web AppBuilder for ArcGIS?

texts:
I've been trying to make a Floor Switcher for my web app in WAB. There is one MapService of my building which i could simple query (layer.setDefinitionExpression("exp")) based on attribute. In-panel widget isn't a good solution because i want to change floors without open any panel. Simple toggle buttons or node list directly in app view will be the best solution.

I found these two topics, one with toggle button, but i don't know is it possible to implement this in WAB app.
How do I use a toggle button with the ESRI javascript api
Second one is that what i'm looking for but unfortunately in panel widget, so isn't very user-friendly.
https://geonet.esri.com/thread/158023


texts: Segmentation Image used to classify RGB image in ERDAS and GIS

texts:
My ERDAS 2013 does not have the advanced Image Objective Extension for feature extraction..
I am using thus the tool Image Segmentation (Unsupervised -> Image Segmentation) to classify my RGB image ...
when I try to use the Segmented Image as a base for classification ( the signature combined of unsupervised and supervised, I also tryed a saparate unsupervised and supervised signature)  I constantly obtain a message: Layer count in input image (1) does not match layer count in
signatures (3)

should I vectorise the segmentated picture, obtain zonal statistics from the original image based on the segmanted polygons (ARCgis)? I heard many use this method, but how to classify the picture afterwards, since the result will be only in the table?

btw. I know SPRING, MONTEVERDI, ORFEO, SAGA, those tools are not able to work on my huge set of data, and I do not have time to split my pictures in tiny elements...but ERDAS was able to segment and classify this size of data....


texts: How to Simplify a Line Network Preserving Topology?

texts:
I have a Shapefile (consisting of European major roads) with about 250.000 Segments which I have to simplify for pgrouting. But I can't seem to find a way to do that properly.

This is what it looks like:

and this is what it should look like:

I somehow have to remove every Point of the Lines which is connected to less than 3 Lines (being not an intersection) while preserving the topological connections between the remaining points. If anybody has an idea, it would be greatly appreciated!

Best regards

EDIT: I tried to implement the idea of @dkastl and managed to get only the unneccessary nodes (nodes with only 2 adjacent linestrings) from my network with the code below (network generation taken from underdark's blog http://underdark.wordpress.com/2011/02/07/a-beginners-guide-to-pgrouting/):

SELECT * FROM
   (SELECT tmp.id as gid, node.the_geom FROM
     (SELECT id, count(*) FROM network
     JOIN node
     ON (start_id = id OR end_id = id) AND (end_id = id OR start_id = id)
     GROUP BY id ORDER BY id) as tmp
   JOIN node ON (tmp.id = node.id)
   WHERE tmp.count = 2) as unn_node;

So, all I now have to do is the merging of the lines. However, I have no clue how. I imagine it has to be a loop which for every row of the result of above query gets the adjacent lines and merges them. Then it would rebuild the network completely and repeat the process until the query above returns an empty result.


texts:
I'm looking for a way to orthogonalize building outlines which I've extracted from a high-res DEM. I'm using the FOSS GIS stack (QGIS, GRASS, SAGA, PostgrSQL, PostGIS) and can code in Python (or Java).

I want any building corners which meet at nearly 90 degrees to be snapped to exactly 90 degrees, by moving nodes.
This needs to work on multiple/selected features
I don't want to manually edit the polygons... there are way too many!

Similar functionality exists in OpenStreetMap, for example in the ID and JOSM editors. If you select a building in the ID editor and press 'O' (letter oh), the corners get tidied up.

This is a blessing for people like me who have poor hand/eye co-ordination when tracing from aerial imagery :) So I know that it's technically possible, and can always reverse-engineer those.

So far I've...

isolated buildings from the DEM using SAGA, using the DTM Slope-based filter
used Raster Calculator in QGIS to get a binary raster (1=building, 0=not building)
converted raster to vector using GDAL
simplified the polygons, using 0.5m tolerance, to get rid of the 'staircase' effect, where the polygon follows the pixel edges

Now, I have my rough building outlines, I just need to orthogonalise them.

Does anyone know of any tools / plugins / libraries in the FOSS GIS stack which provide this functionality already? Have I overlooked something obvious?


texts:
I have 2 layers, one which is a point layer with entertainment locations and the other is a polygon layer of buildings. I'm tasked to find how many buildings are within a 180m radius of the entertainment locations.

Obviously my first instinct was to create a buffer and then use a spatial join however, if I check the Join_Count column afterwards and compare it to the Select by Location result they don't match which means the Spatial Join counts something extra. After fiddling with the data I came to a conclusion that the Spatial Join is counting other near overlapping buffers.

I'm using ArcGIS 10.7

Edit: the question is how do I get to see the correct number of building inside each buffer (or how do i find out how many buildings i have around each entertainment spot in a 180m radius)

Photo for reference:

https://imgur.com/a/s2EunoL - you can see in this photo that there are 126 selected using the Select by Location option for Object 3

https://imgur.com/a/2h3I558 - and here after the spatial join there are 138 for the buffer from that Object.


texts:
I am doing Mann-kendall test in R using follow code but unexpectedly facing an error

setwd("F:/.../rasters/")
getwd()

rlist = lapply(1:111, function(i){raster(matrix(runif(20),5,4))})

rasters=stack(rlist)
rasters
fun_kendall <- function(x){ return(unlist(MannKendall(x)))}
kendall_result <- calc(rasters,fun_kendall)

and error is

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function

If anyone knows why it is giving this error and how to get rid off this


texts: Understanding and using fieldmapping with ArcPy?

texts:
I'm having trouble trying to understand how fieldmapping works in ArcPy (I'm kinda new to Python). In the below example, I want to take a shapefile and append it to a feature class. But I want to use the field map parameter so the fields from the shapefile go to the correct field names in the feature class. So I successfully built this out in ModelBuilder and exported the Python script (see below). But so much of it doesn't make sense, for example...what are all these different codes... -1, 0, #, true, false, etc. And there seems to be no documentation anywhere online on the meaning behind these codes.

Does anybody know what these mean?

Also, I read that I can use the fieldmapping() class as a parameter in this script.

How would I even incorporate that?

If you know a much easier way (via python) to field map then using the above methods, I'm open to those suggestions.

Edit: The code below is only a simple test that I created, so I can hopefully understand fieldmapping better. But in the long run I want to this to be able to handle a LOT of fields from many different shapefiles.

Source shapefile schema: Name, Code, Date

Target feature class schema: Description, Num, Date

# Import arcpy module
import arcpy

# Local variables:
New_Shapefile_shp = "C:\\New_Shapefile.shp"
New_Featureclass__2_ = New_Shapefile_shp
New_Featureclass = "C:\\Target.gdb\\New_Featureclass"

# Process: Append
arcpy.Append_management("C:\\New_Shapefile.shp", New_Featureclass, "NO_TEST", "DESCRIPTION \"DESCRIPTION\" true true false 150 Text 0 0 ,First,#,C:\\New_Shapefile.shp,Name,-1,-1;NUM \"NUM\" true true false 4 Long 0 0 ,First,#,C:\\New_Shapefile.shp,Code,-1,-1;DATE \"DATE\" true true false 8 Date 0 0 ,First,#,C:\\New_Shapefile.shp,Date,-1,-1", "")


texts: QGIS Automatically Generate Atlas for each Column on the Attribute Table

texts:
I want to generate graduated maps in QGIS (using the same style) for each of the columns on the attribute table. Every columns contains the data of one year.

It seems that the Atlas only automatically create maps for each row, it there any way to do it for the columns instead? I saw in another post that says it could be done by transpose, but how to I keep the graduated styles that's applied to each column?

My table look like this:

State          year1        year2       year3

state1         222          222         222

state2         333          333          333

state3         444          444          444


texts: QGIS raster calculator splitting resulting image in two

texts:
I'm trying to calculate vegetation indices via the QGIS raster calculator using two separate bands of a Sentinel-2 tile (both of which have the same extent and CRS) but the resulting image is split in two. The result has the same extent as the input images and reference layer but appears to be split between two horizontally compressed images. I've attached the result below.

This doesn't happen with any other Sentinel-2 tile I have used for the same purpose using the same raster calculator expressions.
Any help would be appreciated!


texts: GDAL point ordering (difference in version 2.2 and 3.0) Python

texts:
I'm currently trying to update some legacy Python code from gdal version 2.2.3 to 3.0.2, (along with python 2 > python 3)
I need to generate contours and unpack the points, along the lines of:
_, output = gdal.ContourGenerate( .... )
points_list = [contour.GetGeometryRef().GetPoints() for contour in output]

Now points_list contains the same very similar points for each contour, regardless of version. However the ordering is changed between 2.2.3 and 3.0.2
Python list comprehension preserves the order provided by the object being iterated. So unless GetPoints is unordered, I don't think this is a python issue.
Can anyone suggest a way to dictate the ordering generated by GetPoints? Or point me at the change log for this behaviour, so I can figure out what's happened?


texts: Changing NODATA value GeoTIFF with ArcPy of ArcGIS Pro?

texts:
I'm using a tool in ArcPy for Python 3 (ArcGIS Pro) which applies an operation to a raster file in GeoTIFF format.

condicion = '("{}"  - {}) / {}'.format(input_image,n,m)
arcpy.gp.RasterCalculator_sa(condicion, output_image)

Where n and m are doubles and input_image and output_image strings referring to the GeoTiff file.

It works fine but it changes the NODATA value from -3,40282306074E+38 (NODATA value used by ArcPy in Python2 and ArcMap) to -1.79769313486231571e+308 (nodata value used by ArcPy in Python3 and ArcGIS Pro).

But I want to preservate the previous NODATA value, I'm trying to use arcpy.SetRasterProperties_management tool. Is it possible to use it this way?

condicion = '("{}"  - {}) / {}'.format(input_image,n,m)
arcpy.gp.RasterCalculator_sa(condicion, output_image)
arcpy.SetRasterProperties_management(output_image, nodata='1 -3,40282306074E+38')

This code does not return any errors, but neither does the raster NODATA value change.


texts: Populating combobox with layers in TOC using PyQGIS?

texts:
I am developing a plugin, and I want user to select one of the loaded layers from a combobox. I've already coded this much:

self.dlg.show()
    layerlist = []
    self.dlg.ui.layerCombo.clear()  # clear the combo
    layermap = QgsMapLayerRegistry.instance().mapLayers()
    for name, layer in layermap.iteritems():
        if layer.type() == QgsMapLayer.VectorLayer:
            if layer.geometryType() == QGis.Line:
                layerlist.append( unicode( layer.name() ) )
    for layer in layerlist:    # foreach layer in legend
        self.dlg.ui.layerCombo.addItem( layer )

First Question: I get the error

  IndexError: list index out of range

which I don't quite know why.

Second Question: Am I doing the right thing? I couldn't find a tutorial on this. Where is the best place to put this code (run, initGui,...)?


texts: Exporting QGIS contour attributes in DXF file?

texts:
I have a .tiff data file that generates contours in QGIS. I save it as .dxf. The contours don't have any attributes or data. That's bad, because I want to make a 3D landscape in Revit using the information.

How do a save/export the contour lines from QGIS so that they have a z-elevation attribute or at least some information?

I ticked the attributes box when extracting the contours, but the attributes doesn't show in the dxf.


texts: How can I migrate from agslib 2.1 to agslib 3.1

texts:
I have an application which is developed using agslib 2.1. I want to migrate t to agslib3.1. When I use the new library, There is an error in MapManager.mxml code as "Call to a possibly undefined method zoomToFullExtent through a reference with static type com.esri.ags.tools:NavigationTool."
I could not find any updated MapManager.mxml in FlexViewer version 3.1. So I was thinking that maybe I must follow a process for migrating from 2.1 to 3.1. does anybody know how I can migrate my code to use the new library, or any workaround for resolving the issue?


texts: Displaying representation symbology (lines with markers) in ArcMap legend?

texts:
I have a polyline feature class with numerous cartographic representations of line types. Some of these representations only have a line component (see top-half image below) while other representations also have a marker component that repeats along the line (see bottom-half image below). In this case, a question mark repeats across the line type of an uncertain geologic boundary.

The representation/symbology appears correct here in the Symbology tab. But when I create a legend for this feature class, the marker part of the representation does not appear.

I thought maybe the default "patch" displayed in the legend was too short, and if a long enough patch were displayed then the marker component of the representation would be placed. But when I try to change the Default Patch Width in the Legend Properties > Layout tab, the line length in the legend remains the same and only the space between the patch and the label increases.

How do I properly display all the components of my representations in my legend (the line and the markers on the line)?


texts: What's the PostGIS equivalent to Extract Nodes in fTools / QGIS?

texts:
I'm looking to convert Park polygons from OSM into points.

The QGIS fTools Extract Nodes tool will convert a polygon to points at the vertices.

What is the equivalent function or work-flow in PostGIS?


texts:
I'm trying use the GDAL library in C++ to read a MSSQL spatial database file (*.mdf).  I see there is a driver for MSSQL Spatial, but there is very limited documentation on how to connect and read from the database.  I would like to read features from the database and tables in general.  I am familiar with reading shapefiles and rasters, but not reading SQL databases.

I would guess you create a connection string, load the GDAL Driver, but then I'm not sure where to go from there. Perhaps execute a SQL query to get a layer.  From the layer get all the features.


texts: Uploading shapefile to GeoServer using cURL and replace file?

texts:
In the sequence of this:

How to upload a shapefile to geoserver using cURL?

which was properly solved, I now need to be able to replace the existing file, cause my process is to be repeated every day.

Any clues?


texts: Converting point vector with many fields to multi-band raster

texts:
I have a shapefile of evenly spaced points which I'd like to convert to a raster.
In QGIS, I know I can do this with Raster > Conversion > Rasterize (Vector to Raster).
However, I see the possibility to use only one of the fields as a "burn-in value", while I'd rather keep all my fields (at least the numeric ones). I guess I'd like a multi-band raster where each band corresponds to each (numerical?) field in my shapefile.
Is there a way to do that?
My aim would be to be able to symbolize this raster based on one of the many variables (bands), in a way that I can easily switch to another variable and symbology as I need, without having to create different raster for different variables.


texts: Custom Overlay Grid - Labeling the axes differently

texts:
i have a map where i want to put the grid axes. i have a separate layer based on which i want to create the axes. now custom overlay feature in arcGIS gives you only one label field to label the axes but i have 2 columns based on which i want to showcase the X and Y axes labels.

How do i go about doing it?

If there is any other way other than custom overlay, please let me know of the same as well.

Note that there is already a grid layer which i am using the information for. I have attached the sample image for reference purpose and the grids shown in the image is the layer i am using as reference for custom overlay.


texts: Getting active PostgreSQL connection username when stored through authentication configuration in QGIS 3.16

texts:
How do I get the username on the current layer PostGIS in QGIS 3.16.2-Hannover?
We tried examples: Getting active PostgreSQL connection username when stored through authentication configuration in QGIS 3, but they don't work.
When running this script, variable username is empty:
layer = iface.activeLayer()
provider = layer.dataProvider()
username = QgsDataSourceUri(provider.dataSourceUri()).username()
print (username)

This is result of provider.dataSourceUri()):

<QgsDataSourceUri: dbname='postgres' host= localhost port=5432
sslmode=disable key='id' srid=4326 type=Point
checkPrimaryKeyUnicity='1' table="locations"."airlines" (geom)>

The value of the username is missing, the question is how to get the username for the current layer?


texts: Removing some bands or layers from a tiff using gdal

texts:
I have been working on a multiband raster image. It is a pix image. I need to just work on the 26 bands that are raster files and not the 4 bitmaps that are are also included in the pix image. I could force the loop that will be going through it to only go through 25 times but that wont work when the pix images change and the number of bands change. Right now I am using gdal and have it all in a 3D array. I am pretty sure the bitmaps would be the last five layers of the array. I do not have anything to tackle this problem so I have no code.

Basically I just need to delete 4 layers from the pix image and work with 26 layers. When I call raster_image.rastercount, instead of it printing 30, I want it to print 26. I have gone through the gdal library, I may have missed something but I couldn't find anything. Is there a way to delete entire bands using gdal?

I have attached a screenshot of how the bands looks and what one of the bitmaps look like when displayed. Is there anyone that can help? Or just tell me it is impossible.

I am trying to use the GetGeomType just to see how it works and if it will work for me. But it does not seem to work:

    raster_image = gdal.Open(in_raster)
    band3 = raster_image.GetGeomType(3)
    print  band3

Above is how I am calling it but it gives me the error:

  AttributeError: 'Dataset' object has no attribute 'GetGeomType'

Then trying it by first getting the raster band and then checking the geom type I get:

    band3 = raster_image.GetRasterBand(3).GetGeomType()
    print band3

  AttributeError: 'Band' object has no attribute 'GetGeomType'


texts: Visualise different types of areas within a city polygon (like a stacked area plot) using ArcGIS

texts:
I would like to visualize different types of areas (e.g. forest, river, industrial areas, ...) according to their proportion within a city using
ArcGIS.

My goal is to "fill" my city boundary polygon from south to north just like a sandglass by simply using the size of each type of area. In the end it should somehow look like a stacked area plot "under the influence of gravity"  considering the city boundary as in the following sketch:

Is this possible? Any hints?


texts: Assignment of defaults to field not working

texts:
I am using ArcGIS 10.2.2 for Desktop.
I have assigned default values to some fields in a feature class, which is sitting in a file geodatabase.

I used the "Assign Default To Field Tool", which claims to have been successful, yet when I go to create a new feature none of the default values populate in the fields.
When I check the field properties it shows the default value that I assigned, it just doesn't use them.
I do have domains assigned to the fields in question, and the default values I have chosen do conform to the coded values within the domains.

I am also using a subtype within the feature class, but I can't see why that would be interfering.
I think I have covered all relevant information, just ask if I have missed something.


texts: Calculating mean from multiple folder using ArcPy

texts:
I am writing a script to loop through multiple folders and calculate mean of the raster of each folder's data sets and save it in same folder. Below script i am working, but its not producing  the mean rasters for each year folder, after running the script its generating only name of "mean" raster in main folder, and overlapping it .

Below images showing the name of the raster in one folder, another folder the raster name is same, only year value is changed in different year.
here how to take the year value in name of mean raster (E.g. for 2001 folder, output mean raster name would be Mean_RAIN_2001.tiff)

Script:

import arcpy, os
from arcpy import env
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
env.workspace = r"D:\Test"
outraster = env.workspace

walk = arcpy.da.Walk(env.workspace, topdown=True, datatype="RasterDataset")
for dirpath, dirnames, filenames in walk:
    print dirpath
    rasterList = arcpy.ListRasters("TIF")
    print rasterList

    for Rasters in rasterList:
        rasMean = CellStatistics(Rasters,"MEAN", "DATA")
        file_name_only = os.path.splitext(Rasters)[0]
        tifname = file_name_only[-4:]
        rasMean.save(os.path.join(outraster,'Mean_{0}.tif'.format(tifname)))
        print rasMean

Printing statement: its printing only folder name,

D:\Test
[]
D:\Test\2001A
[]
D:\Test\2002A
[]
D:\Test\2003A
[]
D:\Test\2004A
[]
D:\Test\2005A
[]
D:\Test\2006A
[]
D:\Test\2007A
[]
D:\Test\2008A
[]
D:\Test\2009A
[]