texts: How you get the square containing a multipolyline

texts:
I use ST_Union to group a few lines. But realize some lines aren't supposed to be joined together after all. I would like to calculate the extent of the multipolyline so if the area is too big probably shouldn't have been joined together.

I need something like ST_Area but for multipolyline instead of polygons.

If I get the square I can use ST_Area


texts: Clip overlapping polygon in equal sharing size

texts:
Within single shapefile I am trying to clip overlapping polygons to their boundaries using centre line of overlapped area using Postgis or ogr2ogr.

I tried to use ogr2ogr command to remove slivers

ogr2ogr dissolved.shp input.shp -dialect sqlite \
-sql "select ST_union(ST_buffer(Geometry, 10)),common_attribute from input GROUP BY common_attribute"

But every polygon in my shapefile got only unique attribiutions so I cannot use ST_UNION.

I tried also to use ST_Difference postgreSQL database with the Postgis extension but also without success.

My main problem is that i can not find method to clip polygons on equal sharing size.


texts: While in the process of digitizing for instance a large body of water or piece of land it takes a long time to process. Is this normal? If not what can I do to speed up the process.


texts: Using the 'Mean' mergeRule with a spatial join

texts:
I have two feature classes. The first is a point file, derived from a raster, with the values for slope of the terrain. The second is a shapefile of polygons. What I want to do is spatially join the points to polygons using the ONE_TO_ONE join option, and for each polygon contain an average value of the slope of all the points that are within the polygon.

The problem is when I execute the following code, the GRID_CODE field (slope values from the points shapefile) is empty for the output feature class.

import arcpy
print "\n\n************************\n"
slopes = r'C:\Users\slope_to_pt.shp'
polygons = r'C:\Users\polygons.shp'
outfc = r'C:\Users\slopeToPolygons'
print "Creating Field Map"

#create the field map
myMap = arcpy.FieldMappings()
myMap.addTable(polygons)
myMap.addTable(slopes)

fIndex = myMap.findFieldMapIndex("GRID_CODE")
print "GRID_CODE index is:", fIndex

#create new field map with the "GRID_CODE" field
NewFieldMap = myMap.getFieldMap(fIndex)

#set the merge rule
NewFieldMap.mergeRule = 'Mean'
print "merge rule is:", str(NewFieldMap.mergeRule)
myMap.replaceFieldMap(fIndex,NewFieldMap)
print "performing spatial join..."
arcpy.SpatialJoin_analysis(polygons, slopes, slopeToPolygons, "JOIN_ONE_TO_ONE", "", myMap, "CONTAINS")

Again, when I run this code the GRID_CODE entry for each polygon is null, when it should be the average of all the values of the points that are within the polygon. Whenever I run it with a "JOIN_ONE_TO_MANY" option, it populates the GRID_CODE field; however, I want to perform a ONE_TO_ONE join.

Additional Information: the datatype for the GRID_CODE field is 'double'

Any ideas on why it wouldn't populated the GRID_CODE field for each polygon with the averages from the points?


texts: create polygon adjacency matrix using python

texts:
I am trying to create an adjacency matrix from a set of polygons. That is, I have a bunch of polygons and I want to identify which polygons have a common edge or "touch" each other. Once I find this information, I want to create an n x n matrix that indicates whether those each polygon either touches or does not touch the other polygon. Something kinda like this below.

      Poly1    Poly2   Poly3

Poly1   0        1       0
Poly2   1        0       1
Poly3   0        1       0

That is the goal anyway. I was wondering if anyone knew of a python package or such that could help me do this. I need to dynamically generate the polygons and then get their adjacency as I go, so using arcgis or qgis not the best solution.

Thanks.


texts: Spatial heterogeneity: What does spatially varying mean of a process exactly mean?

texts:
Imagine the following data set (from http://www.novaims.unl.pt/labnt/geosom/georepository/) which consists of three different regions:

Obviously, blue points are more like to be located close to other blue points. The same is true for red points. Therefore, high spatial dependency can be expected.

But what about spatial heterogeneity? Spatial heterogeneity states that the mean of a process varies with location. Going from left to right, the points (and the mean) shifts from blue to red to blue. So, this certainly is some variation. But on the other hand, comparing only the left and right side of the data, both regions, despite being distant and having different locations, have the same colour (and mean). This is certainly no variation.
So, is this data set subject to spatial heterogeneity or not?


texts: Looking for R-compatible geospatial datasets

texts:
I try to work with geospatial data of North America.
Is there any high-resolution map material available for/in R?


texts: Layer on click event not firing in Firefox but working fine in Chrome,Edge

texts:
I am using ESRI JS API 3.24 in my application.

I am using following code to open specific URL within the application when click on the layer.

var lyrUrl = "Path_To_Layer";
var myLyr = new FeatureLayer(lyrUrl, {
 id: "some_id",
 mode: 1,
 outFields: ['ID', 'Name_English']
});

//Code for the specific URL to open in a new window
myLyr.on("click", function(evt) {
 var enName = myLyr._downGr.attributes.Name_English;
 window.open('main/' + enName + '/index.html');
});

  The issue is that the click event is working fine in chrome, Edge but
  not in Firefox.

I mean not working is that when click on the layer, the link open in new page/tab but in Firefox, nothing happen at all.

Can someone help me in this?


texts: Variable width buffer based on axis directions in QGIS

texts:
I need to buffer a polygon with a different distance in the x and y directions.
I should add that my polygon has the following properties:

closed (obviously)
simple (no edge crossing, no more than 2 edges at any given vertex)
potentially concave
may or may not have hole(s) that may or may not be concave (this being said, 90% of the time there won't be any hole so I'm OK to leave this out)

I understand the buffer of a polygon is basically the Minkowski sum of that polygon and a circle of radius r (r = buffer distance), so I thought my problem may just be computing the Minkowski sum of a polygon with an ellipse.
I have read a lot on the subject of Minkowski sum here, here and there and found some algorithms easy enough to implement in Python (see below) but these are usually for convex polygons which is not my case.
I even set out to decipher CGAL's Minkowski_sum_by_reduced_convolution code but it is way too cryptic for me (I am not familiar with C).

Can you point me to a Python version or anything actually readable (pseudo-code included) of the above or similar algorithm?

For reference only (this is not what I'm trying to implement):


texts: CartoDB - PostgreSQL row_number() using PostGIS ST_Distance

texts:
I'm currently working on a query that should return a subset (i.e. "new") CartoDB table sorted by proximity to a given point. I want to display labels on the map corresponding to the closest, second closest, etc. and thought to do that by using the Postgres row_number() method in a new column:

SELECT
    *,
    ST_Distance(
        ST_GeomFromText('Point(-73.95623080000001 40.6738101)', 4326)::geography,
        the_geom::geography
    ) / 1609 AS dist,
    row_number() OVER (ORDER BY dist) as rownum
FROM locations
WHERE ST_Intersects(
   ST_GeomFromText(
      'Point(-73.95623080000001 40.6738101)', 4326
   ),
   the_geom
)
ORDER BY dist ASC

However, when I try this, CartoDB returns the following error:

Error: column "dist" does not exist

Any suggestions on a better approach or something I'm missing?


texts: Consistent Errors with QGIS using GRASS tools?

texts:
I have recently started using QGIS after using Arc for years. I downloaded QGIS Madeira 3.4.3 and have been working with road layers. After working with my project for a month in QGIS, I started trying to use some of the installed GRASS tools. Every single GRASS tool returns an error (will talk more specifically about what this is later). I searched stack exchange high and low, with no luck.

My first hunch was that there was something funky with the CRS. However, after setting the data frame's CRS to the layers, nothing changed. I then read that GRASS needs access to your .gdal folder for CRS information. I created a system variable to allow for GRASS to access the file, but there wasn't anything in the .gdal file that appeared to have CRS information. Anyway, that changed nothing.

I uninstalled and reinstalled QGIS with GRASS and GDAL using the OSgeo4 installer. Everything works great except the GRASS tools still. I have checked my topology using the check validity tool multiple times and fixed all errors.

I am trying to use v.net.distance to calculate the least cost path between points and polygons via a road layer, as it is the only tool that can use three different datasets to calculate distance and cost. When I run the tool, I get this message:

WARNING: Unable to open vector map  on level 2. Try to rebuild vector topology with v.build

When I run v.build.check, I get:

WARNING: Unable to determine input map's vector feature type(s).
WARNING: No attribute table found -> using only category numbers as attributes

To check the validity of my layers, I opened a blank map with the correct CRS and ran the Shortest Path tool from different points all over the map. It worked just fine. I then ran a whole host of other GRASS tools, and they all returned similar errors to the v.build.check errors. The problem must be something structural in the system setup as opposed to a topology error in my files.

I have been searching and searching for so long on this.


texts: How do different GIS systems determine a Polygon's interior?

texts:
I'm interested in understanding how real-world GIS systems and their data encode Polygons.

Specifically, how do they resolve the ambiguity of a Polygon's interior on a sphere?

Background: in 2D, it's trivial to pick the side of the boundary that has a finite area, since the 2D plane is infinite. However, a sphere is finite, so it's impossible to know which side is inside without making additional assumptions.

Possible approaches that I know of:

Right-hand-rule: exterior boundaries are always specified in clockwise order, and holes are specified in anticlockwise order. (There is of course the left-hand-rule, too).
Smallest area: for any given ring, always pick the side with the smallest area. I'm not sure how you would specify a large-range Polygon: perhaps an empty exterior ring followed by holes?
Equirectangular: Simply consider the equirectangular projection on an infinite 2D plane. However, this assumes features are chopped at the antemeridian, otherwise a fallback would be required to one of the two methods above.

My personal preference is the first approach, but I'm interested in understanding whether this is common in standard GIS systems.


texts: Converting points using FeatureClassToNumPyArray to interpolate cross-section with scipy "interpolate"

texts:
I'm trying to plot elevation cross section and extract elevations along a cross section using scipy's "interpolate" function.  I'm converting a point feature class of profile elevation points to 2 numpy arrays representing "distance" and "elevation", respectively.

Here is my code:

import numpy
import arcpy
from scipy import interpolate
import matplotlib.pyplot as plt

arcpy.env.workspace = "D:\\NJBPN\\volume_calc_scripts\\test_njbpn.gdb"
#interpolate profile
x = arcpy.da.FeatureClassToNumPyArray("test_prof",["Distance"])
y = arcpy.da.FeatureClassToNumPyArray("test_prof",["Z"])
f = interpolate.interp1d(x,y,kind='linear')

xnew = numpy.arange(numpy.amin(x, axis = 0),numpy.amax(y, axis = 1))
ynew = f(xnew)
interpolate.interp1d._call_(xnew)
plt.plot(x,y,'o',xnew,f(xnew),'-')
plt.legend(['data', 'linear'], loc='best')
plt.show()

It keeps giving the following error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 12, in <module>
  File "C:\Python27\ArcGIS10.2\lib\site-packages\numpy\core\fromnumeric.py", line 1987, in amin
    out=out, keepdims=keepdims)
  File "C:\Python27\ArcGIS10.2\lib\site-packages\numpy\core\_methods.py", line 14, in _amin
    out=out, keepdims=keepdims)
