texts: Can't find libpq when installing PostGIS

texts:
I'm trying to install PostGIS 2.2.1 on Mac OS X 10.8.  I've installed all of the prerequisites for PostGIS and am now running the configure command:

./configure CC=/usr/bin/clang CXX=/usr/bin/clang++ "CFLAGS=-arch x86_64 -arch i386" "LDFLAGS=-arch x86_64 -arch i386" "CXXFLAGS=-arch x86_64 -arch i386"

I'm getting the following error:

configure: error: could not find libpq

I have PostgreSQL 9.1.4 installed and I believe I see the libpq libraries which is owned by root:wheel:

/usr/local/pgsql/lib/libpq.5.4.dylib
/usr/local/pgsql/lib/libpq.5.dylib
/usr/local/pgsql/lib/libpq.a
/usr/local/pgsql/lib/libpq.dylib

I did the following but it didn't help:

export LD_LIBRARY_PATH=/usr/local/pgsql/lib

This too in case PostGIS just needs to know where Postgres is installed:

export LD_LIBRARY_PATH=/usr/local/pgsql

Finally, I ran this command with and without sudo:

./configure CC=/usr/bin/clang CXX=/usr/bin/clang++ "CFLAGS=-arch x86_64 -arch i386" "LDFLAGS=-arch x86_64 -arch i386" "CXXFLAGS=-arch x86_64 -arch i386" "LD_LIBRARY_PATH=/usr/local/pgsql/lib"

What am I doing wrong?  I found another post here where someone else had the same problem but the solutions noted there didn't work for me.

Thanks!

UPDATE:
I'm looking in the PostGIS config.log file and can at least see why the error is occurring.  I just don't quite understand what to do to prevent it.

ld: warning: ignoring file /usr/local/pgsql/lib/libpq.dylib, file was built for unsupported file format ( 0xcf 0xfa 0xed 0xfe 0x 7 0x 0 0x 0 0x 1 0x 3 0x 0 0x 0 0x 0 0x 6 0x 0 0x 0 0x 0 ) which is not the architecture being linked (i386): /usr/local/pgsql/lib/libpq.dylib
Undefined symbols for architecture i386:
  "_PQserverVersion", referenced from:
      _main in conftest-sxaGdF.o
ld: symbol(s) not found for architecture i386
clang: error: linker command failed with exit code 1 (use -v to see invocation)


texts: QGIS 3.10 problem with merging layers & multipart promotion

texts:
I would like to merge a few layers in QGIS. However, I can't do this.
When I run the MMGIS option, then I am getting:
Merged input_layers must all be same type of geometry (MultiLineString != Polygon)
According to the MMQGIS author, all layers must be in the same geometry.
In this case, I set all the same geometries for my objects:

But it didn't solve my problem:

So I tried another way:
In the Vector - Data Management Tools - Merge Vector Layers  I selected my layers. Unfortunately, I got sth like this:

According to this question:
Merge layers with different geometry
I used "Join attributes by location" option, although afterward, I got only 1 layer, even when I select "Create the separate features for joint..."

Should I use the GRASS as per below?
Error in merging two shapefiles of different geometry types
I would like to have all the layers in the picture below merged:

Or should I maybe do sth like this:
Merging layers into one- QGIS
UPDATE:
As per the advice below I tried the "Promote to the multipart" option.~
Unfortunately, it doesn't solve my query too, since this tool can promote only 1 layer:

In the effect I have only 1 layer "multipart", but how about other ones?


texts: Why doesn't QGIS georeferencer support gdal scripting for projective transformation?

texts:
I am georeferencing drone images with QGIS, and the transformation method that gives the best results is the projective one.

As I have many images, I would like to write a gdal script to automate the process.

However, when I hit the "generate GDAL script" button, I get the following error:

  GDAL scripting is not supported for projective transformation.

Is it gdalwarp which doesn't support projective transformation, or is QGIS georeferencer that doesn't implement it?


texts: How to generate Tiles with QuantumGIS or Mapnik

texts:
I would like to know how to generate Map Tiles using QGis. I have heard that it is possible to create tiles using QGis, only, I don't know how to do that. Googling did not work for me since I am not really expertised with programming (Python etc.) and GIS.

However, I am experienced with accomplishing tasks that I am not experienced. So, if somebody can direct me to a tutorial/guide (detailed tutorial is preferred), I can manage to generate tiles.

I can also use Mapnik but as far as I understand, Mapnik does not have a GUI so I prefer QGis.


texts: Dissolve Statistics Fields and Statistic Type as Parameters in GeoProcessor

texts:
I have a published service which uses a dissolve.  In arcMap the I get prompted for all the fields, including statistics fields and stat type.  I have the statistics fields marked as parameters.

However, when I publish the service to arcgis for server 10.1 all parameters show except the stats fields in the parameter list on server.  When I try to run the service through the javascript api it uses the field names and types I used just prior to publishing it.  Which, of course, works fine if I pass features with the same name.

Is it possible to publish a service with a dissolve which makes the stats fields and type a parameter that is accessible by arcgis javascript api?


texts:
I want to change the CRS of a shapefile which uses EPSG:6257. Unfortunately, I am using Ubuntu 20.04 which comes with PROJ6. This particular CRS was only made available from PROJ7.2 onwards. So, I am now exploring my options to use a custom CRS in Proj. According to the documentation:

Coordinate reference system xyz is not in the EPSG registry, what do I do?
Generally PROJ will accept coordinate reference system descriptions in
the form of WKT, WKT2 and PROJ strings. If you are able to describe
your desired CRS in either of those formats there’s a good chance that
PROJ will be able to make sense of it.

So, I have been trying to do just that, but I struggle to make this work. I'm using Django's geodjango library, where I can use SpatialReference to set a CRS. According to that documentation:

class SpatialReference(srs_input)¶
Spatial reference objects are initialized on the given srs_input, which may be one of the following:

    - OGC Well Known Text (WKT) (a string)
    - EPSG code (integer or string)
    - PROJ string
    - A shorthand string for well-known standards ('WGS84', 'WGS72', 'NAD27', 'NAD83')

So I have been trying to set this in a number of ways, but nothing works.
If I use the "WKT Version 2" from the EPSG website:
crs = 'PROJCRS["MAGNA-SIRGAS / Medellin urban grid",BASEGEOGCRS["MAGNA-SIRGAS",DATUM["Marco Geocentrico Nacional de Referencia",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7019]],ID["EPSG",6686]],PRIMEM["Greenwich",0,ANGLEUNIT["radian",1,ID["EPSG",9101]],ID["EPSG",8901]],ID["EPSG",4686]],CONVERSION["Medellin urban grid",METHOD["Colombia Urban",ID["EPSG",1052]],PARAMETER["Latitude of natural origin",6.229208889,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]]],PARAMETER["Longitude of natural origin",-75.564886944,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]]],PARAMETER["False easting",835378.647,LENGTHUNIT["metre",1,ID["EPSG",9001]]],PARAMETER["False northing",1180816.875,LENGTHUNIT["metre",1,ID["EPSG",9001]]],PARAMETER["Projection plane origin height",1510,LENGTHUNIT["metre",1,ID["EPSG",9001]]],ID["EPSG",6225]],CS[Cartesian,2,ID["EPSG",4500]],AXIS["Northing (N)",north],AXIS["Easting (E)",east],LENGTHUNIT["metre",1,ID["EPSG",9001]],USAGE[SCOPE["Cadastre, engineering survey, topographic mapping (large scale)."],AREA["Colombia - Medellin city."],BBOX[6.12,-75.66,6.37,-75.50]],ID["EPSG",6257]]'

ct = CoordTransform(SpatialReference(medellin), SpatialReference("WGS84"))

Then I get:

django.contrib.gis.gdal.error.GDALException: Corrupt data.

medellin = 'GEOGCRS["MAGNA-SIRGAS",DATUM["Marco Geocentrico Nacional de Referencia",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7019]],ID["EPSG",6686]],PRIMEM["Greenwich",0,ANGLEUNIT["radian",1,ID["EPSG",9101]],ID["EPSG",8901]],ID["EPSG",4686]]'
ct = CoordTransform(SpatialReference(medellin), SpatialReference("WGS84"))

Returns the same:

django.contrib.gis.gdal.error.GDALException: Corrupt data.

Trying another format from the shapefile itself:
medellin = """GEOGCS["GCS_MAGNA",
    DATUM["Marco_Geocentrico_Nacional_de_Referencia",
        SPHEROID["GRS_1980",6378137.0,298.257222101]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]"""

Unsupported SRS.

Any idea how I can use this custom CRS that is not present in my version of PROJ?


texts:
In QGIS 2.18.14 Processing Toolbox I started to use v.buffer.distance to create fixed-distance buffers around my lines.

I can also make the buffer with square caps (straight ends).

However, I cannot make it without end cap. I mean, the option "Do not make caps at the ends of polylines" does not seem to work.

This is similar to this recent question v.buffer.distance does not create straight buffers in QGIS, but my situation is different, as I can make square caps (-s) while not being able to cut caps (-c).
(Probably @Annekeds is using QGIS - GRASS plugin, while I am working through Processing Toolbox).

Strangely I can make buffer without caps when I run v.buffer.column... so my workaround for now is to use v.buffer.column with same distance in the field.


texts:
My recent assignment is to model a multi-user geodatabase for the whole world using the ESRI platform. This geodatabase contains administrative layers and topographic layers. The schema is to be designed in such a way that users can update different countries simultaneously.

My first problem is that every country has different administrative structures. For example, Country A has 2 administrative levels, Country B has 3 administrative levels and Country C has again 2 administrative levels.

Now when I merge each set of admin levels into layers (i.e. all level 1 areas in one layer, level 2 in a second layer and level 3 in a third layer) and establish a topology there is a gap between Country C and Countries A and B in admin level 3 because that level is not available in Countries A and B. The goal is for the topology to not have a gap.

Am I organizing this data in the wrong way?


texts:
I know osm2po is great for getting an osm file into a database.  I used it to load an entire country in a country_2po_4pgr table.

What I would like to do is then reverse the process, to be able generate an osm2po .gph file based on an SQL query to country_2po_4pgr of a certain polygon inside the country.

Is that possible?


texts:
i have a table with two columns: one named index and the other is called tilenumber and the tilenumber is populated from the index column except that there is _ which needs to be removed from the column.

The index table has values like 2726BA_11 and the tilenumber should result in 2726BA11 but i seem not to know the way forward, to remove the space i have done the following:

select distinct replace (tilenumber,'_', '') from aug2012

and it shows that the result has been achieved but how do i incorporate the select distinct replace in the update statement which is

update set tilenumber="index"


texts: Export georeferenced image from AutoCAD to GIS?

texts:
I have an AutoCAD .dwg file that was assembled from scanned maps and is georeferenced.

I would like to export it to use as a raster in ArcGIS.
What file format should I use, that can be opened in ArcGIS and will maintain the spatial information?
How can I do the conversion?

I have Autocad and Raster Design. I also have the original images, which are .tif files
The "CAD to geodatabase" tool does not work for this.


texts: OGR Filter Vector Layer by a List of Feature ID's

texts:
I often use GDAL/OGR via their python bindings.

I know of the very useful ogr.Layer.SetAttributeFilter() command. However I am hoping to filter by feature ID's, rather than attributes. Also, I don't want to filter by just one value, but rather by a set of values. Ideally I could pass a list of feature ID's to the filter -- is there some method that can allow me to do this (either through SetAttributeFilter or otherwise)?

It occurred to me that I could copy the feature ID's to their own attribute field, but it would be nice to avoid having to do this.