TypeError: cannot perform reduce with flexible type

Here is an updated code to extract the xnew value at ynew=0

xnew = numpy.arange(numpy.amin(x),numpy.amax(x)).astype(numpy.int)
ynew = f(xnew).astype(numpy.float)
shorelinedist = numpy.where(ynew==0.00,[xnew]

plt.plot(x,y,'o',xnew,ynew,'-')
plt.legend(['data', 'linear'], loc='best')
plt.show()


texts: carto db torque slider change color, font

texts:
for a theatre project im trying to make a visualization about protests in the last decade. I managed somehow to do that but unfortunately the  default time-slider of carted is really small and not easy to notice when we are using a screening.

Do you know how to change the color/size maybe the font of the time-slider?


texts: Adding Z (DEM) to my X and Y table (CSV) in QGIS

texts:
I have CSV with Name (ID), X and Y. Also have DEM. Hoping to combine my CSV that has Name with corresponding X and Y, then add the Z from DEM.
"Join Attributes from Location" and "Point Sampling" tool do not seem to do what I need.


texts: Building a QGIS 3.10 script without an output layer

texts:
So I've read about the processing scripts you can do in QGIS 3. However, there seems to be no way of just modifying a layer directly and instead I have to have an output layer, which I do not want.

If you check at https://github.com/qgis/QGIS-Enhancement-Proposals/issues/134, they say that

  it's best to have at least one out always defined to avoid fully black box algorithms if none is set. If a output isn't defined it will raise an exception.

The problem I have with that is because the module I am working on right now only adds fields and changes values in the attribute table. The output layer I would create would need to be deleted right after because there would be nothing in it.

All I want is the Input layer and make my changes within that layer, the output is a hindrance. If there is indeed no way to not have an output I would like to know if there's a way for the output to be just a log for example or something like that that way at least it could be of some use.


texts:
Almost all rasters I have in GeoServer are with float values, which means they contain at least 15 decimals. Is it possible to round these values for the getfeatureinfo?

I was looking in the data settings of my layer, but it doesn't seem like there is an option there. I am sure this must be possible.


texts: Openlayers and Geoserver fail to reproject while editing Vector with WFS-T

texts:
I can't get the Openlayers and Geoserver work with WFS-T editing a vector layer.

I have a layer of points where the Native SRS is EPSG:32639.

After trying a lot of combinations, I ended up with two not working solutions.

1)
a. On the server I set 'Declared SRS' to EPSG:32639 and make SRS hadling 'Force declared'.
b. In Openlayers i set projection: new OpenLayers.Projection("EPSG:32639") for the layer.

it submits

<gml:Point xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326"<gml:coordinates decimal="." cs="," ts=" ">241602.05063766203,4958072.672210057</gml:coordinates>

and gets

'PointOutsideEnvelopeException: 241602.05063766203 outside of (-180.0,180.0)'

2)
a. On the server I set 'Declared SRS' to EPSG:4326 and make SRS hadling 'Reproject native to declared'.
b. In Openlayers i set projection: new OpenLayers.Projection("EPSG:4326") for the layer.

it successfully commits my changes to the database, but the coordinates in the database are wrong - degrees as for 4326 instead of meters of 32639.

The Firebug shows that it submits degrees.

I have proj4js enabled on the client side.

Any other combination including playing with 'displayProjection' of the map leads to not even displaying the points of my layer on the map.

Seems like Openlayers tells projection 4326 all the time AND the geoserver is unable to reproject the coordinates when receiving them from the client.

On the server, I tried unsetting Services / WFS / GML2 / Override GML Attributes checkbox with no success.

I am in a dead end, please help.

Some code snippets:

        Proj4js.defs["EPSG:32639"] = "+proj=utm +zone=39 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";

            new OpenLayers.Layer.Vector('Wells', {
                isBaseLayer: false,
                defaultLabelField: 'NAME',
                strategies: [new OpenLayers.Strategy.BBOX(), config.saveStrategy],
                protocol: new OpenLayers.Protocol.WFS({
                    url: config.vectorWfs,
                    featureType: "Wells_PRS",
                    featureNS: "http://gisxxxx",
                    geometryName: "GEOM"
                }),
                schema: config.vectorWfs + "/DescribeFeatureType?version=1.1.0&typename=gisxxxx:Wells_PRS"
                //,projection: new OpenLayers.Projection("EPSG:4326")
                ,version: "1.1.0"
                ,projection: new OpenLayers.Projection("EPSG:32639")
            }),

            th.map = new OpenLayers.Map(th.body.dom.id, {
                controls: [
                    new OpenLayers.Control.Navigation(),
                    new OpenLayers.Control.PanZoom(),
                    new OpenLayers.Control.LayerSwitcher(),
                    new OpenLayers.Control.ArgParser(),
                    new OpenLayers.Control.Attribution()
                    ,
                    new OpenLayers.Control.MousePosition(
                        {div:document.getElementById("coordinates")})
                ],
                layers: layers
                //, displayProjection: new OpenLayers.Projection("EPSG:4326")
                //projection: new OpenLayers.Projection("EPSG:4326")
            });

            var editingBar = new OpenLayers.Control.EditingToolbar(layers[2]);
            var save = new OpenLayers.Control.Button({
                title: "Save Changes",
                trigger: function() {
                    //if(edit.feature) {
                    //    edit.selectControl.unselectAll();
                    //}
                    config.saveStrategy.save();
                // alert('saved');
                }//,
                //displayClass: "olControlSaveFeatures"
            });
            editingBar.addControls([
                save,
                new config.DeleteFeature(layers[2], {title: "Delete Feature"})
            ]);
            th.map.addControl(editingBar);

This all done reading this link:

http://www.gistutor.com/openlayers/22-advanced-openlayers-tutorials/47-openlayers-wfs-t-using-a-geoserver-hosted-postgis-layer.html


texts: Debugging Error 000725 for buffer and parameter not being valid?

texts:
I am very new to Python.

I built a model in the ModelBuilder in a toolbox in ArcCatalog 10.3.1. I keep getting the following error about my buffer not being able to be executed and that my parameters are not valid. I set a field for a distance value for my buffer, my buffer, and a spatial join output as Model Parameters in my model builder before I exported the script to Python to run in Python IDLE:

  Traceback (most recent call last):   File
  "F:\Projects\Working\HurricanesSCRIPT.py", line 49, in
      arcpy.Buffer_analysis(Hurricanes_in_Counties, Hurricanes_Buffer, v15_Miles, "FULL", "ROUND", "NONE", "", "PLANAR")   File "C:\Program
  Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\analysis.py", line 692, in
  Buffer
      raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000725: Output Feature Class: Dataset
  F:\Projects\Working\ModelProject\ModelData.gdb\Hurricanes_Buffer
  already exists. Failed to execute (Buffer).

I have already put in arcpy.env.overwriteOutput = True just below import arcpy.

Here is the script for my arguments:

# Script arguments
v15_Miles = arcpy.GetParameterAsText(0)
if v15_Miles == '#' or not v15_Miles:
    v15_Miles = "15 Miles" #15 Miles

Hurricanes_Counties_SpaJoi = arcpy.GetParameterAsText(1)
if Hurricanes_Counties_SpaJoi == '#' or not Hurricanes_Counties_SpaJoi:
    Hurricanes_Counties_SpaJoi = "F:\\Projects\\Working\\ModelProject\\ModelData.gdb\\Hurricanes_Counties_SpaJoi" # provide a default value if unspecified

Hurricanes_Buffer = arcpy.GetParameterAsText(2)
if Hurricanes_Buffer == '#' or not Hurricanes_Buffer:
    Hurricanes_Buffer = "F:\\Projects\\Working\\ModelProject\\ModelData.gdb\\Hurricanes_Buffer" # provide a default value if unspecified

And here is the script for my buffer:

# Process: Buffer
arcpy.Buffer_analysis(Hurricanes_in_Counties, Hurricanes_Buffer, v15_Miles, "FULL", "ROUND", "NONE", "", "PLANAR")


texts: ArcSDESQLExecute seems does not work with startTransaction and commitTransaction

texts:
I'm doing some calcultations by ArcSDESQLExecute. The idea is to make some part localy and at the end copy  table/featureclass to SDE and make some extra operations by sql. I want to put all sql to transaction,but  it seems like it does not work as a transaction.
Here is example of code, which works wrong.

sde_connector = arcpy.ArcSDESQLExecute("path to connector")
sql = [ u"ALTER TABLE {tab_mc} ADD bez_UR number(38,8)".format(tab_mc = mc_results),
        u"UPDATE TP_MC_STAT set test_col = 1000",
        # u"DROP TABLE TP_MC_STAT",     ##WHEN TABLE IS NOT DELETED, NEXT SQL RAISE AN EXCEPTION
        u"CREATE TABLE TP_MC_STAT (test_col int)",
        u"UPDATE TP_MC_STAT set test_col = 4"]
sde_connector.startTransaction()
for row in sql:
    try:
        sde_connector.execute(row)
    except Exception as err:
        raise Exception(row," ...>>>... ",err.message)
sde_connector.commitTransaction()

When I run full code without commenting "DROP TABLE" everything is ok and column test_col has value 4. But When I  comment  DROP statement I except no sql is gonna be commited  and test_col has still value 4. But after raising exception there is 1000 in test_col - previous sql were commited. How is this possible?


texts: What are the differences between Anchored, Framed and FramedCloud?

texts: From the API document for OpenLayers, I can tell that a FramedCloud has more functionalities than Framed, which in turn has more functionalities than Anchored. I am hoping for someone to describe these differences in a way which a lay person like me can understand.


texts: Is there a way to set different levels of transparency (opacity) for a polygon's border and its interior?

texts:
I have some polygons which I'd like to display with their interior displayed at a higher transparency than their borders. Right now, I am doing it by creating two separate layers with the same data source, one for the border (20% transparent) and one for the interior (70% transparent). Actually, I usually end up creating 3 different layers in these cases: one each for the polygon's border and fill, and a third to use in the legend, because neither of the other two will work perfectly in the legend.

Here's an example, in this project I have floodplains set to 70% transparent in the fill (because many other items I'm interested in are in the floodplain and I need to see them), and 0% transparent on the border so that its boundaries are clearly visible. Is there a way to set the transparency of the polygon's border separately from that of the interior, so I don't have to duplicate layers?


texts: Getting particular feature using expression in PyQGIS?

texts:
How to get a particular feature using an expression in PyQGIS?

For example, there are 100 houses (with 100 entries in Attributes Table). I want to get only one house with field name ("name"=="John's House")  and make get selected that feature in the map view.


texts: Tool For Generating Formula for Projection

texts:
I am using Proj4Js for northing,easting to lat,long transformation.

My Datum is "AGD84" and projection is "Transverse Mercator".

While doing so I found that if I don't specify "+proj" parameter in the definition file of source coordinate system then it gives me error saying "Projection Code Not Found".

I found that "Proj4Js" is having some ready made library of projection code but nothing out of them is useful to me.

Can anyone suggest any tool in which if we pass datum & projection name then automatically it generates projection specific formula which can be usable with Proj4Js ?


texts: Visibility Analysis Plugin (QGIS 3.18.3) Returning multiple errors

texts:
I am trying to run a viewshed analysis in QGIS 3.18.3 using the Visibility Analysis plugin. I have followed the steps entirely in both of these videos:
https://youtu.be/pt1gFXON7W4
https://youtu.be/fGKgvg0QoXA
but it continued to throw an error. Both the "viewpoint" and the DEM are in UTM coordinates. Any ideas on how to fix this? I ran the same process in ArcGIS Pro and it worked fine.
Here is the error log that I get when I try and run the analysis:
QGIS code revision: 735cc85be9
Qt version: 5.11.2
GDAL version: 3.1.4
GEOS version: 3.8.1-CAPI-1.13.3
PROJ version: Rel. 6.3.2, May 1st, 2020
Processing algorithm…
Algorithm 'Viewshed' starting…
Input parameters:
{ 'ANALYSIS_TYPE' : 0, 'DEM' : 'C:/Users/Arda/Documents/Summer 2021/Viewshed_Paper/New folder/clip_aster.tif', 'OBSERVER_POINTS' : 'C:/Users/Arda/Documents/Summer 2021/Viewshed_Paper/New folder/point_loc_shp.shp', 'OPERATOR' : 0, 'OUTPUT' : 'TEMPORARY_OUTPUT', 'REFRACTION' : 0.13, 'USE_CURVATURE' : False }

Traceback (most recent call last):
File "C:/Users/Arda/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ViewshedAnalysis\algorithms\viewshed_raster.py", line 251, in processAlgorithm
refraction = refraction )
File "C:/Users/Arda/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ViewshedAnalysis\algorithms\modules\Raster.py", line 151, in set_master_window
self.window = np.zeros((full_size, full_size))
MemoryError: Unable to allocate 9.21 PiB for an array with shape (35998455, 35998455) and data type float64

Execution failed after 0.07 seconds


texts: PostGIS: ST_Azimuth between Points in same Table

texts:
I've seen the questions and answers about ST_Azimuth, but didn't find a specific example for this case:

I need to calculate the azimuth from point id 1 to id 2, id 2 to id 3, id 3 to id 4, ...

Table 'points' with id and geom

How can I define the query to calculate the azimuth between id and id+1 (next point)?

Or is there a QGIS field calculator solution?


texts: ogr2ogr: Transformation from NAD27 to WGS84 not working

texts:
I'm trying to use ogr2ogr to transform from NAD27 to WGS84 and check the precision afterwards.
For that, I used a GML file with coordinates -124.0, 48.0

It looks like this:

{
  "type": "FeatureCollection",
    "crs": {
    "type": "name",
    "properties": {
      "name": "urn:ogc:def:crs:EPSG::4267"
      }
    },
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates": [
          -124.0,
          48.0
        ]
      },
      "properties": {}
    }
  ]
}

Then I tried to convert it using ogr2ogr by running this command:

ogr2ogr -f "GeoJSON" -s_srs EPSG:4267 -t_srs EPSG:4326 output.geojson input.geojson

It returned the same coordinates -124.0, 48.0:

{
   "type":"FeatureCollection",
   "crs":{
      "type":"name",
      "properties":{
         "name":"urn:ogc:def:crs:OGC:1.3:CRS84"
      }
   },
   "features":[
      {
         "type":"Feature",
         "properties":{

         },
         "geometry":{
            "type":"Point",
            "coordinates":[
               -124.0,
               48.0
            ]
         }
      }
   ]
}

However, when I tried to convert between NAD27 and WGS84 using epsg.io web page and called this URL:
http://epsg.io/transform#s_srs=4267&t_srs=4326

It returned these coordinates: -124.001314 and 47.9998062

There is a big difference between both results [-124.0, 48.0] and [-124.001314, 47.9998062] - around 140 meters!

Does it mean, the ogr2ogr is not precise or I didn't use it properly?


texts: How to create a circle with specific radius?

texts:
Possible Duplicate:
  Conflicting measurement units with buffer ftools

I have a shapefile as input. How could I create a circle with a radius of exactly 500 meters, using QGIS?


texts: Identify all continuous networks in a dataset

texts:
I have a set of pedestrian streets for a city. I would like to assign all of the streets that can be reached to one another along a continuous network to a unique group. I am wondering if there’s an off the shelf method (QGIS, ArcGIS, PgRouting) that can do this.

I have an image of roughly what I would like to accomplish. Basically, I have a set of roads on an “island” that are broken up barriers that make walking impossible (Interstates, Rivers, etc.).  I want the continuous part of each network to be assigned to its own group. For reference, I am trying to create "T-Communities" described in Grannis, Rick. "T‐Communities: Pedestrian Street Networks and Residential Segregation in Chicago, Los Angeles, and New York." City & Community 4.3 (2005): 295-321.


texts: Calculating number of pixels with different values

texts:
I have a tiff where each pixel has a value in megagrams for biomass.
I'm trying to get a count of the pixels for all the pixels with diffrent megagrams values so that I can get an estimation of biomass.

How do I go about this in QGIS 3.2?


texts: Field calculator not calculating all feature areas: Geometry is NULL

texts:
I had various layers and through geoprocessing tools "Clip" and "Difference" I have created one layer in QGIS 2.16. This layer has 98 attributes - polygons. I want to calculate the area of each of the polygons and only retain those that are larger than 10ha:

Open Attribute table -> toggle edit -> Add field -> Name = "Area" (real).

Field Calculator ->

I do obtain the calculated area for the first 38 polygons, but not for the remaining ones. I tried to locate the ones without calculated area and obtain this error:

What does this error mean?

Is it even possible to calculate the area for the remaining ones? If so, How?


texts: Merging two Landsat TIFF file of same band in QGIS 3.4

texts:
I've been wanting to merge two TIFF files of same band (this case it will be B1 in Landsat). I have followed tutorials which I came across online, but none of them seems to work. Generally, raster > miscellaneous > merge is the procedure. But it only ends up showing one TIFF file. Should I create a virtual raster and then try? Or maybe it's not possible to merge any single band (e.g. B1,B2 etc.)?

Also please tell me if I'm using the right output file type.
Here are few screenshot of the problem:

(In the last screenshot, i clipped the raster with a shapefile to check whether it merged properly.)

This is  kind of similar to my problem which was posted here earlier but no possible solution came out of it: Merging landsat 8 tiff images in QGIS?


texts: Exclude Raster Cells with Value of Zero in ArcGIS 10.2 Raster Calculator Expression

texts:
I need to perform a simple calculation in ArcGIS 10.1 Raster Calculator. I am using a raster layer that contains many cells that have a value of zero. I need to exclude these cells in the expression because they are skewing the results.

Here is the simple raster calculator expression that I am using:
"rho"* 1005*(((("raster")*.02)-272.15)-"temp")/30)

I need to add to this expression so that only values greater than zero in the raster layer "raster" are consider in the expression.


texts: Getting from PDF contour map of Z in area to values of Z specific points in area?

texts:
I have a PDF of a contour map of value Z for an area. For the same area, I have a set of points that I would like to match with Z values interpolated from this contour map. For an individual point, this can be "eyeballed" pretty easily, but I would like to do this for a few hundred points so I am looking for a more efficient way of doing this using for instance Google Earth, arcGIS, r, etc.

I am aware of the ability to overlay images in Google earth and draw lines on a map to export as kml files.

This seems like a reasonable place to start but it seems overly tedious to import these lines into r separately and write some algorithm to interpolate between the lines.

I am wondering what workflows and existing tools people are aware of for this type of problem since I have struggled to find much.

I am very new to working with GIS data.


texts: RGeo::Feature.cast does not convert between projections

texts:
I am attempting to cast a MultiLineString feature in WGS84 to UTM33N using RGeo.

Unfortunately, RGeo::Feature.cast, instead of casting the feature to the desired projection, returns the original feature unchanged, as the checks at the end show. Why does the conversion result like this? What else is needed to convert between these two projections in RGeo?

require 'rgeo'
require 'rgeo-geojson'