CLARIFICATION: I am hoping to apply filters for read-only processes which would involve sorting, selecting & exporting/rasterizing the layer geometries. I primarily use ESRI Shapefile format.


texts: Intersecting with Geodatabase.ISpatialFilter - no results

texts:
I'm trying to intersect two layers and get resulting features with the following code I've written:

   ESRI.ArcGIS.Geodatabase.ISpatialFilter spatialFilter = new ESRI.ArcGIS.Geodatabase.SpatialFilterClass();
   spatialFilter.SpatialRel = ESRI.ArcGIS.Geodatabase.esriSpatialRelEnum.esriSpatialRelIntersects;

   int indOfLayer1 = Int32.Parse(ConvertFromLayerNameToLayerId(m_Map, m_ResourceMapItem, m_firstLayerName));
   int indOfLayer2 = Int32.Parse(ConvertFromLayerNameToLayerId(m_Map, m_ResourceMapItem, m_masterLayerName));

   // Get the map underlying the map service and the IGeoFeatureLayers contained in the map
   ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.MapResourceLocal mfunct =
                    (ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.MapResourceLocal)
   m_Map.MapResourceManagerInstance.GetResource(m_ResourceMapItem);
   ESRI.ArcGIS.Carto.IMapServer iMapServer = (ESRI.ArcGIS.Carto.IMapServer)mfunct.MapServer;

    ESRI.ArcGIS.Carto.IMapServerObjects mapServerObjects = (ESRI.ArcGIS.Carto.IMapServerObjects)iMapServer;
    ESRI.ArcGIS.Carto.IMap map = mapServerObjects.get_Map(iMapServer.DefaultMapName);
    ESRI.ArcGIS.esriSystem.UID layerTypeID = new ESRI.ArcGIS.esriSystem.UIDClass();
    layerTypeID.Value = "{40A9E885-5533-11d0-98BE-00805F7CED21}"; //feature layer
    ESRI.ArcGIS.Carto.IEnumLayer enumLayer = map.get_Layers(layerTypeID, true);
    enumLayer.Reset();

    // first layer
    ESRI.ArcGIS.Carto.IFeatureLayer featureLayer1 = map.get_Layer(indOfLayer1) as ESRI.ArcGIS.Carto.IFeatureLayer;

    // master layer (second)
    ESRI.ArcGIS.Carto.IFeatureLayer featureLayer2 = map.get_Layer(indOfLayer2) as ESRI.ArcGIS.Carto.IFeatureLayer;

    spatialFilter.Geometry = (IGeometry)multipointPoints2;

    // perform the query
    ESRI.ArcGIS.Geodatabase.IFeatureCursor cursor = featureLayer1.Search(spatialFilter, true);

    ESRI.ArcGIS.Geodatabase.IFeature feature = null;
    // Create a ComReleaser for cursor management.
    IGeometryCollection formetotali = new MultipointClass();
    ESRI.ArcGIS.ADF.ComReleaser comReleaser = new ESRI.ArcGIS.ADF.ComReleaser();
    // Create and manage a cursor.
    comReleaser.ManageLifetime(cursor);
    object missing = Type.Missing;
    // Delete the retrieved features.

    while ((feature = cursor.NextFeature()) != null)
    {
        ESRI.ArcGIS.Geometry.IGeometry mageometry = feature.Shape;
        formetotali.AddGeometry(mageometry, ref missing, ref missing);
    }

    DrawIGeometryOnMap(m_Map, (IGeometry)formetotali);

I'm getting no results at all, even if the layers intersect. What's wrong with this code?


texts: Using large data JSON file in Leaflet markercluster

texts:
I want to show the JSON file consisting of many points on the map using the leaflet marker cluster plugin.
My JSON file looks like:
[
{
    "json_featuretype" : "Output",
    "ILCE_UAVT" : 1421,
    "ILCE_ADI" : "KADIKÖY",
    "ADI" : "Marmara Veteriner Kliniği",
    "ALT_KATEGORI" : "Veteriner",
    "ADRES" : "Bahçeler Sk. No:3 FENERYOLU/KADIKÖY",
    "TELEFON" : "0216 347 49 42",
    "WEBSITESI" : null,
    "ACIL_SERVIS" : null,
    "YATAK" : 0,
    "AMBULANS" : "Yok",
    "MAHALLE" : "FENERYOLU",
    "BOYLAM" : 29.05070181700006,
    "ENLEM" : 40.979645135000055,
    "json_ogc_wkt_crs" : "GEOGCS[\"WGS84 Lat/Longs\",DATUM[\"WGS_1984\",SPHEROID[\"World Geodetic System of 1984, GEM 10C\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]",
    "json_geometry" : {
        "type" : "Point",
        "coordinates" : [ 29.050701817, 40.979645135 ]
    }
},
 /////too many points information
{
    "json_featuretype" : "Output",
    "ILCE_UAVT" : 1327,
    "ILCE_ADI" : "FATİH",
    "ADI" : "İstanbul Üniversitesi-Cerrahpaşa Tıp Fakültesi Kulak Burun Boğaz Ana Bilim Dalı",
    "ALT_KATEGORI" : "Poliklinik",
    "ADRES" : "İÜ Cerrahpaşa Tıp Fakültesi Yolu No: CERRAHPAŞA/FATİH",
    "TELEFON" : "0212 414 23 59",
    "WEBSITESI" : "http://istanbultip.istanbul.edu.tr",
    "ACIL_SERVIS" : "Var",
    "YATAK" : 60,
    "AMBULANS" : null,
    "MAHALLE" : "CERRAHPAŞA",
    "BOYLAM" : 28.94229438600007,
    "ENLEM" : 41.00354611200004,
    "json_ogc_wkt_crs" : "GEOGCS[\"WGS84 Lat/Longs\",DATUM[\"WGS_1984\",SPHEROID[\"World Geodetic System of 1984, GEM 10C\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]",
    "json_geometry" : {
        "type" : "Point",
        "coordinates" : [ 28.942294386, 41.003546112 ]
    }
}]

Can I do this demonstration with the leaflet marker cluster with this json file?


texts: Is there a way to show outlines in a multipatch feature in ArcScene?

texts:
When viewed in ArcMap 2D environment, multipatch faces are outlined. When viewed in ArcScene 2.5 D environment, multipatch faces are not outlined, yet the symbol in symbology properties still shows an outline.
Is there a way to show the multipatch faces outlines when viewed in ArcScene?

ArcScene:

ArcMap:


texts: Receiving two different values for same date for Landsat B8 using Google Earth Engine

texts:
I am extracting values of the B8 of Landsat8 in Google Earth Engine and for several month I receive two very different values for the same date and cannot find the issue.
How can I handle this?

My code:
var maskL8 = function(image) {
  var qa = image.select('BQA');
  /// Check that the cloud bit is off.
  // See https://www.usgs.gov/land-resources/nli/landsat/landsat-collection-1-level-1-quality-assessment-band
  var mask = qa.bitwiseAnd(1 << 4).eq(0);
  return image.updateMask(mask);
}

var start = '2020-07-01';
var finish = '2020-07-30';
var pt = ee.Geometry.Point([-49.31582, 69.56833]);
// Map the function over one year of Landsat 8 TOA data and take the median.
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterDate (start,finish)
.filterBounds(pt)
.map(maskL8);

//Map.addLayer(B8);
//Map.addLayer(pt);
Map.centerObject(pt, 16);

print(l8);

// B8 VALUE OF PT

var myB8 = l8.select("B8");
print("myB8",myB8);

var getB8 = function(image) {

  // Reducing region and getting value
  var value_B8 = ee.Image(image)
    .reduceRegion(ee.Reducer.first(), pt)
    .get('B8');

  return value_B8;
};

var count = myB8.size();

var B8_list = l8.toList(count).map(getB8);

print("B8 list", B8_list);

var img = l8.first();

Map.addLayer(img);
Map.addLayer(pt);

var allDates = l8.aggregate_array('system:time_start');

var allDatesSimple = allDates.map(function procDates (ele) {

  return ee.Date(ele).format().slice(0,10);

});

var paired = allDatesSimple.zip(B8_list);

print (paired);


texts: Creating lines from points pair coordinates with ArcPy?

texts:
I've got some points pair coordinates (start and end points) which I have to transform into lines. Until now, I used an append of both the coordinates in a pippo.Point(), a pippo.CalculateGeometry() to define the geometry of each piont, and pippo.append(defined geometry) to identify the pair of points, and then PointsToLine to obtain my line. This is quite time expensive to do for hundreds of lines.

Is there a shorter way to do this?

For example, place starting and ending point of each line in different fields of a single table and import lines directly without passing for points geometry.


texts: How to merge datasets without creating new columns

texts:
I am using Carto Builder and I have two datasets with exactly the same column headings but each dataset has different data under each column heading.

Is there a way to merge the data without creating a new column heading each time. I've been using the Merge Datasets function but can't find a way to avoid creating a new column each time.

As you can see they have the same column names and I just want to merge the data below all into one whole table.


texts: Feature service misalignment in ArcGIS Online

texts:
I have an issue with feature layers in ArcGIS Online. I am adding an ArcGIS Server web service witch creates a feature layer and then I add the feature layer in a web map. There are 11 layers in the feature layer. When I try to open the webmap I get a message that a layer is not responding

This is a water network map. As you can see the pipes are misaligned and some pipes are missing.
When I use the Collector app and download the map for offline use, they get aligned properly and there is nothing missing.

On another instance if I use a map image Layer everything shows up OK.
The spatial reference is Greek Grid (2100)

I am using ArcMap 10.4, ArcGIS Server 10.4, the service details are as following (for one of the 11 layers)

Layer: ΔΙΚΛΕΙΔΕΣ ΠΕΖΟΔΡΟΜΙΟΥ (ID: 0)
View In:   ArcGIS Online map viewer
Name: ΔΙΚΛΕΙΔΕΣ ΠΕΖΟΔΡΟΜΙΟΥ
Display Field: FACILITYID
Type: Feature Layer
Geometry Type: esriGeometryPoint
Description: This dataset contains a collection of features that represent        the water distribution network in a community.
Definition Expression: N/A
Copyright Text:
Default Visibility: true
MaxRecordCount: 5000
Supported Query Formats: JSON, AMF, geoJSON
Min Scale: 2000
Max Scale: 0
Supports Advanced Queries: true
Supports Statistics: true
Use Standardized Queries: true
Supports ValidateSQL: true
Supports Calculate: false
Extent:
XMin: 350329.08951833844
YMin: 4373931.874201424
XMax: 373356.6671507545
YMax: 4398354.849961676
Spatial Reference: 2100  (2100)

Drawing Info:Renderer:Simple Renderer:
Symbol:Picture SymbolLabel:
Description:Transparency: 0
Labeling Info:Label Placement: esriServerPointLabelPlacementAboveRight
Label Expression: [code]
Use Coded Values: true
Symbol:Color: [0, 0, 0, 255]
Background Color: [233, 255, 190, 255]
Outline Color: [0, 0, 0, 255]
Vertical Alignment: bottom
Horizontal Alignment: left
Right to Left: false
Angle: 0
XOffset: 0
YOffset: 0
Size: 8
Font Family: Arial
Font Style: normal
Font Weight: normal
Font Decoration: none
Min. Scale: 500.0
Max. Scale: 0.0
Advanced Query Capabilities:
Supports Statistics: true
Supports OrderBy: true
Supports Distinct: true
Supports Pagination: true
Supports TrueCurve: false
HasZ: false
HasM: false
Has Attachments: false
HTML Popup Type: esriServerHTMLPopupTypeAsHTMLText
Type ID Field:
Fields:
too many to display
Templates:
too many to display
Drawing Tool: esriFeatureEditToolPoint
Capabilities: Create,Query,Update,Sync,Delete,Uploads,Editing
Sync Can Return Changes: true
Is Data Versioned: true
Supports Rollback On Failure: false
Supports ApplyEdits With Global Ids: false
Supported Operations:   Query   Apply Edits   Add Features   Update Features         Delete Features   Validate SQL   Generate Renderer   Return Updates