utm_33n_factory = RGeo::Geographic.spherical_factory(srid: 32633, proj4: '+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
wgs_84_factory = RGeo::Geographic.spherical_factory(srid: 4326, proj4: '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')

geo_json = '{"type":"Feature","properties":{"City":"PRG","Name":"Z","Color":"green","Description":"Pohorelec"},"geometry":{"type":"MultiLineString","coordinates":[[[14.389935,50.087667,0],[14.389219,50.08756,0],[14.389093,50.087525,0],[14.388617,50.087394,0]]]}}'
json_geom = RGeo::GeoJSON.decode(geo_json, geo_factory: wgs_84_factory, json_parser: :json).geometry
projected = RGeo::Feature.cast json_geom, factory: utm_33n_factory, project: true

puts json_geom.equal? projected # => true
puts projected.srid # => 4326


texts: Difference between source SRS and target SRS in Raster Warp function

texts: I am trying to change a raster layer projections from WGS84 to a custom Lambert layer using 'Raster' --> 'Projections' --> 'Warp' function. I have been asked for source SRS and get SRS. what is the difference?


texts: QGIS Atlas map canvas overlap calculation

texts:
My QGIS Atlas uses a fixed map canvas size and a fixed scale. The atlas generation shapefile has polygons in a steady interval. Therefore my maps do have an overlap area. My aim is to calculate said overlap. The intial idea was to get the coordinates of the map canvas extent and form a polygon from that. Now I want to erase this polygon from the $currentfeature polygon and this should be my overlap area. My current non-working and only pseudo code is like this:

from qgis.core import*
from qgis.gui import*
from qgis.utils import iface

@qgsfunction(args=0, group='Custom', usesgeometry=True)
def mapExtent(value1, feature, parent):
  bb = iface.mapCanvas().extent()
  gPolygon = QgsGeometry.fromPolygon(bb)
  clipstuff = clip(gPolygon, $currentfeature)
  overlap = $area(clipstuff)
  return overlap

The result should be shown in a label inside the atlas. I am not so good at coding.


texts: eCognition export samples as a shapefile problem

texts: I need sample shapefile for the other project but when I try to export it, it gives me an error like domain is empty(no images object available). How can I get the samples as shapefile?


texts:
I am trying to load a WFS from GeoServer. My GeoServer service is at localhost:8080 and my web server at localhost:80. Since there is the issue of cross domain and firebug is showing:

XMLHttpRequest cannot load http://localhost:8080/geoserver/wfs.
Origin http://localhost is not allowed by Access-Control-Allow-Origin

I have created a proxy.cgi following http://www.gistutor.com/geoserver/21-intermediate-geoserver-tutorials/38-configuring-geoserver-proxy-for-public-and-remote-data-access.html. But when I go to http://localhost/cgi-bin/proxy.cgi its giving me the OpenLayers site, and when there is no proxy and just leading to http://localhost/xampp, when there is a proxy (here proxy is my institute proxy to connect to internet)..

My problem is I am able to load the data from GeoServer successfully (as if everything is fine) but firebug is showing:

"POST *http://localhost/cgi-bin/proxy.cgi?url=http%3A%2F%2Flocalhost%3A8080%2Fgeoserver%2Fwfs undefined (undefined)".

I am including OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url="; in my JavaScript and also:

protocol: new OpenLayers.Protocol.WFS({
            version: "1.0.0",
            srsName: "EPSG:4326", // this is the default
            url:  "http://localhost:8080/geoserver/wfs",
            featureType: "states",
            featureNS: "http://www.openplans.org/topp"
        }),

I want to know why it's working fine but  still giving me the error, What is the error about and how can I to fix it?


texts: Get Census Tract given a Census Block OR Get Census Tract given Latitude/Longitude

texts:
I'm using Python and am looking to get the Census Tract of a location given either a latitude and longitude or a census block. I can get the census block given latitude/longitude from here: https://geo.fcc.gov/api/census/#!/block/get_block_find but I don't know how to go from census block to census tract.
I need census tract data because I'm working with the Social Vulnerability Index data which is census-tract level.


texts:
I am looking for some free dotNET libraries which can be used to convert an OpenStreet .osm file to .shp shapefile.

I did Google, but mostly found either .exe files, or non.NET files.

Can you point me to the .NET libraries?


texts: R: shapefile extent larger than raster extent

texts:
I am using R for my spatial analysis. I am fairly new to this so I am still trying to find my way around things like changing projections etc...

I have a raster with vegetation index values (0.5 x 0.5 degrees), and a shapefile with polygons. I wanted to use intersect() or something along those lines to get the vegetation index values for my polygons, but I am running into issues with the extent or projection - not exactly sure what is causing the issue.

At first it became clear the .tif didn't have a projection, so I set the projection the same as my shapefile. Testing with compareCRS() gives TRUE. When I try to do intersect() though, I get the error

Warning messages:
1: In intersect(x, y) : non identical CRS
2: In intersect(x, y) : polygons do not intersect

When inspecting the raster and shapefiles, I see indeed that the extents are different - but the extent of my shapefile is actually larger than my raster (and my raster covers the whole world, I plotted it and that's correct). The shapefile is a region in the Amazon, so it should fall inside the raster covering the whole world.
How does that happen? Can I change the extent of my raster without messing up the resolution or the grids?
I've done a lot of searching online, but most of the resources talk about cropping when the extent of one file falls within the other file...
Here are the details of my files (r is the raster, comms the shapefile):

> show(r)
class       : RasterLayer
dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 7200, 0, 3600  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /EVI2_testdata/EVI2_1983152.tif
names       : EVI2_1983152
values      : -32768, 32767  (min, max)

> show(comms)
class       : SpatialPolygonsDataFrame
features    : 99
extent      : 320877.6, 703893, 8549725, 8916662  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables   : 9
names       : Join_Count, TARGET_FID,                    P5_1, TabulateAr,
Shape_Leng,   Shape_Area, Cluster_ID, ClusterID1, Cluster1v2
min values  :          0,          0,               3 Arroyos,          1,    276.5641, 1.010235e+08,          1,          1,          0
max values  :         19,       4282, Virgen de la Candelaria,        103, 217329.5623, 9.944712e+07,          4,          4,          0

Below the code used to get to this point:

library(gdalUtils)
library(raster)
library(rgdal)

# Get a list of sds names out of HDF4 file, make a tif
sds <- get_subdatasets('VIP30P2.A1983152.004.2014270210035.hdf')
filename <- 'EVI2_1983152.tif'
gdal_translate(sds[1], dst_dataset = filename)
# Load the Geotiff created into R
r <- raster(filename)

# read in shapefiles
comms<-readOGR(dsn=path.expand("~/MAP_PhD_research/DATA/Shapes/CommunityPolygons_new/MAP_Tenure_WithComminutyPnt_20130709.shp"), layer="MAP_Tenure_WithComminutyPnt_20130709")
# check projections
identicalCRS(comms,r)
# which is false
proj4string(comms)
# result: +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
# assign same projection to raster file
crs(r)<-"+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(r)
compareCRS(r,comms) # result is true
# still intersect does not work
intersect(comms,r)


texts: Spatial metadata table in SQL Server 2008

texts: I have a db and want to know if it contains any spatial tables. Is there a table in SQL Server that holds spatial metadata i.e. a list of all the spatial tables in the current db? Oracle has the user_sdo_geom_metadata, Postgis the geometry_columns. What about SQL Server? The only spatial metadata I found was the sys.spatial_references.


texts: Join multiple shapes in a workspace based on a common field

texts:
I’m completely stucked!!! I’m working with meteorological data (wind) but I have to arrange the data for a program to be able to process it. I have text files, two per day (one with the wind direction and the other with the speed).

I’m working with ArcGis 10, Model Builder and my little knowledge on python.

I’ve already manage to calculate in a point shapefile the date and the information, if it’s direction or speed, because it was information available only in the filename.
But now, I don’t know how to join the two files together. The fields schema I have it’s like:

POINT_ID / DATE / WIND (dir or spd) / ID_DATE (unique field) // DIR / SPD

I have like 1000 days to process and I have no idea how to continue!! The goal would be to make a join based on the ID_DATE field, so I have a unique file for each day instead of two, with the pertinent values of DIR and SPD (already calculated for each individual file).

I need to implement an attribute join with the two tables (dir and spd). The problem is that I have like 1000 different days, therefore 2000 tables to join.

I was thinking on something like reading the filename, where the date is within the name, and when finding two identical dates performing the join… So I need to read through the folder (or two separate folders with dir and spd) and indicate somehow the common attribute field or filename….

Thanks in advance


texts: Does CQL requires bringing the data on the client?

texts:
I want to use ECQL in order to filter some data on my browser. I have a postgis table which will have around 300k of rows.
Bringing this layer at once might make things a bit slow.
I was wondering how ECQL works. When you filter something, does it bring the whole layer on the client and then filters it? Or it sends requests to the server and returns only what is asked?


texts: arcpy: changing lyr data source (raster layers)

texts:
I am trying to create a lyr file  for each raster file in a folder.
To get exact same display settings I thought I would create a lyr file, then, using arcpy copy-paste it for each file in folder and change the data sources.

The script generates new lyr files but does not change data sources - and gives no error. So, what am I doing wrong?

import arcpy
lyrFile = arcpy.mapping.Layer(r"C:\temp\rst\sample.lyr")
arcpy.env.workspace=r"C:\temp\rst"
for i in arcpy.ListDatasets():
lyrFile.replaceDataSource(lyrFile.workspacePath,"RASTER_WORKSPACE",i)
lyrFile.saveACopy(lyrFile.workspacePath+i)


texts: How to I mask out the black from a GeoTIFF

texts:
New to Google Earth Engine and I'm wondering why my GeoTIFF's come convert the alpha/transparency of the image to black.  I have these weird black outlines from my TIFF files.
See Screen shot.

This is orthomosaic from drone data and this is what the TIFF looks like in photoshop.


texts: Save raster to file when using stack() in R

texts:
I'm working on a function for raster package. In the end of the function I have to use stack() to create a RasterStack.  When I do this the data source is in memory. I want to have the data source on a file.

For example if I use crop() I can do:

asd <- crop(raster, extent, filename="something.tif")
inMemory(asd)
TRUE

I can't do the same with stack(), but I need to have the results on files, because it writes all files on RasterTmpFiles and I end having memory problems.


texts: How do you display related records in QGIS?

texts:
Coming from an esri background I have always used relationship classes in a gdb to display records related to a feature. Eg. Relate 1 survey point to many documents, photos etc.

Does any one know how to do something similar in QGIS?

Ando

P.S I'm using PostGreSQL 8.4.1 and PostGIS 1.4 and my clients are QGIS.


texts: Is this pgRouting perfomance acceptable for my tiny server?

texts:
I'm testing pgrouting functions in a tiny server with 2GB of RAM and a dual core processor.
I wrote a program that create client process and they connect to the server in a specific hour to perform a query, in fact the query is a pgrouting function call. My trouble is that after 50 process, the server crash and report this error in the postgres log:

  terminating connection because of crash of another server process
  2012-09-05 15:59:09 VET DETAIL:  The postmaster has commanded this
  server process to roll back the current transaction and exit, because
  another server process exited abnormally and possibly corrupted shared
  memory. 2012-09-05 15:59:09 VET HINT:  In a moment you should be able
  to reconnect to the database and repeat your command.

and this error occurs for every process.

The table that contains the graph information has about 300,000 rows. The graph has 250 nodes aprox.

My question: is it ok that this behavior occurs with a server like this and using a table with this amount of rows?


texts: I am trying to create new shapefiles from an existing shapefile using QGIS 2.18. On getting to the select by expression box which I am doing by states, a state name appears twice. I have checked the attribute table of the shapefile I am creating the selection from and checked the spelling as well but I seem can not find what the problem is. How do I rectify this?


texts: Are there any open source distribution raster merging functions within JTS or GeoTools

texts: I'm looking for something similar to the combine function from the ISelectionSet interface from ArcObject's SDK 10 .NET framework on a more scaled down level. Ergo, I don't plan on using that interface to combine objects that manage rows or features. If I need to I'll do a cast as some type of table like object, but nothing drawing from ArcObjects.


texts: Automatize adjust polygon bounds based on other polygon

texts:
I have two polygons, as below:

Does exist an automatic way to adjust the vertices of Pink Polygon near of Blue Polygon bounds to Blue Polygon bounds?


texts: Is there anywhere I can download the new USGS Topo TIFS?

texts:
Is there anywhere I can download USGS Topo TIFS?  The new ones?  The USGS site just gives me GeoPDF and their web service restricts the scale you can see them at.
Edit : I am looking for the the new USGS Topo maps (c 2012).


texts: Python Built-in Map and Zip function in QGIS 2.0 Processing Script--syntax error

texts:
I am attempting to run my Python script using the Processing Toolbox in QGIS 2.0. I've used the built-in Python function map() and zip() several times and am now receiving a syntax error even though these functions worked fine in python console. Thanks in advance.

*I was asked to post a little more of my code, so i've edited and added a different chunk as I don't run into these errors when loading the vector layers by pathname…This makes me wonder if there's something wrong with the way these layers input into the Processing script. Mahalo:

from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from processing.core.VectorWriter import VectorWriter
import numpy
import processing

##input1=vector
##input2=vector
##Attribut=field input1
##output=output file

indivLayer = processing.getobject(input1)
grosLayer = processing.getobject(input2)

indivprovider = indivLayer.dataProvider()
grosprovider = grosLayer.dataProvider()

#Break grosLayer dictionary into list of lists (keys represent names of zones)
for key in CG_IDdict.keys():
    listofzones = CG_IDdict.values()

#Get sum of list of zones and put in dictionary with CG_ID keys
sumList = map(sum, listofzones)
sumdict = dict(zip(CG_IDlist, sumList)
#Count in zone
countList = map(len, listofzones)
#Get max and min value of zones
maxValue = map(max, listofzones)
minValue = map(min, listofzones)
#Get average value
meanValue = map(numpy.mean, listofzones)

Error message:

Traceback (most recent call last):
File "C:\Documents and Settings\user/.qgis2/python/plugins\processing\core\GeoAlgorithm.py", line 150, in execute
self.processAlgorithm(progress)
File "C:\Documents and Settings\user/.qgis2/python/plugins\processing\script\ScriptAlgorithm.py", line 213, in processAlgorithm
exec(script) in ns
File "<string>", line 74
countList = map(len, listofzones)
^
SyntaxError: invalid syntax


texts: How do I pass a LatLng Bounds from Leaflet to Mapfish Print?

texts:
I successfully fired up a MapFish print server on Elastic Beanstalk, and it prints out the sample spec files fine, but I'm trying to pass bounds from my Leaflet frontend now, and encountering a lot of difficulty.

The OSM spec provided uses scale and center, but the docs indicate that bbox is another option allowed, so I thought I could use that. But whenever I try to modify said OSM config to use my working bounds from LeafLet, all I end up with is a blank PDF output.

I have been experimenting with Proj4 to convert these latlng points in meter units with the Google Projection used in spec below, but that isnt working out too well for me either.

So my base question is, how would I modify the sample,working, OSM spec, seen below, to allow for lat/lngs bounds?

I figure I would have to modify the SRS to 'EPSG:3857' (since that is what Leaflet is using), and units to 'degrees'...but how do I prep my latlng numbers?

{
layout: 'A4 portrait',
title: 'A simple example',
srs: 'EPSG:900913',
dpi: 190,
units: 'm',
layers: [
    {
        baseURL: 'http://tile.openstreetmap.org/',
        type : 'Osm',
        maxExtent: [-20037508.3392,-20037508.3392,20037508.3392,20037508.3392],
        tileSize: [256,256],
        resolutions: [156543.0339,78271.51695,39135.758475,19567.8792375,9783.93961875,4891.969809375,2445.9849046875,1222.99245234375,611.496226171875,305.7481130859375,152.87405654296876,76.43702827148438,38.21851413574219,19.109257067871095,9.554628533935547,4.777314266967774,2.388657133483887,1.1943285667419434,0.5971642833709717],
        extension : 'png'
    }
],
pages: [
    {
        center: [731000,5864000],
        scale: 100000,
        mapTitle: "OSM Map",
        comment: "This is the first page selected by the user.",
        rotation: 0,
        data: {
            data: [
                {id:1, name: 'blah', icon: 'icon_pan', nameBackgroundColor: 'red', nameBorderColor: 'blue'},
                {id:2, name: 'blip', icon: 'icon_zoomin', nameBackgroundColor: 'yellow', nameBorderColor: 'green'}
            ],
            columns: ['id', 'name', 'icon']
        }
    }
]

}


texts: How to reach the geometry (polygon, polyline) of a dxf in python-ogr module

texts: this might be a basic question, but I am rather new in ogr. I wish to collect all lines from the dxf file and create their boundary. I loaded the layers (ExecuteSQL( "SELECT DISTINCT Layer FROM entities" )); and would like to put every line in a new layer, but it seems that I cannot reach the geometry of the layers. I tried to use layers.GetFeature().GetGeometry() in a for, but it is not working.


texts: Adding mdb geodatabase layer (feature class) to QGIS project using python?

texts:
I want to add a feature class of a personal geodatabase mdb to my qgs project by code (a plugin I am working on).

I am using a command similar to that used for adding shapefiles, but this is not working.
I guess there must be a way, since it is possible to add the layer using the Add Vector Layer button.

This is the code i am using on the Python console:

mylayer = QgsVectorLayer("D:\Test\test.mdb\points", "mypoints", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(mylayer)

The first line runs without error but the second one does not produce any result, contrary to the result is I do a similar command using a shapefile.

I am using QGIS 2.2 Valmiera 32 bits on Windows 7 64 bits. I did a previous research but I did not find any results!!! The feature class is not inside a feature dataset.

I guess I would have to use a QgsDataSourceURI(), but i can not figure out the syntax. I have been trying. Also, I would also want to add that migrating the mdb is not an option, I have requested this before but my superiors have refused. And this geodatabase (points) is updated daily and huge. My goal after loading is clipping it based on a zone of interest.


texts: How can I make GeoServer's WMS layer .jar file in Android development Android studio

texts:
I'm new in Android development in geoserver's WMS map.

Can I get .jar file of GeoServer's WMS layer?

I couln't open it though I tried https://github.com/shalperin/android-wms.git


texts: How can I fix QGIS PostGIS Manager error 'shp2pgsql' is not recognized?

texts:
I have QGIS 1.8.0 installed via OSGeo Installer and work on win 7.
All I need is to upload the shapefile to the database via this plugin. Unfortunately I am receiving an error:

Something's wrong:
["'shp2pgsql' is not recognized as an internal or external command,\n", 'operable program or batch file.\n']

How to install and configure shp2pgsql to make this plugin work fine?


texts: Finding high-resolution land cover data for Taiwan?

texts:
So far, the best data I've found is from Globeland 30 (30 m), but the data for Taiwan is not very discriminant in partitioning the landscape.

Does anyone know of better resolution data?


texts: Lost pixels during a raster merge with gdalbuildvrt

texts:
I'm trying to merge some rasters with gdalbuildvrt in a python script:

rasters = glob.glob('/{0}/ndvi_masked/*.tif'.format(path))

gdal_cmd = 'gdalbuildvrt -q -vrtnodata "0" -overwrite {0}.vrt'.format(ndvi_name)

for pix in rasters:
    gdal_cmd += ' {0}'.format(pix)

subprocess.call(gdal_cmd, shell=True)

gdal_cmd = 'gdal_translate -of GTiff %s.vrt  %s_clip.tif' % (ndvi_name, ndvi_name)

subprocess.call(gdal_cmd, shell=True)

But the result is really weird, cause I losted some pixels as you can see on the 2nde image below.

All displayed rasters

Merge raster


texts:
I need to copy a feature class to a vector database but the message I got is that the feature class already exists.  However, I don't see it in the vector database - it is not there.  When investigating in SQL Management Studio, the table does exist.  So how do I fix up this can anyone please help?

  ERROR 000210: Cannot create output Database
  Connections\Staging_forupdate@vector.sde\vector.GIS.Gov_Housing\Vector.GIS.Stock


texts: Editing shapefile from QGIS in ArcMap causes errors according to QGIS topology check?

texts:
I noticed Polygon shapefiles that were edited with ArcMap yield many errors in the QGIS Topology checker.

The picture below shows a polygon Shapefile. I made one error on purpose so the topology checker shows there is just one error (checking for overlap within the layer)

Next I edit this layer in ArcMap, creating a hole between two of the polygons

Then I use the QGIS topology checker again, yielding 3 additional errors.

Now my questions are:

Can I change some settings in ArcMap so that it will not create this
difference? I think it might be saving the vertices of the edited
polygon with a lower precision.
if it is not possible to fix it in ArcMap, then what is the quickest method to clean these tiny errors quickly afterwards with either QGIS or PostGIS?


texts: Tools for sorting/filtering based on map state

texts:
I'm about to embark on a project that list restaurants on a map. I'd be grateful to find out what the best tools for the job are.

The restaurants will be in a list beneath a map. The list will be refreshed based on the map zooming. The map will also have selectable/hoverable bounding boxes, so the user can select a location (London) and then drill down on the location to get neighbourhoods/boroughs in that area (Westminster, Clapham etc.)

I'm looking at the following tools for the job - Gisgraphy (for data), OpenStreetMap (for maps and data) and leafletjs (for the boxes). I've looked into GeoServer but this seems to be something I can't host locally.

Is there anything I'm missing, or are there any better tools to achieve this? Also are there any examples I could look at where this has been done already?


texts:
I am working in a work space in python on ArcGIS 10.1 on a remote desktop. Here is my code

import arcpy
from arcpy import env
env.workspace = "E:\Example"

PST = arcpy.ListFeatureClasses ("", "Point")
LST = arcpy.ListFeatureClasses ("", "Polyline")

for point in PST:
    for polyline in LST:
        arcpy.SplitLineAtPoint_management (polyline, point, point + "500" + ".shp", 1)

I want the output feature class to have the same file name as the point feature class, but with 500 appended to the end. Currently the way the code is now, it will name the output "Example.shp500.shp" which is giving me an error of failing to execute. How can I name my output "Example500.shp"?


texts: With leaflet-routing-machine, how to bring a route on top?

texts:
On a simple Leaflet map, with leafleat-routing-machine, if two routes pass by the same place, only the top one is visible. I can live with this if I can choose which route is on top, (for example by clicking on it). Go, given a route, how can I put it on top?

Here is a simple code sample I would like to improve:

map = L.map("map");

L.tileLayer("https://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png", {
    attribution: "&copy; OpenStreetMap contributors"
}).addTo(map);

var route1 = L.Routing.control({
    waypoints: [
        L.latLng(-22.985308, -43.204845),
        L.latLng(-22.983449, -43.202773),
    ],
}).addTo(map);

var route2 = L.Routing.control({
    waypoints: [
        L.latLng(-22.983449, -43.204845),
        L.latLng(-22.984708, -43.202430),
    ],
}).addTo(map);

// Now I would love this to exist:
function bringOnTop() {
    ...
}

bringOnTop(route1);
bringOnTop(route2);


texts: Spatial Interpolation Methods like PRISM

texts:
I am looking for some kind of software (free software) that allows me to make an interpolation like PRISM software.

The interpolation used in PRISM is a

  local regression and accounts for spatially varying elevation relationships; also accounts for effectiveness of terrain as barriers, terrain-induced climate transitions, cold air drainage and inversions, and coastal effects

What I am aiming to interpolate are climate time-series datasets.


texts: Where can I download a free dataset containing major ports across the world?

texts:
I am working on an analysis and would like to incorporate major maritime ports from across the world. Some initial searching turned up a dataset produced by General Dynamics, however it will be prohibitively expensive.

I would like to find a free dataset to use, preferably in shapefile or some other Arc friendly format. The only information desired beyond the location (point or polygon) would be the name, however even that is not required. Any help is much appreciated.


texts: Plot a circle with a given radius around points on map using Python

texts:
I have created a map using GeoPandas with a few markers that represent some coordinates on the map. I want to do draw a circle around each marker with the marker placed in the center of the circle.
This is what the map looks like:
This is pretty much what I am trying to do (a circle around each yellow marker):
This is my code so far:
#I have a pandas DataFrame called address_df that has a long and a lat column

address_df['geometry'] = address_df.apply(
lambda x: Point((x.long, x.lat)), axis = 1)

#changing to geopandas
geo_df = gpd.GeoDataFrame(address_df, crs=crs, geometry= address_df.geometry)

#the in_cou is a geoDataFrame that I read_file from a .shp file
in_cou.geometry = in_cou.geometry.to_crs(epsg = 3857)
geo_df.geometry = geo_df.geometry.to_crs(epsg = 3857)

#Joining datasets
cou_hos = gpd.sjoin(in_cou, geo_df, how="inner", op='intersects')

#The figure is create using matplotlib
fig, ax = plt.subplots(figsize=(15,15))
minx, miny, maxx, maxy = in_cou.total_bounds

ax.set_xlim(auto=True)
ax.set_ylim(auto=True)
ax.set_aspect('auto')
#Thicks left empty for cosmetic reasons
plt.xticks([])
plt.yticks([])
plt.title("Providers",fontsize=40,y=1.04, x= .55)

plt.rcParams['font.size'] = 20
cmap1 = ListedColormap(['#2f7d50'], name='allred')
cmap2 = ListedColormap(['#a0a8a3'], name='allblue')

#Plot the first layer of the map.
in_cou.plot(ax=ax, edgecolor='black',cmap = cmap1, legend = False, alpha =0.9)
cou_hos.plot(ax=ax, edgecolor='black',cmap = cmap2, legend = False,alpha = 0.9)

#Plot the second layer of the map. The yellow plus markers.
geo_df.plot(ax=ax, marker="P",edgecolor='k', markersize=150,alpha= 0.8, c= '#ffff00', label= "Providers" ) #geo_df['Type'].apply(lambda x: colors[x]))
plt.legend(prop={"size": 19},bbox_to_anchor=(0.56, -.09), loc='lower center', ncol=3)

plt.show()

I have created this function based on what I found in a previous question (https://gis.stackexchange.com/a/359748/166972). It creates a circular polygon based on a given latitude, longitude, and a radius (The radius is also a column with int values) :
def max_poly(lat, long, dist):
    lon, lat = long, lat # lon lat for San Francisco
    radius = dist  # in m

    local_azimuthal_projection = "+proj=aeqd +R=6371000 +units=m +lat_0={} +lon_0={}".format(lat, lon)

    wgs84_to_aeqd = partial(
        pyproj.transform,
        pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
        pyproj.Proj(local_azimuthal_projection))

    aeqd_to_wgs84 = partial(
        pyproj.transform,
        pyproj.Proj(local_azimuthal_projection),
        pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"))

    center = Point(float(lon), float(lat))
    point_transformed = transform(wgs84_to_aeqd, center)
    buffer = point_transformed.buffer(radius)

    circle_poly = transform(aeqd_to_wgs84, buffer)
    return circle_poly

Then I create a new geoDataFrame with polygons created by max_poly
crs = {'init':  'epsg:3857'}
new_df = gpd.GeoDataFrame(crs=crs, geometry=geo_df.apply(
    lambda x: max_poly(x.lat, x.long, x.max_dis_trav), axis = 1))

But plotting new_df does not work.


texts: I want to convert my shp layers into one DGN file using ogr2ogr (QGIS). When I use this tool, I can only select one layer to convert. My problem is I don't know how to select all three layers at the same time and converting them into the same destination dgn file.


texts: Printing/displaying coordinate system of several raster files using ArcPy?

texts:
How can I extract the Coordinate System of many raster files (GeoTIFF) so the information is displayed, and I can ascertain that all of them are correctly settled in the same projection, using ArcPy?

The ideal would be not to load the files as it is too time consuming.


texts: How to acces features attributes when creating renderer?

texts:
I'm trying to customize the rendering of a point layer in QGIS using python. Thus I have written a renderer class as explained in the documentation (http://www.qgis.org/pyqgis-cookbook/vector.html#creating-custom-renderers).

From what I have understood, I should define which symbol to render for each feature using the method 'symbolForFeature'. My problem is that inside this function when I try to get the attribute map using 'attributeMap' I always get an empty dictionary.

class MyRenderer(QgsFeatureRendererV2):
  def __init__(self, syms=None):
    QgsFeatureRendererV2.__init__(self, "MyRenderer")
    self.syms = syms if syms else [ QgsSymbolV2.defaultSymbol(QGis.Point), QgsSymbolV2.defaultSymbol(QGis.Point) ]

  def symbolForFeature(self, feature):
    attrs = feature.attributeMap()
    # attrs[0] should be the id
    if attrs[0] > 10 :
      return self.syms[0]
    else :
      return self.syms[1]

  def startRender(self, context, vlayer):
    for s in self.syms:
      s.startRender(context)

  def stopRender(self, context):
    for s in self.syms:
      s.stopRender(context)

  def usedAttributes(self):
    return []

  def clone(self):
    return MyRenderer(self.syms)

When I use this class, I get the error: 'KeyError: 0', simply because attrs is empty.

However, the other methods of the feature object seem to provide the expected results, for instance when I replace 'attrs[0]' by 'feature.id()' everything works fine.

What am I doing wrong?


texts: QGIS - Rotating all symbols to face a single point

texts:
I am creating a map with airports, using an airplane symbology (svg).
Instead of having all small plane symbols turned north, I would like to rotate them to be all turned to Paris.
I have 28 airport points in one layer, and the point of Paris in another layer.

I assume I need to draw lines between each point to Paris, then calculate some rotation angle that I can then use in the symbology.

How is this best achieved?


texts: Custom Symbols are too small in KMZ File from ArcMap

texts:
We are working with the local Police department on an iOS and Android app which will have an embedded Google map showing crime data from the past month.

This map is generated via a .KMZ file exported from ArcGIS.

We created a set of custom .EMF icons for this map and the Crime Analyst at the Police Department imported them fine, but they are extremely small in the KMZ layer that gets added to the Google map.

It looks like the icons get saved to the .KMZ file as 16x16 .PNGs.

When the Crime Analyst increases the size of the icons in ArcGIS, it makes them pixelated and blurry.

I've also tried unarchiving the .KMZ file and inserting larger versions of the custom symbols into the folder and then zipping it up and converting it to a .KMZ file, but then the Google Map wouldn't display any of the custom symbols.

Is there a simple way to retain the original size of the imported custom icons or to increase their size without pixelating them?

We have the custom icons as .SVG, .PNG, .EMF, and as a custom icon font. Is it possible to use either the .SVG or the icon font as a custom marker in ArcGIS so that scaling will not be an issue?

Let me know if you need more information or to see some of the files. I don't know much about the Arc suite of software and am just the app developer, but I can get information from the Crime Analyst.

EDIT:  Here are the answers to your questions.

What tool are you using to export the kmz? - In arctoolbox, conversion tools, "layer to kml" tool
What is the projection of your map document? - NAD_1983_StatePlane_South_Carolina_FIPS_3900_Feet  Lambert_Conformal_Conic
What is the unit of the projection? - Feet
Imported them into what? - Select a generic symbol. Edit symbol in Symbol Property Editor. Symbol Type is Picture Marker Symbol. Then select the .EMF file.


texts: Visualization of pictures for points in QGIS

texts:
I'm working with QGIS Desktop 3.2.3.

A Colleague of mine worked on a Google Map and inserted in his Google Map Pictures. He exported the KML and gave it to me so that I can import it to QGIS.

The pictures are part of the attribute table as a URL. Is it possible to download the pictures into QGIS?


texts: Calling specific columns of attribute table using ArcPy?

texts:
I'm trying to figure out how to spline specific columns of a shapefile attribute table, but my question is intermediate to this.

How do I select a range of columns, e.g. the 20th column to the right-most column, as fields to apply a for loop operation to?

My first thought was to try numpy, but I'm open to ideas. For now I'd just like to either print the names of the desired columns or specify the columns to be iterated in the for loop line.

import arcpy, numpy as np
from arcpy.sa import *

fieldList = arcpy.ListFields(r"C:\Users\afullhar\Desktop\Lab6\Data For Lab on habitat modeling\Layers\start_data\10_km_Points.shp")

for field in fieldList:
    print format(field.name[19,:])


texts: ArcPy Calculate Field management syntax error

texts:
I have list of GDBs each containing Mauza and Khasra Feature classes I want to update fields using string
code  is as below
import uuid
import arcpy as ap
import os
DistrictName='Faisalabad'#'SUBMISSION 11-05-2020/'
district='D:/ASADABBAS/Asad/DOA/Outsource/Recevied/'+DistrictName

def getGdbPaths(districtPath):
    gdbNotFound=[]
    directories=[os.path.join(districtPath, o) for o in os.listdir(districtPath) if os.path.isdir(os.path.join(districtPath,o))]
    gdbs=[]
    for directory in directories:
        for villagepath in os.listdir(directory):
            if 'WITH EXCEL' in villagepath or 'WITHOUT EXCEL' in villagepath:
                for subvillagepath in os.listdir(directory+'/'+villagepath):
                    for gdb in os.listdir(directory+'/'+villagepath+'/'+subvillagepath):
                        if not gdb.endswith('.mxd'):
                            if (gdb.endswith('.gdb')):
                                gdbs.append(directory+'/'+villagepath+'/'+subvillagepath+'/'+gdb)
                            elif (gdb.lower()=='mouza'):
                                for subgdb in os.listdir(directory+'/'+villagepath+'/'+subvillagepath+'/'+gdb):
                                    if subgdb.endswith('.gdb'):
                                        gdbs.append(directory+'/'+villagepath+'/'+subvillagepath+'/'+gdb+'/'+subgdb)
                            else:
                                gdbNotFound.append(directory+'/'+villagepath+'/'+subvillagepath+'/'+gdb)
                                print('No GDB found in',directory+'/'+villagepath+'/'+subvillagepath+'/'+gdb)
            else:
                try:
                    for gdb in os.listdir(directory+'/'+villagepath):
                        if not gdb.endswith('.mxd'):
                            if (gdb.endswith('.gdb')):
                                gdbs.append(directory+'/'+villagepath+'/'+gdb)
                            elif (gdb.lower()=='mouza'):
                                for subgdb in os.listdir(directory+'/'+villagepath+'/'+gdb):
                                    if subgdb.endswith('.gdb'):
                                        gdbs.append(directory+'/'+villagepath+'/'+gdb+'/'+subgdb)
                            else:
                                gdbNotFound.append(directory+'/'+villagepath+'/'+gdb)
                                print('No GDB found in',directory+'/'+villagepath+'/'+gdb)
                except NotADirectoryError:
                    print('Invalid directory path',directory+'/'+villagepath)
    return gdbs,gdbNotFound
gdbs,noGDB=getGdbPaths(district)
for gdb in gdbs:
        ap.env.workspace = gdb
        mauzaid=(str(uuid.uuid4())).upper()
        expression='{'+mauzaid+'}'
        ap.CalculateField_management("Mauza",'ID',expression,"PYTHON3")
        ap.CalculateField_management("Khasra",'ID',expression,"PYTHON3")

when I run this code it gives me following error
SyntaxError: invalid syntax (<expression>, line 1)
Failed to execute (CalculateField).

Type of the field I am trying to update is GUID. How can I solve this issue?
PS: I am able to access feature classes but CalculateField_management does not seem to work.


texts: OpenLayers 3 / Cesium polygon projection issue

texts:
I am using OpenLayers 3 and Cesium (via ol3-cesium).

I have created a few rectangles (using a MultiPolygon) stacked on top of each other. Whereas in OpenLayers the rendering (left) is as expected, the output in Cesium (right) is not.

I understand that a rectangle in one projection will not be a rectangle in the other. However, I would expect the software to appropriately project the rectangle.

In this case, what I was expecting is for the Cesium output not to have any gaps between rectangles. On the other hand, the vertices seem to be in the right coordinates, so these polygons are just following the shortest lines between vertices.

What defines which projection is a polygon based upon, and therefore which projection should apply to it? For these vertices, are the right polygons the left or the right ones?

And moreover, how do I get Cesium / ol3-cesium to render the polygons covering the same area as the equivalent polygons in the 2D map?

JSFiddle: https://jsfiddle.net/adlerhn/o1j1fv9x/4/


texts: Changing page layout from portrait to landscape using ArcPy with ArcMap?

texts:
I need my script to access the Page Orientation property of the current MXD, and switch it from portrait to landscape (or back).

I've looked through the ArcGIS 10 arcpy and arcpy.mapping documentation, and haven't had any luck.

Does anyone know how to do this?


texts: Test whether coordinate reference system has units of meters

texts:
Is there a simple way to test whether a spatial object's coordinate reference system expresses coordinates in units of meters? The snippet below shows my initial attempt at a function that provides such a test, but I don't have a great sense for how broadly reliable it is.
Is a better (more reliable and/or more widely tested) approach already available? Failing that, are there commonly used coordinate reference systems for which my approach is likely to fail?
library(sf)

## Function to test whether an sf object's projection expresses
## coordinates in units of meters by checking that the projection:
## (a) is not longlat; and (b) comes with no multiplier factor for
## converting units to meters.
is_proj_in_meters <- function(obj) {
    !is.na(st_crs(obj)) &
        !st_is_longlat(obj) &
         is.null(st_crs(obj)$to_meter)
}

## Create four sample polygons, each with a different CRS

g <- st_as_sfc("POLYGON ((-61.66957 10.69214, -61.565 10.75728, -61.37453 10.77654, -61.40721 10.60681, -61.66957 10.69214))")

## (1) No CRS supplied (units are unspecified)
a <- st_as_sf(data.frame(id = 1, geometry = g))
st_crs(a)$to_meter ## NA
st_is_longlat(a)   ## NA

## (2) WGS 84 projection (units are degrees)
b <- st_as_sf(data.frame(id = 1, geometry = g), crs = 4326)
st_crs(b)$to_meter ## NULL
st_is_longlat(b)   ## [1] TRUE

## (3) 1903 Trinidad grid (unit is Clarke's foot)
c <- st_transform(b, 2314)
st_crs(c)$to_meter ## [1] 0.30479726540000002
st_is_longlat(c)   ## [1] FALSE

## (4) UTM zones applied to Trinidad and Tobago (units are meters)
d <- st_transform(b, 32621)
st_crs(d)$to_meter ## NULL
st_is_longlat(d)   ## [1] FALSE

## Finally, test `is_proj_in_meters()`
sapply(list(a,b,c,d), is_proj_in_meters)
## [1] FALSE FALSE FALSE  TRUE


texts: Multiple entries of same layer in GeoPackage

texts:
I have a GeoPackage of many layers which I have added and deleted over time. I think because of that, I get several entries of the same layer. For example when I use the statistics tool, I get a list of all my layers in my GeoPackage, but most layers appear up to 5 times. It makes it difficult to find the right layer.

I guess each version (new and old) gets stored within the GeoPackage whenever I add a layer with the same name even though the layer exist only one time.

How can I clean up my GeoPackage?


texts: How to calculate the bounding box of a map, given the pixel width/height of the container map & a bounding box which is included in the container map

texts:
How can I calculate the bounding box of a map, given the pixel width/height of the container map and a bounding box which is included in the container map?

I know how to do this with Google Maps JavaScript API... but I want to do this only with a calculation.


texts:
I'm trying to use ArcPy and Python more since the ability to just point and click isn't that marketable. I've been using ArcPy for Select by Attribute. It works fine when I try to select attribute with numbers (float, double etc) like bigger than (>) or smaller than (<). But when I try to include text/string in my expression, I always get the "ERROR 000358: Invalid expression"

For instance:

arcpy.SelectLayerByAttribute_management('Subregional Groups', 'NEW_SELECTION', '"MUNICODE" > 13')

The code above gives me the result I want.

But

arcpy.SelectLayerByAttribute_management('Subregional Groups', 'NEW_SELECTION', '"NAME" = "Allentown City"')

gives me the error mentioned before.

I have also tried

arcpy.SelectLayerByAttribute_management('Subregional Groups', 'NEW_SELECTION', '[NAME] = "Allentown City"')

and it doesn't work either.

I'm using ArcGIS Pro 2.2.


texts: How can I make a polygon feature selectable on the edges, but not in the enclosed area

texts:
My map has markers in a layer and can select any marker by clicking on it.
I also have a layer that contains polygon features, and I can also select these by clicking on them - either the edge or anywhere in the enclosed space.

When a polygon is over a marker, then I can't click on the marker - without hiding or disabling the polygon.

Is there a way to make the polygon feature selectable by clicking on its edges, but transparent to selection in the enclosed area?

Responding to comment from https://gis.stackexchange.com/users/56638/jbalk
This is a different question to How to define layer order in Openlayers?  however, the solution to that question does point to a good solution here. If the markers are on top of the polygons, then I expect that the markers would receive the click events on the markers but the polygons would receive the click events between the markers. This is not the same as 'make a polygon selectable on the edges', but it would solve my issue. Does that make it a duplicate?


texts: Will ArcGIS 10 work properly on Linux using wine or crossweaver?


texts: Text + Chart inside pop up window in Javascript

texts:
A similar question had been posted in this thread :
Chart inside pop up window in Javascript

I'm trying to do the same, but i would like to add others things in the popup (like text, feature properties, pictures, ...).
Unfortunately, when i add elements in the popup, i get [object HTMLDivElement] instead of graph.

Based on this code (http://jsfiddle.net/6UJQ4/), i wrote this :

m.bindPopup('text
' + div);

I think my problem is related to the fact i don't understand the writting of this statement :

var div = $('

')[0];

(i am totally newbie in javascript !)


texts: Count raster value in a bounding box (spatial extent)

texts: Given a bounding box, how can I count the number of specific values (say, I am interested in the number of value == 1 ) in a raster in 1) arcpy, or 2) raster package in R?


texts: Can I automatically save datetime for every entry in my spatialite tabel?

texts: I'm using QGIS and SpatiaLite. Is it possible to store date and time in a field "datetime" for every entry automatically - maybe with some python code? Already existing datetime should not be updated.


texts: What is the order of parameters in +towgs84?

texts:
What is the order of parameters in +towgs84 attribute within a proj string. I assume it is dx,dy,dz,rx,ry,rz,sf
dx,dy,dz are expressed in meters
rx,ry,rz are these in degrees/radians/arc seconds?
sf value units as well would be great


texts:
What is the most pythonic way to create a Word File for a .png image, having only it's bounding box?

For example, given the following bounding box:

<LatLonBox>
  <north>30.46370000000000</north>
  <south>25.15370000000000</south>
  <east>89.06550000000000</east>
  <west>83.06550000000000</west>
</LatLonBox>

the result should be a text file reporting, in order:

pixel size in the x-direction
rotation about y-axis
rotation about x-axis
pixel size in the y-direction
x-coordinate of the center of the upper left pixel
y-coordinate of the center of the upper left pixel


texts:
I routinely have to mosaic multiple GeoTIFFs together, extract a rectangular area, and save that area as a JPG with a world file. Typically I've been using commercial software for that (Geographic Imager) but would like to switch that task over to QGIS.

I can eventually get the task completed using multiple tools, but the process takes much longer and involves creation of intermediate GeoTIFFs along the way. In Geographic Imager I simply load the images into Photoshop, crop, and SaveAs and I'm done.

So what's the easiest, quickest workflow to:

Load multiple GeoTIFFs
Mosaic them together
Select an area to extract (preferably by drawing a box)
Save selected area out to JPG with accompanying world file


texts: How can the gridlines in QGIS be set to be curved lines (neat circles) around polar region (azimuthal projection)?

texts:
I am trying to present a map of our polar regions using QGIS software.  I am using an azimuthal (polar) projection.  I have added a graticule (gridlines) to the map as a background, which are at 10 degrees lat/long intervals.

However, instead of having neat circles around the pole, the mapping output shows straight lines between the 10 degree longitude intervals.  Is is possible to get a curved line between the longitude intervals?  I realize that a workaround is to create a shapefile with the 10 degree gridlines and upload it to my map, but I would prefer a more elegant solution.


texts: Making POINT from longitude/latitude coordinates in existing table using PostGIS?

texts:
I would like to make a POINT using the longitude and latitude columns of an existing table.  Here's a link to the documentation of ST_MakePoint.

Here's my table:

CREATE TABLE sysinst.bio (
    oid INTEGER PRIMARY KEY,
    longitude FLOAT,
    latitude FLOAT,
        geom GEOMETRY(POINT, 26913)
);

Adding Geometry Column:

SELECT AddGeometryColumn ('sysinst', 'bio', 'geom', 26913, 'POINT', 2);

Here's my query:

UPDATE sysinst.bio SET geom = ST_SetSRID(ST_MakePoint(longitude, latitude), 26913);

No returned error though the POINTS are wrongly plotted when viewed in QGIS.

For (-97.5959, 21.1922) , I recieve (-109.4896, 0.0002)

Any suggestions?

I'm open to doing this more efficiently if there's an easier way.


texts: Calculating distance between previous and current addresses for large file in QGIS?

texts:
I'm a fairly new QGIS user.

I have a large dataset (5000+ records) of customer addresses, previous and current, and I want to calculate the distance moved for each customer. Ideally this would be via roads, not as the crow flies, but this is national-level data and I don't yet have any shapefiles for roads covering the whole country (US). I don't have the file yet, but I believe it will have lat/long.

is there any way to do this without locating/downloading a national street file?
How would I go about solving the problem either way?


texts: Openlayers GeoJSON write not maintaining projections

texts:
How do I get GeoJSON.write to maintain the projections used in the construction of the GeoJSON?

points.geojson

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "name": "carey",
        "lure": "",
        "date": "",
        "pic": ""
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -93.9001,
          49.213196
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "dan",
        "lure": "",
        "date": "",
        "pic": ""
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -93.89787,
          49.209946
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "tony",
        "lure": "",
        "date": "",
        "pic": ""
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -93.89032,
          49.21174
        ]
      }
    }
  ]
}

The layer is created with the following

 var vlayer = new OpenLayers.Layer.Vector( "Points" , {
              projection: new OpenLayers.Projection("EPSG:4326"),
              strategies: [ new OpenLayers.Strategy.Fixed(), saveStrategy],
              protocol: new OpenLayers.Protocol.HTTP({url: "./points.geojson", format: new OpenLayers.Format.GeoJSON({ignoreExtraDims: true }) })
         });

Writing

    var format = new OpenLayers.Format.GeoJSON({ projection: 'EPSG:4326', internalProjection: new OpenLayers.Projection("EPSG:4326"),  externalProjection: new OpenLayers.Projection("EPSG:4326") })

    console.log(  format.write( vlayer.features ) );

Produced GeoJSON

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "name": "carey",
        "lure": "",
        "date": "",
        "pic": ""
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -1.0452912E7,
          6311114
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "dan",
        "lure": "",
        "date": "",
        "pic": ""
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -1.0452663E7,
          6310560
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "tony",
        "lure": "",
        "date": "",
        "pic": ""
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -1.0451822E7,
          6310865.5
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "dan",
        "lure": "",
        "date": "",
        "pic": ""
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -1.0461912E7,
          6321375.5
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "dan",
        "lure": "",
        "date": "",
        "pic": ""
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -1.0465505E7,
          6317706.5
        ]
      }
    }
  ]
}


texts: On save add timestamp in attribute table?

texts: I would like to add / update a timestamp value in a column of the attribute table whenever this attribute table is saved. Should I do it with an own python script or is there any in-built functionality which can perform this for me?


texts: Calculating Cut & Fill volume using ArcGIS Desktop?

texts:
I work with ArcMap 10.1.

There is a proposed area for building a dam. I'm just trying to calculate the approximate volume of dam in different elevations.

I created a raster of area and also other several rasters with different fixed elevations (50, 75,100 ...) and tried to calculate cut and fill volume through this order. But after calculation, when I open the attribute table the numbers don't make any sense. I attached the picture of table.

I have these questions:

Is this even a right way that I've chosen to calculate volume of a reservoir? If it is not how can I calculate it?
How can I find out what is the unit of numbers that were given in attributes table (cubic meter?, cubic kilometer, 1000cubic meter?, ..)?


texts:
My desktop app contains a Map using ArcGIS runtime SDK. The app should work in a fully disconnected environment. The app will be installed on a Glider (Aircraft).
I have done all the app on Qt and the map provided via ArcGIS-Server using Map, MapGraphicsView Classes.

I found an example LocalServerFeatureLayer but this is also not working when I am disconnected only (map didn't show)  when connected (the map start working ).
So my question: Is there is a way to get the map of a country downloaded offline and installed into a local server on the same pc that contains the my app. If it is available what is the process to do that.

I am using Qt-C++ , SDK-100.8


texts: Changing shapefile to stereographic projection in QGIS?

texts: I have a shp file of the world's coastlines in QGIS and I was wondering how I would convert it to a Stereographic projection as shown in this map:


texts: Is there an easy common gis abstraction available in .Net

texts:
Is there an easy common gis abstraction available in the .Net programming world that cuts the dependency to a concrete gis implementation, e.g. to the esri silverlight/wpf sdk.

... for example OpenLayers is such a library for several map providers in JavaScript. By using an IoC container (like Spring.Net) the dependency could be hidden. I noticed the GeoAPI buth this seems too complex for an easy abstraction level, I think.


texts: Two raster classes according to cut-off elevation using ArcPy?

texts:
I've defined a python function to take input raster and elevation and output path. The function should return raster with two values 0 for elevation equal to and greater than given elevation and 1 for rest of the pixels. This is what I've done so far:

import arcpy
arcpy.env.overwriteoutput=True
from arcpy.sa import*

def getRaster(input, elevation, outputpath):
    inraster=input
    cutoffelevation=elevation
    arcpy.checkExtension('Spatial')
    outputraster=Raster(inputraster)>cutoffelevation
    outputraster.save(outputpath)


texts: Creating decision rules on pixel level in time series layer in R?

texts:
I am using Landsat 5, 7, 8 images and RStudio. My goal is to compute the disturbance in a given area using the Vegetation Change Tracker (VCT) algorithm. I would like to ask you if you could provide me a general guidance on how to construct it using R Studio. I computed the IFZ of each individual image, I created the LTSS, I have the NDVI and the NBR of each year. My problem is on how to create the rules required form the VCT (see the image below).

The acronyms CLIV, MNCLIV and MNCHIV means:

the pixel was classified as a forest at least in 2 sequentially images,
the minimum years of sequentially images needed for a pixel to be classified as a forest (2 years) and
minimum number of consecutive high IFZ values for a pixel to be classified as a non-forest pixel.

The paper named An automated approach for reconstructing recent forest disturbance history using dense Landsat time series stack.


texts:
I'm trying to install gdal-1.11.0 on CentOS 6 machine.
Following the instructions, this is my config

./configure --prefix /usr/local  --with-fgdb=/usr/local/FileGDB_API --with-pg=/usr/pgsql-9.3/bin/pg_config

As the configure process runs, I see this in the output:

checking for PostgreSQL... yes

However, at the end of the configure output, I see this:

PostgreSQL support:        no

PostgreSQL 9.3 is up and running with postgis on it, the path to pg_config is correct. I know the configure is reading it, because if I put a bad path in --with-pg, it complains...
I had prior installs of earlier versions of GDAL, and GDAL 1.11 without postgreSQL before I realized that I need to build it with PostgreSQL support.
What am I missing?


texts: Save Selected Features in QGIS 3.4 Model

texts:
I'm trying to build some models I have in QGIS 2.18. Each model uses the "Save Selected Features" geoalgorithm so that the model only runs on selected features in the source vector layer.

In 3.4 there is no "Save Selected Features". There is an "Extract Selected Features" but I am unable to get the "polygons to lines" algorithm to point to it as it does in 2.18.

This is the 2.18 model


texts: How can I prevent the scale bar in print composer from bunching up in QGIS?

texts: I'm creating a series of maps in print composer from the same basemap (from the quickmapservices plugin) in QGIS 2.12.2. Sometimes this works, though it only allows the scale to be in metres, feet or nautical miles and not km, which is not ideal, however other times the numbers bunch up over each other such that they are not readable (see image below). How can I prevent this?


texts: Features placed 7.5 thousand kms off correct location using QGIS?

texts:
I have been counting muskrat dens on lakes with aerial photos by creating a feature in a separate shapefile.
I checked my progress after doing about 80 small lakes and suddenly, half of the lakes were off of their correct location (7.5 km south) and half of them are in the right spot.

I can't seem to remember, but I think that I wanted to confirm that the projection was correct and I confirmed the right projection of the layer.

Is it possible that this confirmation changed the first half of my layers and when I continued with the rest, they were placed correctly?

Is there anyway to reproject some of the features in a shapefile?


texts: Anyone who has the ESRI Census 2000 Tiger/Line data?

texts:
The ESRI Census 2000 TIGER/Line Data were free.

http://www.esri.com/data/download/census2000-tigerline
http://arcdata.esri.com/data/tiger2000/tiger_download.cfm

But they vanished. The download link is broken. Is there anyone out there who has the data?


texts: Raster or Brick from nedf with lat lon = matrix

texts:
I have a netcdf file 3 dimensions, lat, lon, time :

lat and long are matrix (nrow=263, ncol=512)
time is vector of date
the variable is precip observed at a given (lat, lon) and date

What I would like is a brick (nrow=263, ncol=512) with the values = precip , one layer per date . That's what I try :

test=brick("netcdf_filename.nc", varname="precip")

And that's what I get :

Error in .rasterObjectFromCDF(x, type = objecttype, band = band, ...)
:    cells are not equally spaced; you should extract values as points

I also try to specify the grid properties like this :

brick("netcdf_filename.nc",varname="precip",nrows=263, ncols=512,
ymn=min(lat),ymx=max(lat),xmn=min(long),xmx=max(long),
crs="+proj=longlat +datum=WGS84 +ellps=WGS84")

but I get the same error. I'm not an expert in handling geographic data , is there a solution to get a brick with this kind of netcdf file ?


texts: How to check if REST plugin is installed for GeoServer

texts:
I am having trouble uploading vector files onto GeoServer 2.4 using Import Vector in the QGIS Processing toolbox. After reading around about the 'REST' plus in I went to: http://<companyname)/geoserver/rest. However I got a 404 error as in it didn't exist even for existing files. I tried to run the same thing on my localhost and it worked.

Also when I run the tool I find an error code 404 in the log.

GeoServer is running on a company server along with Tomcat Apache 6.0

Could someone please tell me how to check the REST plug in, and if it is not installed how to get the REST plug in for GeoServer 2.4.0?


texts: Trying to query based on selected feature, but query is using previously selected feature

texts:
I'm not sure how to replicate this problem with just a code snippet, because it involves so many different elements of my HTML and JS. So here is a jsbin with the whole code.

I use ESRI's popup widget to populate the sidebar with data from the selected feature, and it works. Except for one bit where I'm running a function pctile after selecting the feature (line 1376 in my JavaScript).

pctile should give me a value P that is the percentile rank of the currently selected tract.

If you watch the console while clicking census tracts in the map, you can see that line 1408 writes the value P of the previously selected tract ("P from buildHeader"). After that, the function pctile fires, updating the value of P to the current tract ("P:").

How can I ensure that pctile fires and updates P before I try to use it in my function buildHeader (line 1403)?

UPDATE I think the solution involves using callbacks. I'm trying to make that work now.

Update I have tried using callbacks, but I'm not sure that solution works in this context. Because the function I want to pass as the callback argument is my buildHeader function that returns HTML to insert into a div. From what I understand you can't get a return from a callback function.


texts: Seeking API to check OpenStreetMap for village boundaries?

texts:
I am new to GIS programming and I have started a little project.

I have a lot of gpx tracks, and I want to check which villages i have visited.

Is there an API to check OpenStreetMap for village boundaries?

I would like to get those village boundaries as polygons, because it is easy to check if a track goes through polygons using Python.


texts: QGIS OpenStreetMap plugin: download failed. Server replied: bad request

texts:
I am trying to download Open Street Map data using the OpenStreetMap plugin in QGIS, and getting the following err:

I'm kind of new on geoprocessing, so I have no clue what is going on here.


texts: I'm interested in locating free GIS shapefile data representing the US / Mexican border. I have found a point layer depicting the border crossing points but I'm interested in a polyline feature representing the actual border.


texts: Is it possible to use more than two control point links when georeferencing CAD files in ArcMap?

texts: I'm trying to georeference CAD to a building footprint in ArcMap 10 (although I've had this same problem in earlier versions). The problem is, it won't let me use more than two control point links, so I have to update the display before adding two more. Every time I update the display, the CAD is stretched and placed properly where the control points were, but the other ends of the CAD are off. Is this simply something I'll have to live with? Or is there a way I can use more control points, or possibly convert the CAD to a less restricted file format?


texts: How to digitize 3D geospatial data using FOSS applications

texts:
I need to digitize soil horizon data from photos to be later used in a 3-dimensional GIS workflow.

For each photo I have DGPS coordinates of the top left and right corner and the depth of the profile. Thus, I can place the photos in three dimensional space. I want to digitize the soil horizons on each of these photos, combine it with horizon-specific data and then be able to extract the data per coordinate and depth.

Is there a FOSS application that can drape images in 3d and be used to digitize information in 3d vector format? I usually use QGIS for digitizing and R for all my data processing, but I am fine with other FOSS solutions on Linux as long as I can get my data into these applications later on.