Edit:
The misalignment is not constant in all layers.

For example as shown in the image below the main pipes(the last layer) are misaligned whereas the secondary pipes (second from last) are aligned properly.

Does anyone have a clue?


texts:
I'm looking to create an APP with Navigation Offline.

The problem is my point of interest is not offered by OpenStreetMap (OSM).

I would like to know if can i add navigation to a shapefile (mine) or should i add the new area to OSM?


texts: Allowing others to add points to QGIS Cloud after publishing (or notes)?

texts:
How can I allow others to add points in a layer published on QGIS Cloud?
I have a map of my city and was hoping myself and friends could add notes from our phones at different locations.


texts: Clipping .tif with gdal Warp produces blank image

texts:
I have four band .tif that I am trying to clip to the shapefile extent with:

gdal.Warp(OutDS, inRaster, cutlineDSName = InShapefile1, cropToCutline = True, dstNodata = inNoData)

using Python and GDAL I first open the raster and then try to clip it band by band. The result has only zeros. Clipping shapefile has only one polygon geometry but it is possible that it extends out of raster's bounds in some cases. I've tried to perform clip with color composite but result is the same. Also the output single band rasters are read by, for example, ArcMap as multiband composites.

The link to the sample data is here: https://dl.dropboxusercontent.com/u/83030488/Tiff/Sent10.rar

UPDATE

By implementing the subprocess, call method I managed to identify the error. The problem was that I did not have the system veriable GDAL_DATA set correctly.

Now another error came out:

ERROR 4: `EPSG:3912' does not exist in the file system, and is not recognized as a supported dataset name.

I am calling gdalwarp with:

path= str(outClipedRasterPath)
os.chdir(path)
cmd = 'gdalwarp -dstnodata -9999 -cutline -s_srs EPSG:3912 -t_srs EPSG:3912 ' + InShapefile1 + ' -crop_to_cutline -of GTiff ' + inRaster + OutDS
call(cmd)


texts: Creating a Line at a Length and Angle in WGS84 - Arcobjects

texts:
I need to draw a line at desired length(around 10-20 Km length) and a user entered angle. The map/data uses WGS84 coordinate system. So the measurement is geodesic.

Here are the methods that I tried.

Find the point at an angle and distance.

Obviously this will not work as the computation is planar and this works only at the equator region.

IConstructGeodetic.ConstructGeodeticLineFromDistance method

I called it like below and this seemed to be a promising method as it gave correct distances. But at certain locations, this method simply crashes with outofmemoryexception

ConstructGeodeticLineFromDistance(esriGeodeticType.esriGeodeticTypeGeodesic, FromPoint, lUnit, d,Azimuth, esriCurveDensifyMethod.esriCurveDensifyByLength, 0.01)

Using the Equations found at movable-type website

This equation also has distance issues at different parts of earth. If i draw a 10Km line sometimes it comes as 10, but sometimes its 9 Km. The method I converted to c# is below.

public static IPoint GetAtDistBearing(IPoint co, double distance, double brng)
{
    double lat1 = co.Y;
    double lon1 = co.X;
    double ER = 6374;
    double d = distance;

    brng = ToRadians(brng);
    lat1 = ToRadians(lat1);
    lon1 = ToRadians(lon1);

    double lat2 = Math.Asin(Math.Sin(lat1) * Math.Cos(d / ER) + Math.Cos(lat1) * Math.Sin(d / ER) * Math.Cos(brng));
    double lon2 = lon1 + Math.Atan2(Math.Sin(brng) * Math.Sin(d / ER) * Math.Cos(lat1), Math.Cos(d / ER) - Math.Sin(lat1) * Math.Sin(lat2));

    lon2 = (lon2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI;

    IPoint p = new PointClass();
    p.X = ToDegrees(lon2);
    p.Y = ToDegrees(lat2);

    return p;
}

I think I am missing something obvious. Hope someone can point me in the right direction.


texts: Converting/expanding line into 2D polygon using QGIS

texts:
I am brand new to GIS and am attempting to use QGIS to take an existing line from a shapefile and create a 2D polygon based on it i.e. I want to take a line and fatten it out a bit so it has an area.
Specifically I have a shapefile which represents centrelines of NZ rivers available from NZ Land information data service. I have managed to select the river segments that I want, and have merge them into one line (not sure about that step being useful) and saved as it's own shapefile which I can load into QGIS. For example the lines within the area circled in red in this image:
I would like to take this line (or each separate segment) and fatten it/them out into a 2D polygon area similar to the circled area in this image:
My 2D polygon would have branches though like the first image, and I would like to have one continuous polygon for the entire outline. (I would also like to calculate the area in hectares, but I believe this should be fairly trivial once I have the polygon)
What would be the best method to achieve this? Is there an feature/plugin in QGIS which can automatically fatten a line by specified amount? Is there a programmatic method I can use (I am fairly proficient at Python programming but not very familiar with QGIS)?
I am not expecting a one shot kill and do expect to do some manual drawing to achieve this but am hoping there is some technique or tool to help me out before I go ahead and manually draw the entire thing.
EDIT:
After seeing the answers posted I should clarify that I am trying to basically draw the outline of the "wet bed" of the rivers/creeks and is for my purposes arbitrary to a certain extent. So after consideration this would clearly not correspond to a set distance from the centre lines.
After seeing @UdithaGIS answer and looking at the topographic elevation data it appears what I need is more closely related to the contour lines (but not exactly). I will try to use the contours as a starting point and attempt to modify them for my needs.


texts:
After reinstalling QGIS, the desktop program runs but generates the python error "Couldn't load SIP Module
Pythin support will be disabled"

Because of this error, no plug-ins will load.  I found a couple of descriptions of a similar problem on-line, which recommended:

sudo apt-get install python-sip4
sudo apt-get install python-qt4

Neither of these worked, and in fact neither even returned an installation candidate. Note that the error reports that recommended these were at least 1 year old, or older.

I installed SIP from this source http://www.riverbankcomputing.com/software/sip/download and confirmed its installation, but QGIS continues to generate the same error.

System details:

OS: Xubuntu 12.04 with Xcfe desktop
QGIS: 1.8
Python: 2.7.3
SIP: 4.13.3

Any help would be greatly appreciated.


texts: Sequentially renumber point feature based on distance along line in QGIS

texts:
I have a point layer with features that are along a separate line layer, and have calculated the points' distance along the line using line_locate_point(aggregate('Line','collect',$geometry),$geometry). The values are stored in a virtual field named loc in the point layer.

Now I want to be able to sequentially renumber the points along the line based on loc. How can I do this purely in QGIS, without creating a new layer (i.e. so not using something like Sort and Number)?

I tried using scale_linear() but this is not the appropriate tool (see results below - I would like testno to go 1,2,3,4...)

I'm aware there is an easy solution in PostGIS using rank() OVER PARTITION BY but I'm using spatialite and I don't think rank() exists. I can't get it to work.

The quick and dirty alternative is to save loc as a proper field, save as shapefile and renumber appropriately by opening the .dbf in Excel, but I'd prefer to find a solution that works entirely in QGIS without extra plugins if possible.


texts: The SetNull tool does not correctly change "-9999" into NODATA

texts:
I am using ArcMap 10.5.1 and .tif raster data. When I executed the SetNull tool to convert all "-9999" values into NoData, the range of values in the output raster was greater than that in the original raster.

The input raster's value range is (-9999,896.762), but the output value range is (0, 977.176). I don't know how did this operation generate values greater than 896.762.

My arcpy code stand-alone script:

import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("spatial")

inputRas = "C:\\Desktop\\Data\\input.tif"
outputRas = SetNull(inputRas, inputRas, "VALUE = -9999")
outputRas.save("C:\\Desktop\\Results.gdb\\output.tif")

I also tried the Con statement and the output values had the same messed-up range as well.

Any idea how to fix it?


texts: Create a geoserver SQL view from Oracle SQL using a lat and lng column

texts:
I have a table in a Oracle SQL database in which I have a polygon geometry and two more columns of latitude and longitude. What I need is to show the points of lat and lng in a layer in geoserver.

I have done similar thing using PostGIS but in now with Oracle I face an error.
First I have create a query which selects the two coloumns, transforms them to the desired SRID (EPSG:2100), REPLACE the '.' with the ',' and CAST them to float (cause they are characters).

When I execute the query in the Oracle SQL developer it seems that it works correctly. Also if I copy and paste teh coords they are shown in the correct position on the map.

But when I se my query in geoserver in order to create a sql view layer, although the geometry column is recognized (also the SRID), when I click on the "compute BBOX" option I get a geoserver error.

As I said I did something similar with postGIS and it was working fine. I also tried to use the same table (but its existing geometry; the polygon) and it works also great.

So I guess it must be related with my query but I can not guess what can cause this "unexpected geoserver error".

Thats my query:

SELECT ID, CO_NAME, ADJUSTMENT_FACTOR, PORT_LOCATION, VDSL_ENABLED,
SDO_CS.TRANSFORM(SDO_GEOMETRY(2100,4326,SDO_POINT_TYPE(
CAST(REPLACE(CENTRAL_OFFICES.LAT, '.', ',') AS FLOAT),
CAST(REPLACE(CENTRAL_OFFICES.LNG, '.', ',') AS FLOAT),NULL),NULL,NULL),2100) AS geom
FROM CENTRAL_OFFICES


texts: Are sinks in DEM considered to be errors?

texts:
According to How Sink works in the ArcGIS Desktop 10.3 documentation:

  Naturally occurring sinks in elevation data with a cell size of 10
  meters or larger are rare (Mark 1988), except in glacial or karst
  areas, and generally can be considered errors. (source)

Lakes, ponds, valleys, volcanoes are all naturally occurring "sinks". In fact, any non-monotonic topography has "sinks".

So why are sinks in DEM considered an error?


texts: Create Quick View Orthophotos from aerial images and exterior orientation metadata

texts:
I'm looking for a way to get a quick overview of aerial images taken together with exterior orientation data. Available are horizontal and vertical coordinates of the camera, focal length and also the orientation of the camera at the instant of the exposure. Orientation can be as Roll/Pitch/Yaw or Omega/Phi/Kappa. Also we have available a rough dtm for the mission region.

What I would like is a process to drape the images in correct scale and perspective onto the DTM surface which should result in a somewhat accurate "Quick View Orthophoto". Main purpose are quality checks for example to check for complete coverage of a mission area, sufficient image overlaps, clouds or cloud shadows and so on.

Can anybody point me to a toolset able to do such things?

I'm thinking of maybe using Python in QGIS to process images and metadata into such an project overview.


texts: IApplication in AxMapcontrol (FindExtensionByCLSID in ArcEngine Application)

texts:
I have a small problem. I have begun rewriting code to work with an AxMapControl. Now at first this seemed to work and I was able to get a lot of code converted and working with the AxMapControl, but then I got to the GeoCoding part. I can't seem to figure out what to replace IApplication with to get the FindExtensionByCLSID function like below . . .

pApplication.FindExtensionByCLSID(pUID) as ILocatorExtension

but IApplication isn't available in ArcGIS Engine. How I can get Extension


texts: QGIS raster clip error: Cutline not of polygon type

texts:
I have a raster file and want to clip it, but I'm getting this error:

"ERROR: Cutline not of polygon type"

what does this error mean and how can I resolve it?


texts: How to merge output shapefiles in model builder

texts:
I have a layer with several polygons and a road map. I need to take each polygon and use it to clip the road map, dissolve all of the road records into one, and then merge the polygons into one shapefile. The goal is to have a shapefile with each polygon as one record that contains the total roads within that polygon.

I know the easiest solution would be to dissolve the roads layer before hand and then clip using the polygons but I'm working on a continental scale and my computer keeps running out of memory for that. If I clip the road map and then dissolve everything is fine, the issue is the merge. Using row iterator the model only spits out one shapefile at a time after running a clip and dissolve on it. How do I collect all of those outputs and merge them?

Current Model:
I'm using row iterator to select out each polygon, then use the clip tool to clip the road map, then using the dissolve tool to dissolve all of the roads inside that polygon into one record. Now I'm stuck on how to merge it all back together.

Thanks!


texts: Updating Map Layout Pathways in ArcPy for ArcGIS Pro?

texts:
I have a map which has two insets which looks like this:

The first map frame is in the upper left and the second is all of the background.  Using my exact same map layout as shown here I just want to update the pathway of a layer.

For example, I am just trying to update one layer here, and that is the raster of Alaska shown in the upper left.

I am trying to update the layer like this:

#load in the map
aprx = arcpy.mp.ArcGISProject("my_project")

#get the first map frame
one = aprx.listMaps[0]

#get the layers in that frame
one_layers = one.listLayers()

#list all the pathways
for layer in one_layers:

     print (layer.dataSource)

This returns four layers which are the correct pathways for that map frame.

I want to replace one called "/Users/user/Fall/sen.tif" to "/Users/user/Spring/sen.tif"

I am trying to do this like so:

aprx.updateConnectionProperties("/Users/user/Fall/sen.tif", "/Users/user/Spring/sen.tif")

#save new layout
aprx.saveACopy("/Users/user/new_map.aprx")

I have also tried:

one_layers[2].updateConnectionProperties("/Users/user/Fall/sen.tif", "/Users/user/Spring/sen.tif")

where [2] is the layer location I want to replace but it is the same outcome.

but when I load the new map the layer was not updated.  My desired output is an updated layer with the same symbology as before applied.


texts:
I'm sure there is a really simple explanation to this but haven't found any pointers on this forum or elsewhere. I seem to have lost the 'processing' menu from QGIS and can't figure out how/why!? Any ideas what I might have done to make this happen?

I'm running QGIS 2.8 built against GDAL 1.10 on Ubuntu 14.04.


texts: QGIS label: Capitalize first letter of different words in attribute field

texts: I have a vector layer with an attribute text field with 3 words (word1, word2 and word3). I want to capitalize the first letter of word1 and word3 in label output. For example if a record has  the following words "rio della casa" i want as label "Rio della Casa". I know in QGIS there is the title() statement but it capitalizes the first letter of all the words, so can someone help me to find the right rule for labelling?


texts: Is there an R package for finding the maximum-area rectangle in a non-convex polygon?

texts:
Does anyone know if there is an R package out there with an algorithm which finds the largest rectangle within non-convex polygons? I have Googled, searched the R Sig Geo archives, and looked here on GIS Stackexchange to no avail.

Here is an example of the type of output I am looking for:

This image comes from this blog post which showcases the type of algorithm that would be useful for my project. Unfortunately for me, that algorithm is written in coffeescript and I do not have the time to spend recreating it in R.


texts: Merge overlapping DEMs with irregular boundaries

texts:
I have a 10m DEM covering part of a state, a 20m DEM covering the whole state, and a 90m DEM covering the whole country. How can I merge them together, getting the best quality available in each area?

The 10m and 20m are provided in ADF format, but I've converted them to .tif using ogr2ogr.

Using gdal_merge doesn't work, because the boundaries of each area are irregular. This is the result:

(The northeastern edge of the state is irregular, but there obviously shouldn't be the gap in the middle.)

UPDATE
If it's useful, I also have a shapefile for each 10m DEM, marking its extent.

UPDATE 2

Ok, here's what the whole state looks like with the 20m DEM.

The process I'm going through so far is this:

for f in dtm20m dtm10m_nw dtm10m_e; do
export GDAL_CACHEMAX=1000
echo -n "Re-projecting: "
gdalwarp  -s_srs EPSG:3111 -t_srs EPSG:3785 -r bilinear ./vmelev_${f}/${f} ${f}-3785.tif

echo -n "Generating hill shading: "
gdaldem hillshade -z 5 $f-3785.tif $f-3785-hs.tif
echo and overviews:
gdaladdo -r average $f-3785-hs.tif 2 4 8 16

done


texts:
I'm trying to use a specific marker with my GeoJSON files. They don't show up despite my readings and Firebug.

Can anyone help ?

var map = new OpenLayers.Map('map');
var emprise_style = new OpenLayers.Symbolizer.Point({
    externalGraphic: '../symbol/marker.png'
    });

var empriseStyleMap = new OpenLayers.StyleMap({'default': emprise_style});
var emprise_layer = new OpenLayers.Layer.Vector("GeoJSON", {
        projection: "EPSG:3857",
        strategies: [new OpenLayers.Strategy.Fixed()],
        protocol: new OpenLayers.Protocol.HTTP({
            url: "emprise.geojson",
            format: new OpenLayers.Format.GeoJSON()
            }),
        styleMap: empriseStyleMap
    });
map.addLayer(emprise_layer)


texts: How would you use change the projection within KML file?

texts:
I have created some KML files by converting the projection with both ARC and MapWindow.
I'm looking for ways to automate exporting the data from our application with coordinates in State Plane "ESRI:102272" to something that can be used in Google Earth.

While I'm not a GIS Pro or a Web Developer, I know just enough to export and create maps.

I'm looking to create some kml files to share on our webserver. Can proj4js be used to convert the coordinates within the file?

Any thoughts on this process would be greatly appreciated.


texts: Black edges on overlapped raster images in QGIS

texts:
I'm beginner with QGIS and I have task currently to make some GIS map/data.

I downloaded 3 satellite images from Landsat by USGS (multispectral, georeferenced) and I did overlapping of images using merge tool. I checked with x - no data value and load into canvas when finished. When I used merging tool I've got situation as you can see by the picture.

http://i.imgur.com/xDnaewg.jpg

Images has been overlapped perfectly but with strange black lines on images edges (where overlapping was made). Does anybody can help me regarding this issue? Have in mind that i'm a beginner :)


texts: Modifying the texture/fill of a vector layer in PyQGIS

texts: Currently I am writing a standalone Python script to modify a few vector layers to have a wine fill pattern inside. The geometry I am working with is a polygon. I am aware of how to change outline colors and simple fill colors; but, does anyone happen to know how to add a pattern inside the polygon?


texts: How to recreate unique reference number with QGIS or PostGIS?

texts: I have a shp file with about 500 odd records in it. It has an ID field but I want to add in a new field with a unique reference number. Can I use the field calculator to generate a unique reference? I’m looking for the reference to be alpha numeric, so something like UG_1 or UG1. Alternatively I could load the data into PostGIS and create the field and reference here but I’m still learning the ropes of SQL statements so I’m presuming that this is possible?


texts: Creating grouping for C# ArcObjects UniqueValueRenderer?

texts:
Within a C# procedure I am able to successfully render discrete data values using:

pValRen.AddValue("A11", fieldName, (ISymbol)pExpressWay);
pValRen.AddValue("A15", fieldName, (ISymbol)pExpressWay);

where pValRen is a UniqueValueRenderer() object, and, the strings "A11" and "A15" refer to discrete data values.  But, I need to group "A11" and "A15" together with a single reference value as suggested in the online help.

I have attempted to create a group using:

pValRen.AddReferenceValue("A15", pValRen.ReferenceValue["A11"]);

or

pValRen.AddValue(pValRen.ReferenceValue["A11"], fieldName, (ISymbol)pExpressWay);

but, both of these approaches result in an error.  How does one construct and use a ReferenceValue for the UniqueValueRenderer?


texts: Displaying buffer around feature using population values in ArcGIS for Desktop?

texts:
I am working on a project that is about park accessibility.

I wanted to display a buffer from a park that shows potentially 3 acres/3000 population.  I am using Census tract data and I have a shape of all the parks.

I was able to create a quarter mile buffer and then extract the population data from the shape, but I am asking how do we show how far the reach of a park’s service is for 3000 people?

I am a little familiar with spatial network analyst, but I am not using any streets, I am just relying on the population value in the tracts.

Any suggestions of what I can do?


texts: Finding centrelines from polygons in R sf

texts:
Has anyone found a solution for turning polygons into centrelines using R? It has to be in R as I'm deploying it in parallel on a UNIX cluster. I've got as far as:

library(sf)
centreline <- st_cast(polygon$geometry, "LINESTRING")

But this outputs the entire polygon boundary rather than the centreline.

Here are two solutions in ArcGIS and QGIS.

Creating Centrelines from Road Polygons/Casings using ArcGIS Desktop?

Finding centrelines from polygons in QGIS?

Sample data and code of what I've already tried:

Data

library(sf)
library(dplyr)

# inputs
linepath <- 'river.shp'
pointpath <- 'river_point.shp'
epsg <- 27700

line <- st_read(dsn = linepath) %>% st_transform(epsg)

point <- st_read(dsn = pointpath) %>% st_transform(epsg)

# Buffer point
bufferpoint <- st_buffer(point, dist = 1)  %>%
  st_transform(epsg)

# Cut line with point
# A helper function that erases all of y from x:
st_erase = function(x, y) st_difference(x, st_union(st_combine(y)))

splitline <- st_erase(line, bufferpoint) %>%
  mutate(id = as.factor(row_number()))

buffered <- st_buffer(splitline, dist = 0.5)

dissolved <- st_union(buffered) %>% as('Spatial') %>% st_as_sf()

singlepart <- st_cast(dissolved, 'POLYGON') %>%
  mutate(id = as.factor(row_number()))

centreline <- st_cast(singlepart$geometry, "LINESTRING")


texts: raster-derived KMZ cannot display in google map api v3

texts: I got some raster(tiff)-generated KMZ files, but they cannot be displayed at all in googe map api v3. However, all of my shapefile-generated KMZ files can be displayed properly. Anyone faced the similar situation before?


texts: Seeking C# code to convert Easting/Northing to latitude/longitude in GDA94 / MGA zone 56?

texts:
I am very new to GIS.

I'm pulling data from an ArcGIS Rest API service for a city in Australia to display as pins on a map in another system. The data I'm getting back has all the actual coordinates in x / y coordinates.

I've found that the x,y coordinates I'm getting correspond to the projected coordinate system with the "well known ID" of 28356, which is apparently GDA94 / MGA zone 56 (see https://epsg.io/28356).

This seems correct and covers exactly the area I'm interested in. This web page on the Australian Government's Geoscience Australia website has an online form that confirms I've got the right coordinate system, http://www.ga.gov.au/geodesy/datums/redfearn_grid_to_geo.jsp

How do I convert these x/y coordinates into latitude/longitude in my own code?

Ideally, I'm after a function as simple as:

Coordinate ConvertFromGDA94Zone56ToLatLon(Coordinate source)
{
    // Magic code here
}

I'm writing in C#, if that's any help, but even if you can just give me the maths I can try and work it into code.

[ ...edit... ]

I've found the Proj4 library and worked out that it can do what I need from the command line, but perhaps there's a .NET Core wrapper I'm not aware of...?


texts: URL for querying OSM server, getting xml response, for roads only in large areas

texts:
I've already Googled extensively around this subject, and can't quite get it right.

Can anyone give me example URLs which return XML data consisting of just roads, for example, over areas larger than the 0.5 degrees limited by api.openstreetmap.org? Not had much luck with XAPI.

I'm just looking for ways to get minimal amounts of useful spatial data for regions enclosed by bounding boxes of various sizes.


texts:
I'm attempting to produce interpolated plots of water quality time-series data from ArcSDE.   This is the data format:

OID    WELL_ID  WELL_NAME  SAMPLE_ID  SAMPLE_DAT    PARAM   RESULT

90-------1---------WELL A-----19-24231----5/4/1966----CHLORIDE--------130

91-------1---------Well A-----19-24295------7/7/1966----CHLORIDE--------130

92-------2---------Well B-----19-24478----1/6/1967-----CHLORIDE---------122

93-------3---------Well C-----19-24597----5/5/1967-----CHLORIDE----------115

and onward for thousands of records.  Each well has many records, and the sampling is erratic, some not being sampled for years.   I've tried a few things, but I want to see if anyone else has encountered something similar and how they approached it.


texts: Documenting inability to Start Editing SQLite geodatabase layers in ArcMap?

texts:
I've recently started creating SQLite geodatabases using ArcGIS 10.2 for Desktop with both ST_GEOMETRY and SPATIALITE storage data types and was impressed by how easy this was to do.

I then tried to Start Editing them from ArcMap only to be blocked by an error:

  No editable layers.

and a warning

  The workspace containing this data cannot be edited.

I looked and could not find this limitation mentioned anywhere, and originally mistook it as SQLite being read-only to ArcGIS for Desktop.

Is the non-editability of SQLite within ArcMap edit sessions documented anywhere?


texts: Update the_geom field with polygons from a different table using one common field

texts:
Using CartoDB SQL API I'd like to: Update the_geom field in a table that has no spatial data by selecting the polygons from a table with spatial data and joining based on one common field..

In this particular case, the polygons are property boundaries or parcels and the common field in both tables is the "parcel_id" which is a unique identifier as a string


texts: Select vertices, including multi-part number

texts:
I have an Oracle 18c table with an SDO_GEOEMTRY column:
create table a_test_sdo_geom (shape mdsys.sdo_geometry);

The table has a row that is a multi-part feature.
insert into a_test_sdo_geom (shape)
    values(sdo_util.from_wktgeometry (
    'MULTILINESTRING ((671834.096 4861699.7127, 671836.5099 4861701.9158),'||  --part 1
                     '(671838.2206 4861700.7607, 671842.2311 4861703.3157))'   --part 2
    ));
commit;

Get Vertices:
I can successfully select the individual vertices as rows using the GET_VERTICES() function:
select
    v.id as vertex_num,
    v.x,
    v.y
from
    a_test_sdo_geom a
    ,table(sdo_util.getvertices(a.shape)) v

VERTEX_NUM          X          Y
---------- ---------- ----------
         1 671834.096 4861699.71
         2 671836.510 4861701.92
         3 671838.221 4861700.76
         4 671842.231 4861703.32

Question:
When selecting the individual vertices, is there a way to also select the feature part number?
Example:
VERTEX_NUM          X          Y   PART_NUM
---------- ---------- ---------- ----------
         1 671834.096 4861699.71          1
         2 671836.510 4861701.92          1
         3 671838.221 4861700.76          2
         4 671842.231 4861703.32          2


texts: Is there a OSM data server that returns vector data based on a bounding box?

texts:
I would like to know if exist some OSM data server, not tiled, just data. I know it is possible to import OSM data in your own server, but because you need robust server architecture, it wont be possible for now.

I am looking for a server that returns OSM data using a bounding box. A server that can be used with Mapnik OSM Plugin.


texts: Grass v.Centerline installation with "No Metadata Available" error

texts:
I have Grass 7.8.5 installed and I am trying install the v.centerline extension. I just get a warning with the following, and I see it hasn't installed when going to the Manage Extensions tool.
(Mon Mar 15 14:50:32 2021)
g.extension extension=v.centerline
Downloading precompiled GRASS Addons <v.centerline>...
Fetching <v.centerline> from <http://wingrass.fsv.cvut.cz/grass78/x86_64/addons/grass-   7.8.5/v.centerline.zip> (be patient)...
Updating extensions metadata file...
Updating extension modules metadata file...
WARNING: No metadata available for module 'v.centerline'.
Installation of <v.centerline> successfully finished
(Mon Mar 15 14:50:42 2021) Command finished (9 sec)

I installed this as part of QGIS 3.16. What could be the issue/fix here?


texts: QGIS 3 is there a way to calculate the merging of vector lines

texts: I am using QGIS 3.6.2, I have a vector layer with thousands of vector lines on it. I don't know if this is possible but I want to see if these vector lines merge at all. To be more specific I have trails that I want to see if they merge to one or few points to cross roads. Is there a plugin I can use or a method to find this out?


texts: How do I create rendered vector field in qgis

texts: Perhaps I am misunderstanding the "vector field renderer" and "VectFieldCalc" plugins, but I have no idea how to make a vector field of arrows, with directions and magnitude. I have the data loaded currently with two seperate .vrt layers, and have tried many other ways. The data originally came from two netcdf files. It would be really nice if I could make the layer simply from those two rasters.


texts: How to characterize polygons that are created when you join 2 vector layers?

texts:
I am using QGIS for the first time.
I have a map of an aquifer with 2 different vector layers that split the aquifer differently into several polygons. In the first layer, I used some points(measuring stations) to create Voronoi polygons, in the other layer I got the polygons from vectorizing a Raster that has some land attributes, now I want to know how much area of each of the colorized layer fits into each of the Voronoi polygons.
For example, how much green, blue and light green are in each of the Thiessen polygons below. I need to know for each one, not the total sum, because each Voronoi polygon has some different values that I will have to use later.
Thanks in advance.


texts:
So I'm looking to write a script that can be called by model builder.  What I would like this script to do is to select three records from a feature class that have the highest field values ('AverageGasPrice').  Basically, I want the script to select the three most expensive gas prices. I am not sure where to start, however.

Is there also a way to have the model launch and allow the user to enter which Feature Class they would like to specify as the input FC? If this is not possible, I will also require assistance with how to 'hardcode' the input file using relative path (everything will be stored in one file geodatabase).

Thanks for the help!


texts: Extracting information from Esri Grid Format (ADF) file using R?

texts:
I have an adf file that has a variety of "attribute" or field layers, and I am trying to extract one layer from that adf file, and convert it to a tiff image for further processing in python. Here is an overview of what the data looks like from R:

> r1
class       : RasterLayer
dimensions  : 7272, 10970, 79773840  (nrow, ncol, ncell)
resolution  : 300, 300  (x, y)
extent      : 911643, 4202643, 234254.4, 2415854  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102    +x_0=2500000 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : modef_proj1
values      : 1, 9374  (min, max)
attributes  :
           ID Rowid COUNT    Hio     NIV_Hio
   from:    1     0    31 0.0001   Muy Bajo
   to  : 9374  8491     1 0.9374   Muy Alto

Does anyone know an easy way to extract the Hio attribute layer.m?

When I convert the adf file to a tiff I always only extract the Rowid column. I am comfortable using ArcGIS, R, Python, or GDAL. But, I am not succeeding in converting this file format.


texts: Passing Filter Parameters to GeoServer WFS via URL?

texts:
I'm trying to use PHP's CURL function to fetch geoJSON information from a GeoServer instance. I'm doing this using url variables rather than trying to construct a full xml GetFeature request.

However, I would like to be able to get a subset of the results based on some of the property fields' contents.

So, while I can get all of the results using this url:

http://www.myURL.com/geoserver/namespace/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=layername&outputFormat=json&BBOX=1,2,3,4

can I also limit the results to "Where field LIKE value" or "field = value"?

I've spent an hour trawling through the Geoserver/ECQL documentation and none of it clearly states "append the functions to your URL in this format". As a result, I'm not clear if it's possible to use url variables to perform these tasks, but some people seem to suggest that it is.

Can anyone help with a clear example of a working filter using GeoServer's WFS with URL parameters?


texts: Verifying whether point lies on line or not when both given in GCS?

texts: Googled and tried to look on Esri's forum but couldn't get anything. Rather got things which were verifying if the point is within a polygon or not. But I want to know, if I have starting point and end point of a line on a map with datum as WGS84 and geographic coordinate system, and I have another point in the same system, how can I verify if the point lies on the line or not? I am pretty new in the field of GIS, so maybe it's a very trivial question.


texts: Importing features from MIF-file without their prefiltering

texts:
I wanted to import a MIF layer into QGIS and found out in the attribute table that it was already "prefiltered", meaning that from a total of 7,000 features it was filtered to around 60.
I didn't find an option to deselect the filter. I then attempted to safe the layer as Shape or GeoPackage and reload it, but it didn't work.
Does someone know the reason for it and how to fix it? I have never had these problems with QGIS before.


texts: Duplicate an existing geoserver style using gsconfig

texts:
I am playing around with gsconfig and trying to assign a new style on a layer I have created.
What I want to do is to actually create a new style, which is the same with an existing one.
I searched online and show that they propose to read the sld from a file:

with open(path_to_sld_file) as f:
   cat.create_style('newStyle', f.read(), overwrite=True)

Would it be possible to get the existing style from geoserver and then pass it as a variable in the "create_style" function?

f.read() returns a string so I guess it should be possible if I just pass the variable with the SLD instead of f.read().


texts: What does independently decodable mean in the context of the OSM PBF Format?

texts:
The primitive blocks of the OSM PBF Format are stated to be independently decodable. What should be the extent of this property?

For example, should all the nodes referenced by a way be stored in the same block or does this property only applies to the string tables, lat/lon offsets and granularity attributes?


texts: Adding Visibility Analysis algorithms to standalone QGIS script

texts:
I'm writing a script to compute the intervisibility between two sets of points using the Intervisibility algorithm from the Visibility Analysis plugin which (to my knowledge) is installed by default in QGIS and available in the processing toolbox.

It needs to be standalone because I intend to pipeline and parallelize the process and don't want the interface/overhead of the application. Using the in-application python console to list the available algorithms,

for alg in QgsApplication.processingRegistry().algorithms():
    print(f"{alg.id()} --> {alg.displayName()}")

I get the following options at the end of the list.

...
qgis:vectorlayerscatterplot --> Vector layer scatterplot
qgis:voronoipolygons --> Voronoi polygons
visibility:Intervisibility --> Intervisibility network
visibility:Viewshed --> Viewshed
visibility:create_viewpoints --> Create viewpoints

I understand that to use these algorithms in the standalone script I need to add the algorithms provider. I do so using the following code.

import os
import sys

sys.path.append(os.path.abspath('/home/USER/.local/share/QGIS/QGIS3/profiles/default/python/plugins/'))

from qgis.core import QgsApplication

from ViewshedAnalysis.visibility_provider import VisibilityProvider
from qgis.analysis import QgsNativeAlgorithms
qgs = QgsApplication([], False)
qgs.initQgis()

import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(VisibilityProvider())
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

However, the same code to list the algorithms now ends with,

...
qgis:vectorlayerscatterplot --> Vector layer scatterplot
qgis:voronoipolygons --> Voronoi polygons
: -->
: -->
: -->

And, calls to visibility:Intervisibility fail with Algorithm not found error. What am I missing?


texts: PyQGIS: Programmatically Adding Point Features via QgsMapToolEmitPoint

texts:
I've created a map tool for QGIS 2.0 Python plugin that emits a point with the click of the mouse. I'm trying to create an edit tool, adding a point feature at the emitted point while selecting another layer's feature at the same time.
The map tool was created under __init__ like so:
def __init__(self, iface):
    …
    self.canvas = self.iface.mapCanvas()
    self.clickTool = QgsMapToolEmitPoint(self.canvas)
    …

Connected the map tool under initGui:
def initGui(self):
    …
    self.clickTool.canvasClicked.connect(self.addAndSelect)
    …

The function addAndSelect() should add a point and select features. I think my problem lies here, but I'm not sure what it is. Currently, the elements in my other layer are selected and the feature form opens, however, no feature is added to the attribute table and all info filled in on the form is lost:
def addAndSelect(self, point, button):
    layer_caps = self.myLayer.dataProvider().capabilities()
    self.feat = QgsFeature()
    pntGeom = QgsGeometry.fromPoint(point)
    self.feat.setGeometry(pntGeom)
    self.myLayer.addFeatures([self.feat], False)

    pntBuff = pntGeom.buffer((self.canvas.mapUnitsPerPixel() * 2), 0)
    rect = pntBuff.boundingBox()
    Nature_idx = self.resLayer.fieldNameIndex('Nature')
    Diam_idx = self.resLayer.fieldNameIndex('Diametre')
    self.selectList = []
    self.natureList = []
    self.diametreList = []
    iter = self.otherLayer.getFeatures()
    for feature in iter:
        if feature.geometry().intersects(rect):
            self.selectList.append(feature.id())
            self.natureList.append(feature.attributes()[Nature_idx])
            self.diametreList.append(feature.attributes()[Diam_idx])
    self.otherLayer.setSelectedFeatures(self.selectList)

    self.iface.openFeatureForm(self.myLayer, self.feat, False)

    self.myLayer.updateExtents()
    self.myLayer.reload()
    self.canvas.refresh()

    self.myLayer.editingStopped.connect(self.unsetMyMapTool)

The map tool is set under the run method:
def run(self):
    …
    if result == 1:
        self.canvas.setMapTool(self.clickTool)
        self.myLayer.startEditing()

The map tool is unset with unsetMyMapTool method:
def unsetMyMapTool(self):
    self.canvas.unsetMapTool(self.clickTool)


texts: OpenLayers 6 - reprojected GeoJSON data not displaying on map

texts:
I'm using latest OpenLayers (6.1.1) and proj4js 2.6.0 to display a multipoint geometry.

I first define my input data project with proj4 then create an open layer vector source as follow:

proj4.defs("EPSG:2154","+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
ol.proj.proj4.register(proj4)

var profile = '{ "type": "Feature", "geometry" : { "type": "MultiPoint", "coordinates": [ [ 308007.223185, 6652665.6637, -0.167366668582 ], [ 308008.22140099998796, 6652665.6039899997413, -0.131412133574 ] ] } }'

var vectorSource = new ol.source.Vector({
    features: (new ol.format.GeoJSON()).readFeatures(profile, {
        dataProjection: 'EPSG:4326',
        featureProjection: 'EPSG:2154',
    }),
});

var vectorLayer = new ol.layer.Vector({
    source: vectorSource,
});

var map = new ol.Map({
  target: 'map',
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM(),
      name: 'source-osm',
      visible: true,
      opacity: 0.7,
    }),
    vectorLayer,
  ],
  view: new ol.View({
    center: ol.proj.fromLonLat([-1.95, 47.2]),
    zoom: 10,
  })
});

Problem is the data doesn't appear on map. If I have a look at the vector source object in the console, it tells me that projection_: null.

If I do not specify a reprojection, points are displayed (but at an incorrect position obviously).

I checked using epsg.io that my points coordinates are correct.

Any idea?


texts: SQL Server Spatial Enable Spatial column and Analysis

texts:
I have:

a table called ABC with columns (ID, County_Name, Lat, Long) in SQL
server pre-filled with the values under its respective columns,
except for County_Name column which is currently null.
an ArcGIS SDE feature class - called US_County with columns (id,
County_Name, Shape).

I am trying is call the SQL Server Spatial native STintersects() to get the county name from the US_County feature into my other non-spatial ABC table.  But I can't seem to get it working.

First I add in the Shape column to the ABC table like so:

ALTER TABLE ABC
ADD Shape AS geometry::Point(Long,Lat,4326);

Then I do the STIntersects() like so:

UPDATE ABC
SET County_Name = a.County_Name
FROM ABC c
JOIN US_County a
ON c.Shape.STIntersects(a.shape) =1;

Was wondering if the above SQL statements is correct?


texts:
I am Esri-trained, trying to move into QGIS using Python.  Simply trying to add a vector layer in QGIS returns the below error:

TypeError: QgisInterface.addVectorLayer(QString, QString, QString): first argument of unbound method must have type 'QgisInterface'

Here's my code:

layer = qgis.core.QgsVectorLayer("C:\mypath\mypath\GTM_adm1.shp", "GTM_adm1", "ogr")
qgis.gui.QgisInterface.addVectorLayer(layer, "GTM_adm1", "ogr")

Not sure what I'm missing here but the Vector layer object is in the registry.


texts: error message: "Cannot split, erase or overlay with points, collection, Text or layout frame objects" in Mapbasic

texts:
In order to automate procedures using MapBasic, I wrote the code for overlapping two objects -- one a polygon object and the other a point object -- using Objects Intersect  statement written below.
It gives me the above mentioned error.

But actually, when I manually do the intersect via Mapinfo application (via Table\Update) it gives me a result which is a modified target table (BB.TAB).

'Open the point table
Open Table "C:\Temp\AA.TAB" Interactive

'Setting the Target
Select * from BB
Set Target On

Select * from AA

Objects Intersect Into Target Data CC=CC,DD=DD

How can I rectify Mapbasic error ?


texts: Get coordinates of all pixels in a raster

texts:
I have a georeferenced raster and I need to extract a list with the coordinates of pixel in reference system (EPSG:32631) and the coordinate of pixel in raster reference.

I want a list like this (just an example):

x_raster,y_raster,Nord,East
0,0,700000,4780000
0,1,700000,4780010
...
...

How can I do using Python and GDAL?


texts: Mapbasic set mapwindow for hotlink code to run

texts:
I need to run a hotlink process on a dataset but I keep getting a "Could not convert data"

Map From DSC_DCDB
Browse * from DSC_DCDB
Set Map Window DSC_DCDB  Layer 1 Activate Using PATH_Survey+HOTLINK_Survey On Labels Objects Relative Path Off Enable On ,Using PATH_JUMPUP+HOTLINK_JUMPUP On Labels Objects Relative Path Off Enable On ,Using PATH_ASCON+Hotlink_ASCON On Labels Objects Relative Path Off Enable On

What should I do? I want to run this as part of a larger mapbasic process.


texts:
Is there an open source software that I can use to collect data with the Juno 3D GPS?

I'm hesitant to buy a Terrasync or ArcPad license.


texts:
I want to transform lon/lat point coordinate to Albers X/Y cordinate. So I defined a function as follows:
This code modified base on section of Coordinate conversion  of answer

import numpy as np
from osgeo import gdal,osr

inprojF=r'\WGS 1984.prj'          #needed to tranform projection
outprojF=r'\WGS_1984_Albers.prj'  #wanted projection

def convertXY(x,y):  # function to convert coordinates
    xy_source=np.array([x,y])
    inprojds=open(inprojF,'r').read()
    outprojds=open(outprojF,'r').read()

    inproj=osr.SpatialReference()
    inproj.ImportFromWkt(inprojds)
    outproj=osr.SpatialReference()
    outproj.ImportFromWkt(outprojds)

    ct = osr.CoordinateTransformation(inproj, outproj)
    xy_target = np.array(ct.TransformPoints(xy_source))
    xx = xy_target[0]
    yy = xy_target[1]
    return xx, yy

When I type:

convertXY(110.0,40.0)

Hint error:

in method 'CoordinateTransformation_TransformPoints', argument 1 of type    'OSRCoordinateTransformationShadow *'

Debugging the code I found that the error is in ct.TransformPoints(xy_source)

How to improve above code?


texts: programmatically create and publish geotiff layer in geoserver

texts:
We need to automatize process of creating datastore and publishing geotiff layer. Also, before publishing we need to change Declared SRS and to save Native Bounding Box to our database.

Geotiff file will be uploaded to the file system by our customer.

Yes, I`ve seen Create a Layer in GeoServer using REST, also I see that there is publishGeoTIFF(...) methods in Geoserver-manager in GeoServerRESTPublisher.java.

Another solution that I see - is to make all the operations with bot written with selenium or htmlUnit.

Could you please suggest us, what approach is more elegant and suitable for us. Is there any actions in our workflow that we will not be able to do with Geoserver-manager.

P.S.: we also have seen at gsrcj project, but it can`t create programmatically geotiff layer.

Thank you in advance,
Dmitry.


texts: Using polygon information to compare/contrast with point information on same map in MapInfo?

texts:
For my project at work, I am trying to compare and contrast property information with points I have mapped (from my workplace)

Example:
I have points with information on water samples, the results of their chemical testing, well depth and date sampled.

I have tax parcel polygons for each county (U.S.) and that table can tell me whether that property has a private or public well

I'm interested in finding what points are on properties that have a private or public well, or no well.

I was thinking of trying to make a thematic map for each, but it gets confusing fast. I was thinking of applying X color for private wells, Y color for public supply, and Z for none or other.
After that, I was going to apply a range of low to high levels of toxicity for the chemical our workplace is investigating, that has 4 different colors.

Is there a less complicated way of comparing/contrasting polygon information with point information in two separate tables that can't be joined?

I am new to MapInfo, and have not had an easy time finding guidance on how to interact with the software. This is my first GIS software program that I've jumped into using, with just a few hours of training and a lot of trial and error.


texts: Export topology checker results as table/shapefile in QGIS

texts:
I have used the plugin "topology checker" and I got table with all the topology errors :

I would like to save this table as new layer/ csv or any way  I can export it.
I saw there were disccussion about it few years ago in github but couldn't find any feature to do that.
Is there any way to export this table?


texts: What is the default projection in Leaflet?

texts:
I've read an API documentation for Leaflet on its site, but I can't realize: what is a default projection to view maps? It was written about spherical and elliptical Mercator and Equirectangular - so which of this three? Or does it depend on tile server I use, or it doesn't depend since leaflet reproject tiles instantly on-the-fly before adding on the screen?

Actually i can't understand what projection i need to save geoJSON files in QGIS. To display correctly on the leaflet map, json coords should be lat-long, not x-y. But if i try to save it in Equirectangular - coords will be x-y..?


texts:
I am currently working with spatial data in R using sp When I used WGS84-referenced coordinates, I used to specify it by setting the proj4string to:

+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84"

But I don't really understand how this works internally. Now I want to use data that is referenced to ETRS89 datum. How should I set my proj4string? Should I use a reprojection?


texts: Spatial join without temporary layer in QGIS

texts:
I can't find how to perform a spatial join in QGIS without a temporary layer created.
What I want is to update directly the table in the database (Oracle in my case), for a specific column.

Any idea how to do so ?

Note : For this test I have a point layer in which I want to put a field value from a polygon layer which carry this value. So intersect or within is enough to test.


texts: Checking if borders of layers match using ArcGIS Desktop?

texts:
I have two shapefiles which are made out of polygons within one city.

Districts
Parishes

I need to check if their borders within the city match- if the parishes are within districts or if they somewhere extend beyond one district to another. The most outter borders match the city border so it's ok.

Is there any way to do so using any ArcGIS Desktop 10.3 or 10.5 toolbox (I have both versions available to me) besides going through it manually?


texts: Using travel costs ($) as impedance in ArcGIS Network Analyst?

texts:
I am building a multimodal network with 4 modes (Walking, bike, bus and metro) using ArcGIS Desktop 10.5.

Now i would like to use costs in $ as an impedance for the best (thus cheapest) route.
However, I have a metro that uses a fare/km but a bus system with a flat rate.
This means u pay for example 2$ to get onto the bus system and from then on can go wherever you want to go within the system.
Now how do i model this in a correct way, without charging the fee multiple times for each transfer?


texts: Problems with latitude/longitude grid in QGIS 2

texts:
I try to produce a thematic map in QGIS, like this one done in MapInfo:

I want to do this in QGIS. I work in projection "Albers-Siberia" projection ( PROJ.4: +proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105 +x_0=18500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs). The latitude/longitide grid I made with help of this article. I have got the following result:

I have a problem:

As you can see grid lines don't intersect in the North Pole. How
    can I solve this problem?


texts:
I am using ArcGIS Pro and cannot find a way to undo field calculations.

With ArcGIS 10.x for Desktop, field calculations could be undone when an edit session was active, but Pro seems to be in a perpetual edit session that does not include the ability to undo field calculations.

Is it possible to undo a field calculation in ArcGIS Pro?


texts:
How can I create "buffers" or vector from values from a layer of points, where the "buffers" are "tied" to the streets?

Example: A bus station, I would like to create vectors within 100 m range with streets, from the point, regardless of the direction, following the road system.

I need to make a linear "buffer", following the system road, to identify the distance the bus station will cover on every street.


texts: How to make polygon from cells of a raster?

texts:
I have 100 rasters, each raster has the same cell size and extent and i want to make a polygon of each cell.

Can anyone help me doing this?


texts: How to apply compression (e.g. LZW) when exporting a map from QGIS composer to TIF?

texts:
QGIS > Composer > Export as image > *.TIF

I need to apply lossless compression (e.g. LZW or ZIP), when exporting a map from QGIS (v2.18) composer to TIF image. This settings are missing in the export dialog. Where/How can I set it up?
(I know that I can export and re-save with external software like Imagemagick, or export to PNG and re-save in external software... this is just an unnecessary step).


texts: Creating model to copy feature classes from ArcSDE to FIle Geodatabase?

texts: I am trying to create a model in ArcGIS that copies the feature classes which are stored in SDE geodatabase “SQL” to file geodatabase “in order to sync my SDE data to FGDB“, here what is the best practice to do this, I mean I am trying to use copy features tool, but It seems not efficient because of the following reason which is that I have 10 feature classes stored in my sde database so I need to use copy feature 10 times with 10 parameters and if one of my F.C  was renamed then I should edit the model again , If there is any tool or method that copies the whole sde to FGDB


texts: Edit a raster (DEM) to create a constant slope angle

texts:
I'm trying to edit a raster (DEM) that has square sides (90 degrees) to create a surface with a constant slope angle of 18 degrees.

I'm working in QGIS and have not been able to find any help for this issue.


texts: Failed to create dynamic vector tiles using PostGIS on existing table

texts:
I am unable to run PostGIS ST_AsMVTGeom. I read Serving Dynamic Vector Tiles from PostGIS by Paul Ramsey and the PostGIS documentation.
The code is below:
SELECT ST_AsMVT(mvtgeom.*) FROM (
   SELECT ST_AsMVTGeom(point_geom) AS geom
   FROM ancillary_data
) mvtgeom

The error message:
ERROR:  function st_asmvtgeom(point) does not exist
LINE 1: SELECT ST_AsMVTGeom(point_geom) AS geom
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
SQL state: 42883
Character: 8

Surprisingly the code below works.
SELECT ST_AsText(ST_AsMVTGeom(
ST_GeomFromText('POLYGON ((0 0, 10 0, 10 5, 0 -5, 0 0))'),
ST_MakeBox2D(ST_Point(0, 0), ST_Point(4096, 4096)),
4096, 0, false));

Edit:
I realized that I was using PostgreSQL point column instead of PostGIS geometry.But the problem is still unsolved.
I am using PostGIS 3.0.1 in PostgreSQL 9.6.


texts: What possible explanations are there for a dramatic change in vegetation appearance in Sentinel 2 data?

texts:
I'm looking at land cover in an area of South East Asia. I pulled down Sentinel 2 imagery for Feb 2017, and used the QGIS SCP plugin to complete the land cover analysis, as per this tutorial: https://fromgistors.blogspot.com/2016/09/basic-tutorial-2.html. Wanting to better understand which areas of lesser growth were plant growth recovering from cultivation, and what were otherwise sparse forest, I took a look at some older satellite data, only to find it looked very different. I've included examples below using using Sentinel bands 4,8 and 12, all of which were processed in identical ways.

If this is indeed just a difference in the visual output of SCP/QGIS, I would love to know how to get them more comparable.

Jan 2016

Jan 2017

Feb 2017


texts: Mapbox Supercluster on the server side, how?

texts:
My tile server pipeline is like this: GeoJSON -> MBTiles -> PBF. The tools that I used are tippercanoe and mbutil.

I'm not quite familiar with the vt-pbf, geojson-vt tools. Could anybody help explain the relationship of these tools compare to the mbutil,tippercanoe? I want to use the supercluster on the server side. How to achieve that? What to change on the mapboxGL-JS side?

Reference:
https://www.mapbox.com/blog/supercluster/


texts:
I am trying to add rotation to a NetCDF File

One Idea was to add a additional parameter to my current grid mapping that describes the rotation point and the rotation:

PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",
  GEOGCS["GCS_WGS_1984",
    DATUM["D_unknown",
      SPHEROID["WGS84", 6378137.0, 298.257223563]],
    PRIMEM["Greenwich", 0.0],
    UNIT["degree", 0.017453292519943295],
    AXIS["Longitude", EAST],
    AXIS["Latitude", NORTH]],
  PROJECTION["Oblique_Mercator"],
  PARAMETER["longitude_of_center", 7.0],
  PARAMETER["latitude_of_center", 51.4],
  PARAMETER["azimuth", -10.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  PARAMETER["rectified_grid_angle", -10.0],
  UNIT["m", 1.0],
  AXIS["x", EAST],
  AXIS["y", NORTH],
  AUTHORITY["EPSG","8011112"]]

Another try with Affine Rotation:

FITTED_CS["BPAF",
  PARAM_MT["Affine",
    PARAMETER["num_row", 3],
    PARAMETER["num_col", 3],
    PARAMETER["elt_0_0", 0.16867854068466956],
    PARAMETER["elt_0_1", -0.2699420929659323],
    PARAMETER["elt_0_2", 918079.626281209],
    PARAMETER["elt_1_0", 0.2699420929659323],
    PARAMETER["elt_1_1", 0.16867854068466956],
    PARAMETER["elt_1_2", 6445039.217828758]],
  PROJCS["WGS84 / Google Mercator",
    GEOGCS["WGS 84",
      DATUM["World Geodetic System 1984",
        SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
      PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
      UNIT["degree", 0.017453292519943295],
      AXIS["Longitude", EAST],
      AXIS["Latitude", NORTH],
      AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["semi_minor", 6378137.0],
    PARAMETER["latitude_of_origin", 0.0],
    PARAMETER["central_meridian", 0.0],
    PARAMETER["scale_factor", 1.0],
    PARAMETER["false_easting", 0.0],
    PARAMETER["false_northing", 0.0],
    UNIT["m", 1.0],
    AXIS["x", EAST],
    AXIS["y", NORTH],
    AUTHORITY["EPSG","900913"]],
  AUTHORITY["EPSG","8011113"]]

Another idea was to use two dimensional coordinates that are pointing directly to the georeferenced coordinates:

netcdf test {
dimensions:
        time = UNLIMITED ; // (3 currently)
        y = 10 ;
        x = 27 ;
        z = 2 ;
variables:
        double time(time) ;
                time:long_name = "time" ;
                time:units = "hours since 2018-06-23T05:00:00Z" ;
                time:calendar = "standard" ;
                time:axis = "T" ;
        double lat(y, x) ;
                lat:long_name = "latitude" ;
                lat:units = "degrees_north" ;
        double lon(y, x) ;
                lon:long_name = "longitude" ;
                lon:units = "degrees_east" ;
        double TSurf(time, y, x) ;
                TSurf:_FillValue = -999. ;
                TSurf:valid_min = -50.f ;
                TSurf:valid_max = 150.f ;
                TSurf:long_name = "Surface Temperature" ;
                TSurf:units = "degree Celsius" ;

// global attributes:
                :Conventions = "CF-1.7" ;
}

All three methods don't seem to work properly in QGis or Geoserver

Anyone has already some practice with NetCDF raster data that is not north up oriented?


texts: Using max as resampling method with gdal_translate?

texts:
My Issue

According to the docs it appears that gdalwarp has many more resampling methods than gdal_translate. Can anybody explain why the resampling methods can't be used with both functions? Alternatively can you explain how to achieve the same result?

I have a big vrt that I'm trying to crop into a raster with a lower resolution covering a smaller geographic region. gdal_translate has the '-projwin' option which allows me to select smaller region without reading the entire dataset.

I see that gdalwarp has the '-te' option which would give me a similar result, except that it looks like it tries to read the entire dataset into memory prior to cropping it to the extent that I specify. (ERROR 1: Integer overflow).

I would be happy using either method as long as the result that I get comes in a somewhat performant manner for a small region.

Things I've tried

gdal_translate -r max … - GDAL_RASTERIO_RESAMPLING = max not supported
gdal_warp -te xmin ymin xmax ymax … - Integer overflow .
A combination of translate to crop the image without downsampling, then gdal_warp to downsample. While this works the two operations take too long for my use case.

Exact commands tried

gdalwarp -ts 10 10 -te -122.4396 47.5189 -122.1892 47.6830 dem.vrt test_wc.tif -overwrite --config GDAL_CACHEMAX 100 -wm 100 -multi -wo NUM_THREADS=ALL_CPUS
gdalwarp -ts 10 10 -te -122.4396 47.5189 -122.1892 47.6830 dem.vrt test_wc.tif -overwrite --config GDAL_CACHEMAX 100 -wm 100 -multi -wo NUM_THREADS=ALL_CPUS

What can I do to use the max resampling method, while cropping a large dataset?

Gdal Info

Driver: VRT/Virtual Raster
Files: dem.vrt
       LongFileList…
Size is 1296003, 421203
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000416666666666,61.000416666666666)
Pixel Size = (0.000277777777778,-0.000277777777778)
Corner Coordinates:
Upper Left  (-180.0004167,  61.0004167) (180d 0' 1.50"W, 61d 0' 1.50"N)
Lower Left  (-180.0004167, -56.0004167) (180d 0' 1.50"W, 56d 0' 1.50"S)
Upper Right ( 180.0004167,  61.0004167) (180d 0' 1.50"E, 61d 0' 1.50"N)
Lower Right ( 180.0004167, -56.0004167) (180d 0' 1.50"E, 56d 0' 1.50"S)
Center      (  -0.0000000,   2.5000000) (  0d 0' 0.00"W,  2d30' 0.00"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  NoData Value=-9999999


texts: Inverted shading with polygons in QGIS 2.8

texts:
Is there a way to modify the style of a poygon .shp file so that everything except the polygon is rendered as shaded?

EDIT: see these images if you have trouble understanding. I am going for the effect in the second image.


texts: Can QGIS resample Rasters during display?

texts:
When I was using ArcGIS a user could resample a raster during the display (I think via layer properties under the display tab). This would smooth the appearance of the raster so it wouldn't appear pixelated when the user zoomed in too far, I found this useful when using hillshades.

I haven't found a way to do this in QGIS, but can it be done? If not is there a GDAL tool I can use to smooth the appearance of my hillshades so they don't look pixelated when they are getting near the end of their optimum scale.

Thanks
Ando


texts:
I am a new in PyQGIS, I want to create a Point layer from a raster layer. In table of point layer created, there is 3 fields by name of "X_coord", "Y_coord" and "Azimuth", each point in Point layer shows x ("X_coord" field), y ("Y_coord") center coordinate of each pixel from raster layer and "Azimuth" field shows azimuth from another single point layer to each center pixels. When I run this script for raster layer with 67*52 pixels in 30 meter resolution time for running is 5 seconds, for 133*105 pixels time is 15 seconds and for 291*220 pixels time is 22 minutes!
I think this is unusual.How can I decrease time for run this script for large data?

My script is:

from qgis.core import *
from osgeo import gdal
from PyQt4.QtCore import *
import time

startTime = time.clock()

pnt = QgsVectorLayer(r'E:/Sample_Dataset_Lorestan/Point.shp', 'POINT','ogr')
for feat in pnt.getFeatures():
   geom = feat.geometry()
   pntxy=geom.asPoint()

# Open tif file
ds = gdal.Open(r'E:\Sample_Dataset_Lorestan\subset_a1.tif')

numpy_array = ds.ReadAsArray()
# GDAL affine transform parameters

geotransform = ds.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
rotationX = geotransform[2]
rotationY = geotransform[4]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]

cols=ds.RasterXSize
rows=ds.RasterYSize

def pixel2coord(x, y):
   xp =  (pixelWidth * x) + originX + (pixelWidth/2)
   yp =  (pixelHeight * y) + originY + (pixelHeight /2)
   return(xp, yp)

#Create temporary vector layer and add to map
vl = QgsVectorLayer("Point?crs=", "temporary_points", "memory")
pr = vl.dataProvider()
QgsMapLayerRegistry.instance().addMapLayer(vl)
# Add points:
pr = vl.dataProvider()
# Enter editing mode
vl.startEditing()
# add fields
pr.addAttributes( [ QgsField("X_coord", QVariant.Int), QgsField("Y_coord", QVariant.Int), QgsField("Azimuth", QVariant.Int) ] )

# add a feature
fet = QgsFeature()

# get columns and rows of your image from gdalinfo

for row in range(0,rows):
   for col in range(0,cols):
    rspnt = pixel2coord(col,row)
    rspnt2 = QgsPoint(rspnt[0], rspnt[1])
    fet.setGeometry(QgsGeometry.fromPoint(rspnt2))
    fet.initAttributes(3)
    fet.setAttribute(0, rspnt[0] )
    fet.setAttribute(1, rspnt[1])

    az = pntxy.azimuth(rspnt2)
    if az < 0.0:
         azim = 360 + az
    else:
        azim = az

    fet.setAttribute(2, azim)
    pr.addFeatures( [fet] )

# Commit changes
vl.commitChanges()
print "completed within:", time.clock() - startTime, "seconds", "\n"


texts: deCarta Android API geocode throws Protocol not found exception

texts:
I am using the deCarta Android API to find the GPS Lat/Lng of a street address.  However, whenever I try to use the geocode method of the Geocoder class it always throws the exception:

Protocol not found

I am passing in a StructuredAddress object populated with the street address for Google:

address.toString():

1600 Amphitheatre Pkwy
Mountainview, CA 94043

ArrayList<GeocodeResponse> results = Geocoder.geocode(address,false);

Please advise.


texts: field map string values in one field to domain field

texts:
Is there a method in ArcGIS 10.3.1 to calculate coded value domains in a domain field based on a string field source? Please see image.

The solution appears to be to use the CalculateField GP in ArcCatalog as opposed to calculating within an ArcMap session. While this method works, is it possible to display the coded values in the field without the descriptions?


texts: Contouring in ArcGIS from raster file yielded rectangular contours

texts:
Why does Contouring in ArcGIS using Contour tool in ArcGIS from raster file yielded rectangular contours?
The attached image of contour is shown below

The contours shows a rectangular gridded contours which doesn't look like contours
I have also tried the same with QGIS and get the similar results.


texts: WMTS: declare "GetTile" operation for statically hosted tiles?

texts:
Supposedly, one can turn a collection of static tiles, hosted according to the ZXY schema also used by OpenStreetMap tiles and Google Maps tiles, into a standard-conformant OGC WMTS by accompanying it by an (also static) WMTSCapabilities.xml file for satisfying the GetCapabilities WMTS operation. (See also question GET encoding for static WMTSCapabilities.xml.)

Amongst others, that XML document should describe how tiles can be obtained with parametrized URLs. If the tiles are only available at URLs like (say) https://example.com/map/1/2/3.png (this one being for zoomlevel 1, tile column 2, tile row 3) but not at any URLs like

https://example.com/map.cgi?service=WMTS&request=GetTile&version=1.0.0&layer=example&style=default&format=image/png&TileMatrixSet=WholeWorld_CRS_84&TileMatrix=1&TileRow=2&TileCol=3

then should the ows:OperationsMetadata have an entry for the GetTile operation like

<ows:OperationsMetadata>
   <!-- ... -->
   <ows:Operation name="GetTile">
      <ows:DCP>
         <ows:HTTP>
            <ows:Get xlink:href="https://tile.osm.ch/switzerland/">
               <ows:Constraint name="GetEncoding">
                  <ows:AllowedValues>
                     <ows:Value>KVP</ows:Value>
                  </ows:AllowedValues>
               </ows:Constraint>
            </ows:Get>
         </ows:HTTP>
      </ows:DCP>
   </ows:Operation>
</ows:OperationsMetadata>

or should only the ResourceURL elements of the respective Layer elements specify how to obtain tiles? The latter looks like

<Contents>
   <Layer>
      <ows:Title>My Map</ows:Title>
      <ows:Identifier>mymap</ows:Identifier>
      <ows:WGS84BoundingBox crs="urn:ogc:def:crs:OGC:2:84" dimensions="2">
         <ows:LowerCorner>-180 -85.051128779807</ows:LowerCorner>
         <ows:UpperCorner>180 85.051128779807</ows:UpperCorner>
      </ows:WGS84BoundingBox>
      <Style isDefault="true">
         <ows:Identifier>default</ows:Identifier>
      </Style>
      <Format>image/png</Format>
      <TileMatrixSetLink>
         <TileMatrixSet>SOSM standard</TileMatrixSet>
      </TileMatrixSetLink>
      <ResourceURL format="image/png" resourceType="tile" template="https://example.com/map/{TileMatrix}/{TileCol}/{TileRow}.png" />
   </Layer>
   <!-- ... -->
</Contents>

, so template="https://example.com/map/{TileMatrix}/{TileCol}/{TileRow}.png" specifies where to find the tile, anyway.


texts: How to display latest sample value for a point in QGIS?

texts:
We have a dataset consisting of points (x,y). At each of these points, analytical data has been collected at different dates. The dates are not the same for all points. An example of the dataset is shown below:

Point; x; y; date; analytical data
B1; x1; y1; 1-1-17; 10
B1; x1; y1; 1-2-17; 100
B1; x1; y1; 1-3-17; 5
B2; x2; y2; 1-1-17; 1
B2; x2; y2; 1-3-17; 10
B2; x2; y2; 1-4-17; 100
B3; x3; y3; 1-2-17; 1000

The question is now: In QGIS, how can we plot the 'analytical data' from the latest sampling date? So, one value for each point.

We assume this can be done by using 'an expression', we just cannot figure out how?


texts: How to add TIN surface Z values into Z field of line vertices in ArcGIS Desktop?

texts: I want to transfer Z values of TIN surface into Z field of line vertices. This line feature class is covered by this TIN surface. I used interpolation shape tool but this tool is disrupting line during interpolating process (Linear Method). Is there any other way to write Z values of this TIN surface into the z field of line vertices?


texts: PostgreSQL/PostGIS - group by series of values

texts:
I need to sum up linear objects' length and to group them by intervals of several years. I have one table storing my objects like this :

- gid serial NOT NULL,
- year INTEGER,
- the_geom geometry(MULTILINESTRING) ;

I need a result like this :

period          | length
----------------+-----------
 2005 - 2014    | 18.6
 1995 - 2004    | 16.1
 1985 - 1994    | 7.6
 1975 - 1984    | 19.0
 1965 - 1974    | 28.2
 1945 - 1964    | 10.2
 before 1945    | 0.1

I can't find out on the Web how to perform this, except by producing each line with a different query and to use UNION ALL to merge them together, which is not so good...

Many thanks for help !

wiltomap


texts: Creating view on ArcCatalog using tables from two different databases residing on different servers?

texts:
I have some features (Roads, Drainages and Parcels - Line, Point and Polygon) data. The spatial attributes of these features are captured by drawing them into arcmap and then stored feature class in a feature dataset  in a geodatabase that is on MSSQL2012 database on Server A. The non spatial attributes of these features are existing on another DEVSQL2012 on Server B.

I am able to connect to the different databases using ArcCatalog Database connection and so able to view both the spatial and non-spatial data separately.

I am now trying to create a database View from my ArcCatalog to join the spatial table on Server A  to  the non-spatial table on Server B.  But I keep getting various errors.

Please does anyone has an idea how I can get this resolved.

SELECT [DEVSQL2012].TVR_CIPS.dbo.Projects.Id,
    TVR_CIPSdboProjectsName
    TVR_CIPS.dbo.Projects.Name,
    TVR_CIPS.dbo.Projects.ProjectTypeId,
    TVR_CIPS.dbo.Projects.ProjectManagerId,
    TVR_CIPS.dbo.Projects.IsFunded,
    TVR_CIPS.dbo.Projects.Scope,
    TVR_CIPS.dbo.Projects.FundSource,
    TVR_CIPS.dbo.Projects.TotalEstimatedCost,
    TVR_CIPS.dbo.Projects.Consultant,
    TVR_CIPS.dbo.Projects.Contractor,
    TVR_CIPS.dbo.Projects.LocationLimits,
    TVR_CIPS.dbo.Projects.GISId,
    TVR_CIPS.dbo.Projects.StatusId,
    TVR_CIPS.dbo.Projects.PermalinkMapURL,
    TVR_CIPS.dbo.Projects.CIPDocURL
FROM TVR_CIPS.dbo.Projects INNER JOIN
     [Frussmolt\TVR3PRD02].PublicWorks.TC_USER.CIP_Point_Only
ON
    [DEVSQL2012].TVR_CIPS.dbo.Id = [Frussmolt\TVR3PRD02].PublicWorks.TC_USER.CIP_Point_Only.Project_ID

Here is the error I'm getting

I don't understand why it is seeing that column as invalid. The 'IsFunded' column is supposed to be a Yes or No (0 and 1) field, though the field datatype for this column is 'Short integer', Could this be a problem?


texts:
Does anyone know if there are any sample questions available for the ESRI Web Application Developer Certification exam?

ESRI has released 3 sample questions for this test:

http://downloads2.esri.com/campus/downloads/certification/samplequestions/EWDA10_Sample_Questions.pdf


texts: color gradient (or color ramp) in Arcgis and Qgis

texts:
I am currently working on my thesis and I would like to ask you:
does anybody know the mathematical formula that calculates the color gradient in Arcgis and Qgis? In other words, which mathematical equation is used to represent the color transition from one color to another in a buffer? And also: is there a different equation that calculates the color classification in other shapes, such as square, circle etc? A good example can be found in How to Calculate Color Gradient (www.ehow.com), but this is rather general. I would prefer something that concerns Qgis and Arcgis.

EDIT: 20/01/2013

The arrow doesn't mean anything...