texts: How To Determine Percent of Matched Data?

texts:
In my attribute table, I have 2 columns with qualitative (categorical) data. Both of them are from different sources so I am trying to find the percent of matched data from these 2 columns. I created a new column labeled 'Matched', and wrote the following codeblock in the Field Calculator to determine the values, but I am getting an error on line 2-3. Do I have to convert them to string using str()?

def matchedRecords(x):
    x = (column.name)
    x2 = (column2.name)

    if x == x2:
        resultIs = "Matched"
    else:
        resultIs = "Unmatched"
    return resultIs

And then matchedRecords(!Matched!) in the second area of the Field Calculator. I want to determine the total matched values versus unmatched values to find the percentage. Is this the best way to do it? Thanks.


texts: Arcgis ModelBuilder - OD cost matrix on network with iteration on attributes

texts:
I am totally new to ArcGis' ModelBuilder, and I need some advice as how to build an OD cost matrix for large sets of origins and destinations.

I want to compute an OD matrix where I can have as many as 20,000 points as origins and as destinations. Using the normal procedure with Network Analyst results in an "out of memory error". That's why I turned to modelbuilder.

I managed to (re-)build the core model that computes the distance between OD pairs with the network dataset (niger_ND) and the points (p100). The output is written in the database reseau.mdb. This works fine.

Now, I want to add some features to my model, but I do not know exactly where I should put them:

Each point in the layer containing the points has a unique identifier. As origins, I want to select only 1000 of them at a time. When the analysis is done, I want to relaunch the process on the next 1000 points until every point has been used. How should I design the model to use my points one batch at a time?

With each batch of origins, I save the result in a database. How can I avoid overwriting the previous result by giving a different name to the output at each iteration?

I guess it is not too complicated when you know ModelBuilder, but I do not fully understand its logic yet.


texts: Intersecting Floating Points within a value range to parcel boundary

texts:
I have an NDVI with values of -0.99 to 0.88. The raster is floating points with no attribute table. I also have farms as parcel polygons.
Using ArcMap or ArcGIS Pro, how can I intersect only the floating pixel values of 0.08 to 0.20 that are on my specified parcels? I want to identify the farms that are changing crops, and last year that was the value range.


texts:
What would be a way to compute moving averages of a time series in GEE? Consider we imported an image collection X and filtered by start and finish dates.

var Y =ee.ImageCollection(X.filterDate(start,finish));

I want to be able to compute moving averages of the time series


texts:
Anyone know of any open source .NET or Javascript libraries for creating GeoPDFs?

I know of QuickPDFLibrary, but i'm preferably looking for something open source.

I want to programmatically generate some GeoPDFs.


texts: Create Points at 10 Meter intervals along line with QChainge

texts:
I have a route network and want to add every 10 meters a point with the 'QChainage' plugin.

I first used 'join multiple lines' to make the route a single section.
Then I transferred my CRS to my corresponding UTM-zone, to work with meters.

Changed 'Chainage every ' to 10 Meters and pressed Ok.

As result, however, I get only one Point.

The 'Devide Feature into x Parts' under 'Advanced ' seems to work, but it is not what I am looking for, since I want specific Distance.

Does anybody have an Idea what I can do, to make it work? I searched similar question regarding my task, but the answer always suggests to use the plugin, which I can't get to work.


texts: Calculating "travelled distance" of a GPS track without including "stops"

texts:
I have a set of GPS tracks from a cycling event and I want to correlate their distance along the route.

Problem is: when a rider stops, generally by entering a covered place, GPS keeps logging (low precision) points around its location, and by summing those noisy segments I get a false increase in travelled distance.

If I look at the plotted tracks, it's easy to "see" that actually the rider didn't move, but how could I create some cumulative distance algorithm that "doesn't count" some points because the rider is "obviously not moving" there?

One interesting thing I noticed: if I consider only one coordinate, for example, only latitude for a road heading straight north, the fuzzy points at a stop keep wandering back and forth around a definite latitude, the same applying for longitude (although in this case that would be the same longitude of every point of that road segment, because it is "vertical" on the map).

But if I use euclidean distance, then this random movement is converted always into "forward" movement and keeps increasing the distance.

If the stop were in the middle of a straight road, I could get the average direction and "normalize" and average the positions along the road direction, but oftentimes the stop point is at a sharp turn (a cross-road) or it is the turning back point of a route.

Is there any known method for this, or does anyone have a suggestion?


texts: SCP V.4.9.6 - after creating ROI, Warning: The following signature will be excluded if using Maximum Likelihood

texts:
I am creating ROI in Clouds using "create a ROI" (polygon) and this warning often shows after i save the ROI. Since the warning says that it will affect the Maximum Likelihood, that is why i delete the ROI that i created. But sometimes there is no warning.

Any advice why this happened? is it the size of the polygon i created? is it necessary to delete it?


texts: Using the Split tool on a single polyline with strange results

texts:
I am having issues with the Split tool.  When I try and split a line at a given point my result is not  a new feature where the break was meant to occur.

This is my line (one feature)

This is where I place the cursor for the split to happen.

The result is two lines but they are not as expected....I want one line North of where I put the split tool and another South of where I put the split tool....instead I have a small segment of line two in the middle of line one.

Has anyone seen this before or know what I might need to change in order for this split to work as I want it to?


texts: How to check if a feature class is being used in MXD?

texts:
Trying to remove any unused feature classes contained throughout multiple geodatabases.  Want to check and see if each feature class is being used in an MXD.  Is there any tool/script/methodology to do this?

I think it is a one way street where feature classes do not carry MXD usage information.  If that is the case, any other ideas on how this might be accomplished?


texts: How to zoom to selection on ArcGISDynamicMapServiceLayer?

texts:
I am using the latest ArcGIS Javascript API. All the examples and code I've found so far only deal with geometries of a FeatureLayer. How do I zoom to the visible features of a ArcGISDynamicMapServiceLayer or ArcGISTiledMapServiceLayer layer?

By visible I mean that some of the features of the layers are filtered out thanks to the layerDefinitions property.

I wish I could keep using FeatureLayers but I can't because I am running in to a wall due the sheer amount of geometry data downloaded, so I have to let the server do the rendering work and use a ArcGISDynamicMapServiceLayer to get images instead.

Thank you.


texts:
I am trying to build a raster representing the output from a logistic regression model. i keep getting syntax errors. heres the script:

-2.41-(15.97 * "Topographic_Predictors\wadem10m_annsolrad") - (9.56*("Topographic_Predictors\wadem10m_annsolrad"^2))+(.01*"Topographic_Predictors\wadem10m_slope")-(11.25"Topographic_Predictors\TPI\wadem10_tpi20")-(12.55*("Topographic_Predictors\TPI\wadem10_tpi20"^2))


texts: Is it possible to change GoogleEarth folder icon to reflect style of hidden polygon children within?

texts:
In my GoogleEarth kml file I have multiple polygon groups, each "like" group gathered under their own folder. I have hidden the individual polygons in their respective folders by "checkHideChildren" so as to reduce clutter. However, by doing so, the shared polygon style within each folder is also hidden - i.e., all one sees at the top level are plain folder icons.

Is there a way to change the respective folders icons to reflect the style/color polygon patch within whilst keeping the children hidden? I do not wish to have an external link to a custom icon. I would like to keep the kml file self-contained and utilize the existing style or icon system, as appropriate.


texts: Is there any valid way to include a Point feature in geoJSON if the coordinates are unknown?

texts: The Point I'd like to include is a marking on a historical map. Most of the items on the map have been located, so they make valid geoJSON Points. A few locations can not be identified, but I'd like to include them in the geoJSON file for completeness. As best I can tell, coordinates: [] and coordinates: null are rejected by geoJSON readers. Is there another option?


texts: Duplicate polygons into new polygons based on values in columns

texts:
I have a polygon shapefile which has polygons that have single mill name attributed to it and also some polygons which have multiple names in other columns (Mill1, Mill2,Mill3).

What I would like to do is be able to split the polygon with multiple Mill name attributes into new polygons so that they are like the below.

I am struggling to figure out a way to do this easily in ArcGIS. I am also open to options of solving this using arcpy.


texts: Creating prediction grid within certain boundary for surface interpolation in R?

texts:
I am trying to create a prediction grid within my boundary polygon.

Firstly I need to make sure that my boundary shapefile is in UTM coordinates:

 NZGD <- CRS("+init = epsg: 2193")

Load my file:

boundary <- readOGR(dsn = "C:\\Users\\gjiang\\Documents", layer = "boundary")
border <- spTransform(boundary,NZGD)
plot(border, axes=TRUE, main="Projected", cex.axis=.95)

The file has been converted and the plot looks ok.

So I start to create a prediction grid:

vals <- border@bbox
deltaLong <- as.integer((vals[1,2]-vals[1,1])+1.5)
deltaLat <- as.integer((vals[2,2]-vals[2,1])+1.5)
gridRas <- 5
gridSizeX <- deltaLong/gridRas
gridSizeY <- deltaLat/gridRas
grd <- GridTopology(vals[,1],c(gridRas,gridRas),c(gridSizeX,gridSizeY))

pts <- SpatialPoints(coordinates(grd))
pts1<- SpatialPointsDataFrame(as.data.frame(pts),data=as.data.frame(rep(1,nrow(as.data.frame(pts)))))
proj4string(pts1) <- CRS("+init=epsg:2193")

Check the results:

plot(pts)
plot(pts1)

They show me the entire grid!

It also turns out error in the following step:

pts1$border = over(pts1, border)

# Error: identicalCRS(x, y) is not TRUE

Can anyone tell what have I done wrong?


texts: ArcMap 3D profile graph gives a level line

texts:
ArcMap's Profile Graph in 3D analyst works great for drawing a line on the map and seeing the profile.  I am doing this part easily; single segments, multiple segments it works fine.

However, it is not working for me with existing features.  I follow directions from here (https://dartgeoblog.wordpress.com/2009/07/20/creating-topographic-profiles-in-arcgis/) for a vector line layer.   I have loaded gdb feature class and shapefile.  I select one (or many) features, then click 'Profile Graph'.  The result with existing features is a Profile Graph window showing a level line, elevation = 1, that is an appropriate length.  The DEM is visible underneath the entire segment, (and is not flat).  I have tried providing gdb feature classes before and after 'Add surface information' to make it 3D enabled.  I have rebooted ArcMap.  What am I doing wrong here?


texts:
ArcGIS is one of the most popular software-packages available for GIS.

What are the lesser-known but useful features of the ArcGIS? Which little/hidden tricks do you known that make your experience with this software better?


texts: Extract random coordinate points using classes (ranges) of values from a raster image

texts:
I want to use a raster to stratify an image using different texture. I have a raster image (ndvi_2019) for which I want to extract random coordinates as points from each classes. First, I want to form these create these classes I want to extract these coordinates by classing the values in the ndvi raster object as below;

class 1= <0
class 2= 0:0.1
class 3= 0.1:0.2
class 4 = >0.2

The assumption is that the classification forms the different strata which can then indicate different strata.

ndvi_2019<- (choma_clip[[4]]- choma_clip[[3]])/(choma_clip[[4]]+choma_clip[[3]])
r<- ndvi_2019

Attempt:

r<- setValues(r, rnorm(ncell(r)))
pts<- spsample(as(extent(r@extent, "SpatialPolygons"), 30, type="random")) class1<- r/r <0
class2<- r/r[r== 0:0.1]]
class3<- r/r[[r== 0.1:0.2]]
class4 <- r/r>0.2


texts: Why QgsFeature.geometry() returns an empty geometry in PyQGIS?

texts:
I am a new for QGIS, working in Python.

I add a buildings(.shp) vector layer into QGIS, just like this:

I'd like to get the geometrys of some features, then I can do analysis on the center point.

But when I use API geometry(), it returns me a empty geometry and I cannot get the center point. When I tried .hasGeometry(), it returns me a False.

My code is following:

path = '/Users/vctcn93/Desktop/export/canton.shp'

def SelectRange(px1, py1, px2, py2):

    layer = qgis.utils.iface.activeLayer()

    p1 = QgsPointXY(px1, py1)
    p2 = QgsPointXY(px2, py2)
    rec = QgsRectangle(p1,p2)

    layer.selectByRect(rec, 0)

    ids = layer.selectedFeatureIds()
    features = [QgsFeature(id) for id in ids]

    return features, rec, layer
    #retruns list<QgsFeature>,QgisRectangle,QgisVectorLayer

def SelectFea(features, rectangle):
    for a in features:
        geometry = a.geometry()#here returns me a empty geometry
        bbx = geometry.boundingBox()

        if bbx.xMaximum() >= rectangle.xMaximum() \
        or bbx.yMinimum() <= rectangle.yMinimum():
            features.remove(a)

    return features

Please help me on this problem, I tried a lot but cannot figure out.


texts: Apportioning two layers using mapbox studio classic

texts: I am using Mapbox Studio Classic to take two shapefiles--one with census tract geometries and the other with my own custom defined geometries. I want to apportion the two geometries so that I can get population data for my custom defined geometries. Is there anyway to do this using Mapbox? This is for a political science project and I don't have access to more robust programs like ESRI Arcgis.


texts:
I loaded OSM data into PostgeSQL with Osmosis&Osm2PgSQL, and see many tags are missing. I need to make it import lanes tag. How can I make it save this tag too?

edit: my solution was to copy /usr/share/osm2pgsql/default.style into a new file, edit it adding a line:

way   lanes         int4         linear

then call osm2pgsql with parameter:

osm2pgsql ... --style my.style


texts: Identifying contiguous polygons using ArcGIS Desktop?

texts:
I am trying to identify areas of contiguous conservation easements using the NCED data.  Unfortunately, the ArcGIS Polygon Neighbors tool only analyzes first-order contiguity while I'd like to analyze recursively.  This python script for QGIS worked fine but also only returns first-order contiguity.

As per the top comment under this question (@Hornbydd) I think the best approach would be to utilize the iterators in ArcGIS ModelBuilder but I'm not sure it's possible.

Any ideas on other ways to conceptualize this so that I don't have to use rely on the polygon neighbors tool?

Or is the best method really just to dig in with custom python scripting?

I'm platform agnostic but more experienced with ArcPy than python with qgis.

The dissolve method works but still keeps some areas separate that in practical terms are contiguous even if the polygons are separate... in the screenshot below the two red outlined polygons should be one contiguous area... likely an artifact of disparate data sources but it would be great to set a distance on this (0.5 miles).


texts: I want to change the name of an attibute field in QGIS 3. The duplicate referred to seems to require the plug in table manager which is no longer available. I have looked at properties and have spent much time trying to work out how to do this.


texts: Using proximity tools in ArcGIS Desktop?

texts:
I want to find the minimum distance of each point to its enclosed boundary, I have been trying to use the near tool, but the field is giving 0 as values.

This is a picture of the attribute table.


texts: How can I import data from Excel to SAGA?

texts: I have just started using SAGA and I want to plot some geographic coordinates (latitudes and longitudes) and I want to import them from Excel. Please tell me what to do from Excel to when the points are plotted.


texts: AGOL Python Api - Modify Symbology of a Feature Layer

texts:
I want to modify the JSON properties of a Feature Layer hosted online to display a certain field in symbology. I am able to access the JSON properties, even the property itself but when I call an .update_definition on it I get :

AttributeError: 'PropertyMap' instance has no attribute 'update_definition'

Here is my code :

item = gis.content.get('XXXXXXXXXXXXX')
print(item)
flc = FeatureLayerCollection.fromitem(item)
print(flc)
print(flc.manager.properties)
#print(flc.manager.properties.maxRecordCount)
#flc.manager.update_definition({'maxRecordCount':2000})#works
print(flc.manager.properties.layers[0].drawingInfo.renderer.field1)
flc.manager.properties.layers[0].drawingInfo.renderer.update_definition({'field1':'AQROUTES_2'})

What is weird to me is that I am able to modify the maxRecordCount attribute (commented out) but not the drawingInfo.

Has anyone ever had success doing this?


texts: Regex to get unique values in field using QGIS field calculator

texts:
I need to get the unique values from a field.

(?<=,|^)([^,]*)(,\1)+(?=,|$)

Works in an online test but I get invalid expression in QGIS Field Calculator. How can I rewrite this in QGIS Regex implementation or is there another maybe pythonic way to get the same result? In the end it needs to work in QGIS Modeller.

regexp_matches( "Locality_L","(?<=,|^)([^,]*)(,\1)+(?=,|$)")

Update 2
As per @she_weeds suggestion I tried with '' around the regex statement, it runs but there's no value in the table even though a sample result is shown in the field calc pane.


texts: How to tell if a GeoJSON path intersects with another feature in Leaflet?

texts:
I have an application where the user draws a path (a series of connected straight lines) and this path may not intersect any feature in a particular GeoJSON layer.

I need to check that no point along these lines intersect the GeoJSON layer, not just the endpoints.

How can I perform this check?


texts: ST_MakePoint(long, lat) HINT: No function matches the given name and argument types. You might need to add explicit type casts

texts:
I am trying to convert a normal table containing coordinates value to a geospatial data table .Which will contain  geometry column too .

For additional information I created table as

 create table latLng(
 lat varchar,
 long varchar);

  alter table latLng add column geom geometry(Point,4326);
  UPDATE latLng SET geom = ST_SetSRID(ST_MakePoint(long, lat), 4326);

ERROR:  function st_makepoint(character varying, character varying) does not
exist
LINE 1: UPDATE latLng SET geom = ST_SetSRID(ST_MakePoint(long, lat):...
                                        ^
HINT:  No function matches the given name and argument types. You might need
to add explicit type casts.


texts:
I get a lot of shapefiles with no metadata.  Grrr.  Often there are tons of table columns (fields) with no description of what the column is for or what the coded values mean.  Is there an easy way to export a table so that the domain of unique values are listed per column?  In other words, if a table has two columns, say, Color and Shape, is there a way to export a table that lists all the unique values in Color and all the unique values in Shape?  I know a quick way is to open "Select by Attribute" in ArcMap, select a field and click "Get Unique Values" to see the list of unique attributes within the field.  It would be great if there was an easy way to export a table of these unique values by column (field).

To keep it simple, I'm just interested in text or string fields.


texts: Extracting nightlight data using Google Earth Engine

texts:
I would like to extract Night-time light data (VIIRS) for a set of points in Australia (imported from csv file). The output, however, does not show the nightlight information. I am not sure if there is something wrong with the points list (the link available below). Also, is it possible to get the monthly data instead of yearly data?
https://code.earthengine.google.com/?asset=users/tronganhtrinh/points
// Load collection and geometry
var australia = points;
var collection = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG').select('avg_rad');

// Declare years of interest
var years = ee.List.sequence(2015, 2016);

// Map a function to select data within the year and apply mean reducer
var byYear = ee.FeatureCollection(years.map(function(y) {
      var image = collection.filter(ee.Filter.calendarRange(y, y, 'year'))
              .mean();
        // get the mean value for each year
      var newFC = image.reduceRegions({
                              // use reduceRegion for single geometry
                              // use reduceRegions for multiple features
        collection: australia,
        reducer: ee.Reducer.mean(),
        scale: 30,
        tileScale : 1, // increase tileScale if running out of computation or export
      }) // NOTE THAT THIS RETURNS A FEATURE COLLECTION, thus we map over the features to add date
        .map(function(feat){
          return feat.set('year', y,'system:time_start', ee.Date.fromYMD(y,1,1).millis());
        });
        // now we return the featurecollection, not the image
      return newFC;
    })
  ).flatten(); // flatten() the feature collection of feature collections

// to test code printing, increase the scale or decrease size goemetry
// print(byYear)

// export the image collection as a table (CSV)
Export.table.toDrive(byYear,"avg_rad","avg_rad","avg_rad_Australia");


texts: Is it possible to search "My Places" in Google Earth?

texts:
The search feature in Google Earth appears to only search the Google database.  Is there a way to make it search the My Places content on ones own computer?

I read that a My Places search was possible in older versions (circa 2010) using CTRL+F, but this does not appear to work in the current version, 7.1.2.2041, which was released 10/7/2013.


texts: Creating Data Frame or Group Layer using ArcPy?

texts: I am attempting to write some Python code which creates a blank data frame or group layer which can then be populated with other layers. I'm at a loss as to how to do this since I'm somewhat new to Python coding.


texts: Copy and paste a DXF Polygon from Deed Plotter to a Shapefile in ArcViewer troubles

texts: I am creating an easement layer in ArcViewer by using Deed Plotter to create the shape of the easement and then adding the dxf file as a layer in ArcViewer.  Then I start an edit session, editing the Easement.shp layer so that I can select the dxf layer, copy, and paste it into the Easement.shp.  However, recently I started to get the error: "The features could not be copied because the target layer is the wrong geometry type.  Item not found in this collection."  The Easement.shp is a polygon layer and the dxf is a polygon.  Why am I getting this error?


texts:
I've got a SPROC in SQL Server 2008 where I'm trying to calculate the distance between a point and a polygon, and I'm getting unexpected results. I am passing a parameter (@Radius) in miles to the SPROC, but the distances don't seem to be in KM. I can't figure out what they are in. My SPROC is below:

DECLARE  @Point geometry;
SET   @Point = geometry::STGeomFromText('POINT(' + CAST(@Longitude AS nvarchar(32)) + ' ' + CAST(@Latitude AS nvarchar(32)) +')', 4326);

SELECT  [FS_Geometries].[Description], ([Geometry].MakeValid()).STDistance(@Point)
FROM  [FS_Geometries]
WHERE  (([Geometry].MakeValid()).STDistance(@Point) <= (@Radius * 1609.344))
AND   [ParentId] = @CustomerId
AND   [ParentType] = 80
ORDER BY (([Geometry].MakeValid()).STDistance(@Point)) ASC

In all my other queries, I need to convert my radius from miles to meters. When I do that here, the results are not correct. I have a result set of about 750 polygons that are separated by about 25 miles of geography, so when I pass "1.0", I'd expect to get less than my 750 polygons. Unfortunately, I get the full 750 in the results.

Is the STDistance() method in SQL Server 2008 based on something other than meters? Also, is the STDistance() method in this case measuring from the edge of the polygon or the centroid of the polygon?

Thanks!

UPDATE ---

My updated conversion function...

public static double ConvertMilesToDegrees(double miles, bool isLatitude)
{
    return isLatitude
        ? miles * (((double)1) / 69)  // 1 Deg. Latitude = 69 Miles...
        : miles * (((double)1) / 61); // Avg of 69 and 53...

    // NOTES: A degree of longitude is widest at the equator at 69.172 miles (111.321)
    // and gradually shrinks to zero at the poles. At 40° north or south the distance
    // between a degree of longitude is 53 miles (85 km).
}


texts: OpenLayers - How do I change the WMS URL in Code

texts:
I want to be able to change the url for a layer in code.

The reasons is that I am implementing local caching.

When local caching is Not enabled, we use 4 servers to deliver tiles to run parallel http requests.

However, when the caching is enabled, it caches tiles for the four servers. There is no guarantee that the request for the same tile will be for the same server, so, potentially, the local cache may end up four times the size it needs to be. (and four copies of the same tile will be in the cache)

To counter this, we have set up another server to be used only for users who are using the local cache to guarantee that the cache will always return a tile if it is there (simply because the server is always the same)

I have experimented with layer.getFullRequestString( { }, url );

I am able to actually retrieve the request string if I use this syntax:

layer.getFullRequestString( { } ) (note, it errors if you do not put the {} )

but have not been able to actually cause a change in the layers fetch url. The function seems to ignore the supplied url completely.

I have tried a single url and a url array to no avail.

Openlayers Documenmtation is as follows:

getFullRequestString:function(  newParams, altUrl   ) Combine the layer’s url with its params and these newParams.

  Add the SRS parameter from projection -- this is probably more
  eloquently done via a setProjection() method, but this works for now
  and always.

  Parameters newParams  {Object} altUrl {String} Use this as the url
  instead of the layer’s url Returns {String}

Any thoughts or better ways to address this very much appreciated

Thanks


texts: Deleting rows from attribute table for large raster file using ArcGIS model builder

texts:
I'm not super familiar with raster datasets, so I may be making this too complicated. I have a large vegetation type raster that I need to convert to a shapefile. When I do raster to polygon, the file size exceeds the 2GB limit.

I can separate out the vegetation types I need in the attribute table window by turning on editing, selecting by expression, switching the selection, and then deleting the ones I don't need. I'm trying to automate this in model builder (so I can easily repeat with different veg type combinations). I can use select by attribute to get the selection I need, but I can't find any way to delete rows for a raster dataset within model builder.


texts: How is ray tracing different from ambient occlusion when hillshading with SAGA?

texts: SAGA offers both raytracing and ambient occlusion in its Analytical Hillshading module, but doesn't explain how these modes differ. Is one significantly more expensive in computation, and is the output noticeably different?


texts: How to remove default map label of places in Mapbox GL JS?

texts:
I've added Mapbox map layer in my HTML page using following script.

mapboxgl.accessToken =
    '***********************';
var map = new mapboxgl.Map({
    container: 'map',
    style: 'mapbox://styles/mapbox/light-v10',
    zoom: 6,
    center: [84.1015344, 28.3911862]

});

Everything is great but there is default map label which is cluttered and I don't want them in my analysis part. So, how can I remove those labels?

Here is the image of the generated map with too much labels in it.


texts: Losing numeric table precision during shapefile export in ArcGIS Pro

texts:
I'm trying to export a GDB file version 10.5 in ArcGIS Pro 2.0.0 to a shapefile, but every time I do, I lose precision with the unique feature IDs. For example the IDs will go from 554423216 --> 554423210.000

This taints the data b/c now each feature does not have a unique ID, plus I can't accurately join attribute tables.

Why does this happen?


texts: How to get mean coordinates of features from a point layer?

texts: I have a vector layer containing several points and I would like to generate the mean coordinates of these points, using QGIS only if possible. I'm looking for something similar to centroids, but for points. Any hint?


texts: Adding additional geometry column in PostGIS?

texts:
I'm importing many sets of geodata into PostGIS, and they have different SRID's. (Some have EPSG:3857, some EPSG:4326, some something else).

I'd like to create an additional geometry column, eg. the_geom_mercator with SRID EPSG:3857, and also keep the original geom column in whatever SRID it came in.

How can I do this with a PostGIS function?


texts: Change LIDAR classification of points inside polygons

texts:
I'm trying to create a pdal pipeline to adjust the classification values for all points inside polygons, but I can't seem to get it to work. The filters.range works correctly because the output file contains the filtered values. It's the filters.overlay that doesn't seem to be doing anything. I have an attribute called CLS on all the polygons in the shapefile.

    {
        "type": "readers.las",
        "filename": "area_small.las"
    },
    {
        "type":"filters.range",
        "limits":"Classification[0:0], Classification[3:6]"
    },
    {
        "type":"filters.overlay",
        "dimension" : "Classification",
        "datasource": "buildings_small.shp",
        "column": "CLS"
    },
    {
        "type" : "writers.las",
        "filename" : "output.las"
    }

I have also tried to use a VRT file

<OGRVRTDataSource>
    <OGRVRTWarpedLayer>
        <OGRVRTLayer name="buildings_small">
            <SrcDataSource>buildings_small.shp</SrcDataSource>
            <LayerSRS>EPSG:4326</LayerSRS>
        </OGRVRTLayer>
        <TargetSRS>EPSG:4326</TargetSRS>
    </OGRVRTWarpedLayer>
</OGRVRTDataSource>

That doesn't work either. The las file does have a different projection, but they two files line up perfectly in GlobalMapper

Is there another way to do this? Using a filters.assign changes the Classification value, but I can't figure out how to make it only do that for points inside the polygons.


texts: Creating spatial Reference to ASCII file

texts:
I have produced ASCII files in MaxEnt. These files are able to be loaded into Arcmap10.1 and display as an image, but I cannot interpret them or create them into a polygon for further analysis. I have read some forums that say I need to add spatial reference and calculate statistics (http://gis.utah.gov/wp-content/uploads/ASCII_Elevation_Data.doc) but when I do this to the original file and then import it will not create the resulting raster into a polygon, stating that "The environmental window is ignored"

Ideally, I need to add all the different environmental considerations to the ASCII file and then import it as a raster. This will allow me to analyse it and then create a polygon of the extent of the layer i am looking at.

Can anyone suggest a series of steps to overcome this and produce a raster and then a polygon?


texts: How to make a series of time per class using a feature collection?

texts:
I am working with severity of the fires, for this I vectored the severity sections and exported them to my asset, then I tried to use my feature collection to make a series of time per severity section (low, medium, high), however, I only throws a result
Could you help me, I do not know if the error is when I export or when I call the feature collection

// Run reduceToVectors per class by masking all other classes.
var classes = ee.List([0, 1, 2, 3])
  .map(function(n) {
    var classImage = severidad.eq(ee.Number(n));
    var vectors = classImage.updateMask(classImage)
      .reduceToVectors({
        reducer: ee.Reducer.countEvery(),
        geometry: incendio,
        scale: 30,
       // eightConnected: false,
        maxPixels: 1e8})
      .geometry();
    return ee.Feature(vectors, {"class": n,"size":vectors.area(1)});
  });
var result = ee.FeatureCollection(classes);
Map.addLayer(result);
Map.centerObject(result)
print(classes);

// Export an ee.FeatureCollection as an Earth Engine asset.
Export.table.toAsset({
  collection:  result.filterMetadata('size', 'greater_than',0),
  description:'cochamo',
  assetId: 'severidad_Cochamo',
});

Link first part  here
https://code.earthengine.google.com/5ea068880b269f30b16afcb5e618ced7

Time series.
I try to make the time series with the vectors that I export in part I, but it does not work for me. I do not know how to select the three features that make up the feature collection, which is what I have to do with the time series.

var NDVITimeSeries = ui.Chart.image.seriesByRegion({
  imageCollection:NDVIbyyear ,
  regions:  severidad,
  reducer: ee.Reducer.mean(),
  band: 'NDVI',
  scale: 200,
  xProperty: 'year',
  seriesProperty: 'label'
});
NDVITimeSeries.setChartType('ScatterChart');
NDVITimeSeries.setOptions({
  title: 'NDVI',
  vAxis: {
    title: 'NDVI'
  },
  lineWidth: 1,
  pointSize: 4,
  series: {
            0: {color: 'FF0000'}, // severidad baja
            1: {color: '00FF00'}, // severidad media
            2: {color: '0000FF'}  // severidad alta

    }});

print(NDVITimeSeries);

link second part here
https://code.earthengine.google.com/3f54352b1bbe46a5982f1cf49e79c3f8


texts: I had set up my local cartodb instance in server. I had searched specifically for geocoding API in cartodb, but couldnot find any. I would like to know is there any API for geocoding in cartodb?


texts:
When using the QGIS SpreadSheet Layers plugin, how do I update the SQLite query without closing and reopening QGIS?  I cannot see a "watch file" equivalent.

The plugin creates and uses a *.vrt file of the same name.  This file lists each of the columns and has a separate SQL query for each column imported from the spreadsheet.

So a more general question is: How to refresh or update a vrt layer in Qgis?


texts: Problem accessing EVLR slot in LASheader with lidR version 3.1.2

texts:
I'm using lidR 3.1.2. This problem has been occurring since very recently with the following error
Error in eval(call("@", object, slot)) :
  no slot of name "EVLR" for this object of class "LASheader"`

The error occurs when I try to view the header, or with writeLAS().
On another computer that initially had version 3.1.1, there was no problem dealing with the exact same .Rdata containing the object. I updated it to 3.1.2 and the error occurred there as well.
The .Rdata file was created earlier containing an object which is a list of point cloud subsets. The subsets were created using lidR functions.
Furthermore, I reinstalled the previous version and there was no problem viewing the header.


texts: Projection system for India or Indian subcontinent with metric system

texts:
I am creating a buffer around my geometries using QGIS.
The unit of the buffer is shown in degrees because the projection system is WGS84.

I need to change the projection system so that unit for buffer comes in metric system (I want to add 10-meter buffer).

Which projection system should I use? My shapefile contains geometries are from India.


texts: Writing a script that identifies polygons within a certain distance

texts:
To start with I have no scripting, coding or programming knowledge.

I have access to Mapinfo and QGIS. The work space will contain point and polygon data for many different things such as Tree Preservation Orders (Points), Charging Zones (Polygons), etc.

The user wants to know if a point they manually enter is within set distances of the aforementioned points and polygons e.g. Is Co-ordinate within 500m of a Tree Preservation Order (and print results); is Co-ordinate within a Charging Zone (and print results).

Is this something I can achieve in QGIS Process Modeler, or MapBasic? Overall there will be about 10 layers to query, and the user queries many co-ordinates (addresses/buildings) a day.


texts: Registering multiple servers using ESRi IdentityManager.registerToken

texts:
I have multiple secured service that needs to be added to my map for which I need to register multiple servers using ESRi IdentityManager.registerToken method.
  IdentityManager.registerToken({
        'server': "https://testurl1",
        'token': "abcasdasasd"
       });

       IdentityManager.registerToken({
        'server': "https://testurl2",
        'token': "MV5VF5SeK9pigO5094rC9MMEyY39SnmB"
       });

The first IdentityManager.registerToken for the testurl1 is able to register the token successfully but the IdentityManager.registerToken for the testurl2 throws the following error:
Cannot read property '3' of null at new k (init.js:457) at b.hasSameOrigin (init.js:731) at init.js:731 at Object.some (init.js:70) at Object.b.canUseXhr (init.js:731) at Object. (init.js:1999) at Object.forEach (init.js:71) at Object.registerServers (init.js:1998) at Object.registerToken (init.js:1999) at Object.onReceiveData (Widget.js?wab_dv=2.10:48) "TypeError: Cannot read property '3' of null at new k (https://js.arcgis.com/3.31/init.js:457:489)
How can I fix this issue?


texts: Starting ArcGIS Server cached MapService without database connections

texts:
I have a cached ArcGIS 10.0 MapService that I need to start but the DB which was used as the datasource for creating the tiles is no longer available.
Is there a way to start the service without datasources?


texts: Algorithm for finding population for a given center point and radius in US

texts:
Here is what I have

The population data aggregated by zip code given by US Census.

Here is what I am doing right now:

Reverse geo code the center point to find the zip code of the center
point.
Associate a lat-long to each zip-code(by reverse geocoding, I know
this will be error prone).
Find the lat-longs(zip-codes) which are inside the circle, using
Geokit-Rails plugin.
Finding the sum of the data for those zip-codes(zip-code=zipcode of the center point or distance(between lat-long of the zip code and center point) < radius) and showing.

Am I doing this right? Is there methods efficient than this to find the population? Is there any algorithm already available?

I am using geokit-rails plugin for reverse geocoding and for finding the zip-codes(with associated lat-longs) within the circle.

The problem with this setup is that the zipcodes cannot be actually represented by the lat-long,it is an area.So, I am not sure the lat-long pair given by Google Map API is actually lies in the center of the zip code or something like that. Also,for radius 0 to almost 2.4 miles radius,the population remains same with this setup. That is, even for a radius of 0.1 miles, it will show the whole population of the zip code.I am looking for a way to approximate this.

I have asked same question in SO. Please have a look at this question also.


texts:
I had a pure Python 3 configured for myself and recently installed QGIS (just QGIS, not OSGeo4W).
I need to keep my initial Python 3 but also add qgis.core module.
Can it be done with no harm to my previously installed Python?
I saw that installing OSGeo4W should help but anyway there is some paths need to set in the user environment settings.


texts: Importing Geopandas gives OSError: Could not find libspatialindex_c library file

texts:
I would like to begin with geopandas but I have an issue when I import it.
Terminal return me this message :

>>>import geopandas
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/geopandas/__init__.py", line 1, in <module>
    from geopandas.geoseries import GeoSeries
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/geopandas/geoseries.py", line 12, in <module>
    from geopandas.base import GeoPandasBase, _series_unary_op, _CoordinateIndexer
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/geopandas/base.py", line 14, in <module>
    from rtree.core import RTreeError
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/rtree/__init__.py", line 1, in <module>
    from .index import Rtree
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/rtree/index.py", line 5, in <module>
    from . import core
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/rtree/core.py", line 125, in <module>
    raise OSError("Could not find libspatialindex_c library file")
OSError: Could not find libspatialindex_c library file

I tried to install spatialindex with pip, brew or directly with github source, and to install again Rtree but it didn't work.

Can somebody help me?

I am on macOS


texts: Generate PNG/JPG maps from multiple ASCI/GeoTIFF in Python/GDAL

texts:
I need to generate maps (saved in PNG, JPG or PDF) from raster files (for each file - one map) which meets following criteria;
* each map has exactly the same color pallette, so that it will be easy to compare them,
* I wish to add any label (for example source file name) at the top of each map,
* (optionally) I wish to include legend (color pallette)
* (I will be in heaven) if there will be coordinates in map included.

It is quite easy to save GeoTIFF as PNG (gdal_translate -of PNG), but I can't find a solution to set mentioned options. Is it possible to generate such maps from Python/GDAL? How? Should I call from Python QGIS map composer or there's another solution?

Edit

I started to write my Python script to resolve my issue.
I assume that:

I have opened a Composer View where I have all elements already added
For each raster file, I need to:

remove old raster file
add new raster file
(if needed) refresh map in composer view. Legend and all other elements should probably leave the same
export as PDF/PNG/any other graphics format.

from PyQt4.QtGui import QImage, QPainter
    from PyQt4.Qt import QSize
    from PyQt4.QtCore import QRectF

def StringToRaster(raster):
    # This function loads raster from provided path http://gis.stackexchange.com/questions/144058/how-to-load-a-raster-layer-using-pyqgis
    if isinstance(raster,basestring):
        fileInfo = QFileInfo(raster)
        baseName = fileInfo.baseName()
        path = fileInfo.filePath()
        if (baseName and path):
            raster = QgsRasterLayer(path, baseName)
            if not raster.isValid():
                print "Layer failed to load!"
                return
        else:
            print "Unable to read basename and file path - Your string is probably invalid"
            return
    return raster

def GenerateRasterMaps(rasterList,rasterStyle):
    # I wish to use exactly the same composer view - only raster file should be changed

    myComposer=iface.activeComposers()[0]

    for rasterFile in rasterList:
        if (newRaster):
            QgsMapLayerRegistry.instance().removeMapLayer(newRaster.id())
        newRaster=StringToRaster(rasterFile)
        newRaster.loadNamedStyle(rasterStyle)
        QgsMapLayerRegistry.instance().addMapLayer(newRaster)

        # Now I should only refresh current composer view and export map as a PDF. Unfortunately, after few hours on reading (and testing): http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/composer.html and QGIS API, I still have no effect. Why it is so hard to export current Composer View to PDF/PNG. Why it is not simply one command?)


texts: `Process finished with exit code 139 (interrupted by signal 11: SIGSEGV)` after `gdal.RasterizeLayer`

texts:
I get this error when I run this:

vector_layer = "mask.shp"
raster_layer = "raster.tif"
target_layer = "mask.tif"

# open the raster layer and get its relevant properties
raster_ds = gdal.Open(raster_layer, gdal.GA_ReadOnly)
xSize = raster_ds.RasterXSize
ySize = raster_ds.RasterYSize
geotransform = raster_ds.GetGeoTransform()
projection = raster_ds.GetProjection()

# create the target layer (1 band)
driver = gdal.GetDriverByName('GTiff')
target_ds = driver.Create(target_layer, xSize, ySize, bands=1, eType=gdal.GDT_Byte)

target_ds.SetGeoTransform(geotransform)
target_ds.SetProjection(projection)

# rasterize the vector layer into the target one
ds = gdal.RasterizeLayer(target_ds, [1], ogr.Open(vector_layer).GetLayer(), burn_values=[1])

target_ds = None

It happens when I call target_ds = None after RasterizeLayer. If I don't run RasterizeLayer it does not fail. So what I am trying to do is create a mask of the polygons in mask.shp using the projection and geotransform information from raster.tif.

Note I am using gdal version 1.10.1.


texts: Holes between imagery tiles projected using gdalwarp

texts:
I am having the following issue when reprojecting imagery tiles using gdalwarp.
These are 4096x4096 imagery tiles, which I'm gdalwarp'ing from WGS84 to a State Plane projection. When reprojected, the tiles, naturally, overlap, and so have added an alpha channel which takes care of it nicely.
However, when viewed together, after reprojection, the tile edges don't quite match up. Note also that this is NOT an issue or artifact with the alpha channel---I can use nodata instead of alpha, or just leave the area black, and the edges remain broken, with holes, just the same. Below is an example of a 1:80 scale crop of the seam between two tiles after projection. It maintains this pattern of holes along the entire seam between two tiles.
The basic command I'm using looks something like this:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:0000 /path/to/source.jp2 /path/to/dest.tif

I've exhausted most of the options I can think of from gdalwarp ( http://www.gdal.org/gdalwarp.html ) and the warp options ( http://www.gdal.org/structGDALWarpOptions.html#a0ed77f9917bb96c7a9aabd73d4d06e08 ), including different resampling methods, increasing the SAMPLE_STEPS density, and playing with additional SOURCE_EXTRA pixels. All with no change whatsoever.
I realize creating a mosaic first, reprojecting, and then retiling would be one way around this issue, but I'm dealing with tens of thousands of these ~16 megapixel tiles, so as far as I can tell, a mosaic is not feasible.
Any ideas on an approach to deal with these edge inaccuracies? Thanks!


texts: Combined WMS request issue (one layer has black background)

texts:
We are using MATLAB to request a combined WMS from GeoServer. The request is auto-generated as this:

http://geo-pprod.test.cresis.ku.edu/geoserver/arctic/wms?SERVICE=WMS&LAYERS=greenland_natural_90m,rds_line_paths&CRS=EPSG:3413&FORMAT=image/geotiff8&HEIGHT=527&BGCOLOR=0xFFFFFF&REQUEST=GetMap&WIDTH=620&BBOX=-1465312.4055326465,-3439603.5364015717,1669546.4179967656,-774973.5364015717&STYLES=,&VERSION=1.3.0

It's requesting a GeoTIFF store with a linestring (PostGIS store) both published in EPSG:3413.

greenland_natural_90m: Greenland GeoTIFF (declared native as EPSG3413)

rds_line_paths: PostGIS linePath geometry table (re-projected from WGS84 to EPSG3413)

The result from this request looks like this:

http://i.imgur.com/FSC6fcQ.png

However when I load the same two layers in our OpenLayers web client I get the following (correct):

http://i.imgur.com/rsr8XTF.png

Note that we have tried every image format and get the same results. We have also tested on multiple versions of GeoServer. Any help would be appreciated. Unfortunatly our server in not public yet so please keep that in mind if additional information is needed.

What we have tried:

View in OpenLayers (works)
Get each layer individually (works) **We really need it combined
Get layers (combined/individual) in multiple formats (Works individual, fails combined in MATLAB)
Revert to old versions of GeoServer (this has worked before, but continues to fail now)

Notes:

We are running Geoserver 2.3.2 WAR served via Tomcat6 through Apache on an RHEL linux server.


texts: How to use a IF statement to change colors for layers using PyQgis?

texts:
How to change the colors according the values in 'C_index' field?
I tried a FOR loop, in order to read all values in this field, and an IF STATEMENT to set all 'red' and 'blue' values to the same colors. But it keeps returning me just one color, the first one in the loop, as shown in the image below.
Here's the entire code:
iface.addVectorLayer('E:/Bkp/Programacao/Python/Declividade_linhas/Linhas_pontos.shp', '', 'ogr')
layer = iface.activeLayer()

exp = 'case\r\nwhen  \"MDT\" > \"Z1\" then angle_at_vertex($geometry, 1)- 0\r\nwhen  \"MDT\" < \"Z1\" then angle_at_vertex($geometry, 1)- 180\r\nend'

simpleMarkerSymbolLayer = QgsSimpleMarkerSymbolLayer()
simpleMarkerSymbolLayer.setShape(QgsSimpleMarkerSymbolLayer.Arrow)
property = QgsProperty()
property.setExpressionString(exp)
property.setActive(True)
simpleMarkerSymbolLayer.setDataDefinedProperty(QgsSymbolLayer.PropertyAngle,
property)
markerSymbol = QgsMarkerSymbol([simpleMarkerSymbolLayer])
markerSymbol.setSize(7)

markerLineSymbolLayer = QgsMarkerLineSymbolLayer()
markerLineSymbolLayer.setSubSymbol(markerSymbol)
markerLineSymbolLayer.setPlacement(4)

#set colors
for f in layer.getFeatures():
    if f['C_index'] == 'blue':
        markerLineSymbolLayer.setColor(QColor('blue'))
    else:
        markerLineSymbolLayer.setColor(QColor('red'))

lineSymbol = QgsLineSymbol([markerLineSymbolLayer])
layer.renderer().setSymbol(lineSymbol)

layer.triggerRepaint()


texts:
I am attempting to convert a GeoJSON Feature Collection into multiple shapefiles. The shapefile a feature ends up in depends on a specific property of the GeoJSON, layer_id.

My desired output is a single zipped file, containing all the different shapefile folders.

So far I have converted GeoJSON with uniform layer_id to a single shapefile, using the ogr2ogr npm module. The output of the ogr2ogr command is a zipped shapefile.

I would like to either:

Learn how to designate a property for ogr2ogr to "group by" when
creating a shapefile, creating multiple shapefiles
Learn how to create unzipped shapefiles with ogr2ogr, which I can then compress into a single zip


texts: How to effectively merge SRTM 1 arc data?

texts:
We have around 1800 1 arc-second SRTM tiles (Geotiff). How to effectively merge all the tiles into one single file ? We are thinking about splitting 100 files per folder and merge them together using QGIS merge. We don't have to create VRT as the client insists on a single tile.

My question is what are the parameters we have to tweak in QGIS Merge option ? There are few parameters in "Merge raster layers (through Processing tool)-Preferred data storage, Interpolation, Overlapping cells. Read somewhere, in case of SRTM tiles the first pixel line will overlap with the adjacent one. So do we have to use "First value in order of grid list" in case of overlapping cells.

Any kind of suggestions are welcomed.

Edit:

Finally i've managed to merge 2600 files. Divided the files into five folders and created VRT for each folder of 515 files. Then converted VRT to Geotiff using gdal_translate command in OSGEO shell (with DEFLATE, PREDICTOR = 2 parameters). Original file size was around 6.25GB. But the compressed output was 3.25GB each.

Then i again created VRT for those five files and did the same translation.


texts: Using Google Maps Tile Service with OpenLayers 5

texts:
How to pass Google Maps API key to access tile service using OpenLayers 5.

Some google searches suggested me to set the url of tile layer as

http://mt.google.com/vt/lyrs=m&x={x}&y={y}&z={z}

Is this a legal way to access Google tile server data directly without an API key?

As per this forum, its a violation of Google Maps TOS.

If this is a violation, what should I need to do to meet their TOS.

I need to overlay some WMS data hosted on GeoServer over Google Map.
I can do so using Google Maps tile server, but without using a Google API.


texts: Writing raster field values into attribute table using ArcPy?

texts:
Is it possible to use ArcPy (10 or 10.1 ) to write raster field types into a featureclass in either ArcSDE or a file geodatabase?.

I am hoping to use arcpy as means of automating the process of adding rasters as an attribute to a feature class.


texts: Find number of weighted point clusters & centre of clusters

texts:
I have a points layer with a frequency field indicating a population at that point. There is also a field indicating the broader area by which I need to analysis broken down into. I am expecting the points to be mostly grouped in a number of clusters. I need a way of determining how many clusters there are, probably by density, and to return a point for the centre of each.

I've made a heatmap in the display section of the points layer by setting a radius and weighting by the frequency field and in this created three clusters in my trial area, that look right, but I don't see how to get a coordinate for them. Any suggestions of what tool(s) to look at.

Identifying clusters in vector point data using QGIS?  is basically the same, but as none of the answers explain how I can do weighted cluster analysis in QGIS.  I don't know postGIS or R, and the plugin manageR that was cited no longer seems to exist.


texts: How can I constrain this filtering even more?

texts:
I'm using QGIS and I am trying to filter a dataset using the advanced filter. My dataset contains earthquakes in the U.S., Canada and Mexico, however I would like to exclude Canada & Mexico for now.

The problem arises when I filter using:

  "place" LIKE '%Mexico%'

as this now includes Mexico AND New Mexico.

I've been trying to experiment with SQL commands but I've had no luck so far.

EDIT: Added screenshot of part of my attribute table, the column labeled "place" is what I'm trying to filter by.

"earthquakes_2005.csv" dataset


texts:
I like to convert a georeferenced tif-file into a dwg-drawing.
At first I convert the raster into a polygon. The result is Grey-coloured picture
When I "Export to CAD" the polygon it is still Grey.
How can I secure the original color scale?

I am using IBSCAD and ArcGIS.


texts: QGIS Calculate Attribute Value from sum of dependent Point features

texts:
I have two layers one with line geometry and other with point geometry. Features from both the layers are perfectly snapped.
In the line layer, I want to calculate the field which is depend on features in point layer. It is like calculating value on the basis of downstream trace from line end point to the further till last end point.
The line attribute value is depends on the accumulating sum value of point feature attribute towards downstream. (From line end point to actual last line end point. Flow direction added.)
Like in the image, the selected line should get the value from downwards point (as per the flow direction) attribute as a sum of their values (10+9+0+7+8). I am not getting any idea about how can i make it.
Please suggest.


texts:
I am looking for information regarding error reporting with WMTS. I have consulted the official documentation at http://www.opengeospatial.org/standards/wmts and I found several chapters discussing exceptions. All the information is there, but it is a large document and I am having trouble quickly determining where to look.

I also found this website https://github.com/nasa-gibs/onearth/wiki/WMTS-Response-Codes which gives a clear overview of how an exception should be formatted.

I have two questions. Is the second link the correct way of doing it according to the standard? And does anyone have any more resources explaining WMTS error reporting (Google so far has not helped me).


texts: How to customize ArcGIS Silverlight API 2.4 legend?

texts:
I am using Esri's ArcGIS Silverlight API 2.4 and below is my legend.

My map has one feature layer and one dynamic service layer with nested layers like below

Feature Layer A
Group Layer
Group A
Layer A1
Layer A2
Layer A3
Group B
Group B1
Layer B1
Layer B2
Group B2
Layer B2
.....

'Feature Layer A' is the same as 'Layer A1' and I want to bring this 'Feature Layer A' to position of Layer A1. But this layer should be remain as feature layer only.

How can I do this? i.e feature layer inside the dynamic service layer group?


texts: Having trouble with 'and' 'or' statements and negative numbers

texts:
I have set up a program that calculates the relationship between a point and the first point of a polyline. The aim of the program is to run this calculation and then copy the result to a certain field in the attribute table. However, if the cell/field is already occupied it copies to the next field in line and so on.

This runs pretty smoothly in most cases. However, if the result/calculation happens to return a negative number all hell breaks loose. When the negative number is returned the loop registers the result in the correct up/down link column but seems to ignore the command to check if the cell is already occupied. This leads to other valid data being overwritten.

I'm not sure whether this is a result of the way I have set up my 'and'/'or' operators or whether I modify the way I use negative numbers in python.

Below is an example of the trouble I have been having. If anyone has any ideas please let me know.

q_fid = 1
calculation = -30
linkup = "EDRN21"
linkup2 = "EDRN03"
linkup3 = " "
linkdown = "EDRN10"
linkdown2 = " "
linkdown3 = " "
link_id = "EDRN44"

arcpy.MakeFeatureLayer_management("uk_cities", "node_feature")
arcpy.SelectLayerByAttribute_management ("node_feature", "NEW_SELECTION", '\"FID\" =  {0}'.format(q_fid))

rows = arcpy.UpdateCursor("node_feature")
for row in rows:
    if calculation < -0.001 or calculation > 0.001 and linkup == " ":
        row.setValue("Up_LinkA", link_id)
    elif calculation < -0.001 or calculation > 0.001 and linkup2 == " ":
        row.setValue("Up_LinkB", link_id)
    elif calculation < -0.001 or calculation > 0.001 and linkup3 == " ":
        row.setValue("Up_LinkC", link_id)
    elif calculation > -0.001 and calculation < 0.001 and linkdown == " ":
        row.setValue("Down_LinkA", link_id)
    elif calculation > -0.001 and calculation < 0.001 and linkdown2 == " ":
        row.setValue("Down_LinkB", link_id)
    elif calculation > -0.001 and calculation < 0.001 and linkdown3 == " ":
        row.setValue("Down_LinkC", link_id)
    rows.updateRow(row)


texts:
I am trying to load WMS layer from geoserver in QGIS.
The Error in QGIS is :
Failed to Download Capabilities
Download of Capabilities failed: Connection Refused

So I went and found this question:
Failure to connect QGIS with GeoServer WMS?

but the solution in the first answer when I go to the following link:

http://localhost:8080/geoserver/wms?SERVICE=WMS&REQUEST=GetCapabilities

I get an XML Parsing Error, how can I fix this error?

  XML Parsing Error: not well-formed Location:
  http://localhost:8080/geoserver/wms?SERVICE=WMS&REQUEST=GetCapabilities
  Line Number 6244, Column 50:
  -34.8977165222168http://localhost:8080/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd">

Ok Update:
Some of the sld styles did not have workspaces assigned to them, once they were sorted I still cannot load th WMS in QGIS.

When I load the localhost get capabilities link same as above ^^^

I no longer get the same error message(progress?)

But now the error states:
This XML file does not appear to have any style information associated with it. The document tree is shown below.


texts: How to change the border of CircleMarker in Leaflet?

texts:
I'm using  Leaflet 1.3.4 and I want to change the border of CircleMarker (make it bigger).
I created it in this way:

      var geojsonMarkerOptions = {
        radius: 25,
        fillColor: "#24D238",// "#28ea3f",//"#0163FF",
        color: "#A9F6B2", //"#0163FF",
        weight: 2,
        opacity: 1,
        fillOpacity: 1,
        // className: 'marker-cluster'

      };

      geoJson = new L.geoJson(buildingPoints, {
        pointToLayer: function(feature, latlng) {
          return new L.circleMarker([latlng.lat, latlng.lng], geojsonMarkerOptions).bindTooltip('53', { className: 'text', direction: 'center', permanent: true, interactive: true });//{color:"#0163FF", radius: feature.properties.radius});
        },
        onEachFeature: function(feature, layer) {
        }
      }).addTo(map)

I know how to do it via DivIcon, but I need to use CircleMarker.
I can't find any border param in leaflet.css corresponds to CircleMarker.

I've tried specifying className in options, but it doesn't help (even if I use SVG renderer).


texts:
I am very new to WPS and currently playing around with the WPS Request-Builder.
I would like to receive the top x highest or lowest features by a given attribute.
I have been trying to achieve this by vec:Query and adding filters. I imagined a combination of sortBy and maxfeatures - but these seem to be strict WMS vendor parameters and not available for this kind of WPS requests.


texts: Interpolating raster from points from a vertical datum conversion point feature class

texts:
I have a shapefile with points distributed uniformly representing the difference between two vertical datums. These are spaced in 4 corners of a 2950 x 3700 square.

I would like to create a continuous raster at a lower resolution by interpolating the nearest point values, so I can convert between datums more accurately.

What would be the best raster interpolation tool for this job and what parameters should I use?


texts: Resolving Arcpy "GetObject not point" Error

texts:
I'm new to Python (v 2.6.5) and am trying to write an ArcMap (10) script that reads the vertex coordinates of a polygon shapefile which includes a small hole (doughnut polygon).  I've been working from the ESRI Help page on reading geometries.  The script I've developed outputs coordinate information until it hits coordinate 1025 of 1088 at which point it fails with the error "Array: GetObject not point".  I suspect this is the 'null point object' returned as the separator between my outer and inner rings (described in the Help), but I'm stumped on how to account for this in a way that my script doesn't fail (I thought the "if pnt:" line would handle that...).  Here's my script - thanks for any guidance!

import arcpy
infc = arcpy.GetParameterAsText(0)
desc = arcpy.Describe(infc)
shapefieldname = desc.ShapeFieldName
rows = arcpy.SearchCursor(infc)
for row in rows:
    feat = row.getValue(shapefieldname)
    arcpy.AddMessage("Feature %i:" % row.getValue(desc.OIDFieldName))
    partnum = 0
    for part in feat:
        pointnum = 0
        arcpy.AddMessage("Part %i:" % partnum)
        for pnt in feat.getPart(partnum):
            if pnt:
                arcpy.AddMessage("Point number " + str(pointnum) + ": " + str(pnt.X)+ ", " + str(pnt.Y))
                pointnum += 1
            else:
                arcpy.AddMessage("Interior Ring:")
        partnum += 1


texts:
The Python documentation for OGR lists a Polygonize() method belonging to the Geometry class, as does the C documentation.  The docstring suggests it can be called on a geometry object of type "MultiLinestring", however doing so results in an AttributeError.  Other methods such as Buffer() work fine.

Did they move or remove this in GDAL 2.x?  Lib/site-packages/osgeo/ogr.py does not contain a Polygonize() function.  It does say "OGR 1.9.0" in the documentation, but I would assume this normally means OGR >= 1.9.0.

Currently running gdal version 2.2.4.

Code sample:

multiline = ogr.Geometry(ogr.wkbMultiLineString)
multiline.AddGeometry(line)

polygon = multiline.Polygonize()

Results in:

AttributeError: type object 'object' has no attribute '__getattr__'


texts: What is the right syntax to the decision tree using the Con statement or any other method to classify?

texts:
Can you please show me how to construct the decision tree for the attached land cover classification.
I assume it is nested Con statement, but I am confused on how to write the statement.


texts: Programmatic authentication to ArcGIS Server secured layers via RESTful API

texts:
I have an ArcGIS 10.1 Server instance exposing secured mapservices on the Internet.
My need is to code down a client application (which I am currently building using the 3.3 version of the ArcGIS Javascript API) enabling the user to view those secured webservices:

I think this online ESRI example is a good start.

My will is not to have the user prompted for EACH mapservice authentication because I already know that ALL of the mapservices belong to her/him and therefore are accessible by her/him under the same username&password.
In my idea, the prompt for credentials should appear just ONCE and therefore the JS code should feed credentials to each of the mapservice via some sort of RESTful login call.
To me it doesn't seem like the ArcGIS Server REST API is providing such a call...maybe I am wrong.

Is therefore this "RESTful" way to login to secured mapservices feasible with ArcGIS Server (making it possible to access secured services programmatically)?
If so, can any of you provide examples or links to web resources explaining this?


texts: How to show aggregated values in a CARTO BUILDER popup?

texts:
Is it possible to have some column calculations displayed in popups?

Given this example map:

https://atomotic.carto.com/builder/8238cb36-0f15-4fa4-9c80-cdce7c2bd040/embed

I would like to:

display a popup on a city that shows a value calculated using sum on a
column,
access the value selected on category widget in order to
compose an URL.


texts: Web Mapping Services from Europe besides OSM

texts:
Wondering if there are any other Web Mapping Services based in Europe or even Asia?
Fed up with the Bings and Googles and even OSMs of the WMS world.
This is 2012 not 2002, shouldn't we have a much better FREE view of the world by now?


texts:
I have a question in order to developing a web app using dynamic maps with JSON files. I'm using Visual Studio 2015 with MVC 5.

At this moment i have developed complete app but I'm unable to draw points (with Lat and Lng) in a Leaflet map from a JSON File.

My JSON file is similar to:

[{
"latitud": 43.526523590087891,
"longitud": -5.6150951385498047
}, {
"latitud": 43.511680603027344,
"longitud": -5.6671133041381836
},
 .
 //More Lat-Lng Points
 .
{
"latitud": 43.526451110839844,
"longitud": -5.6140098571777344
}]

Currently, I'm using a example of heatmap but is imposible for me represent markers when i try to get coordinates from my JSON file.

Any help?


texts: PyQGIS 3.4: otb:unsupervisedkmeansimageclassification -> AttributeError: module 'processing' has no attribute 'runalg'

texts:
I'm trying to classify a raster with PyQGIS:

import processing
processing.runalg("otb:unsupervisedkmeansimageclassification","FalseColor.tif",128,None,100,3,1000,0.95,"class.tif",None)
The python console indicate :

  AttributeError: module 'processing' has no attribute 'runalg'

It is a little bit surprising because I've followed the QGIS documentation. https://docs.qgis.org/testing/en/docs/user_manual/processing_algs/otb/learning.html#unsupervised-kmeans-image-classification

I've seen this post: QGIS 3.0 Error when calling processing.runalg()

I've tried with processing.run:
processing.run("otb:unsupervisedkmeansimageclassification","FalseColor.tif",128,None,100,3,1000,0.95,"class.tif",None)

The python console indicate :

  TypeError: run() takes from 2 to 5 positional arguments but 10 were given

Is there an updated documentation about the algorithm otb:unsupervisedkmeansimageclassification?


texts: Similarity between two sets of lat/long points

texts:
I am trying to measure the similarity between two trajectories (two sets of lat/long points) which are the outputs from two APIs. I am trying to measure the difference between the routes given by these two APIs and give an estimate of how much percentage of routes match. I have seen about Frechet Distance and Hausdorff Distance but couldn't understand the details. My understanding is Frechet Distance gives even better results than Hausdorff Distance and my idea is to take the two polylines (outputs from two APIs) and take their Frechet Distance. But how can I get the percentage of match between these two routes. Does
densifyFrac parameter in these methods can give any such information. Or are there any better approaches for this problem.


texts: Transform scale of polygon pattern fills with Carto/Mapnik

texts:
In CartoCSS you can specify an SVG pattern fill for a polygon as so:

#layer {
    polygon-pattern-file: url(images/fill.svg);
}

However, the size of the fill is too big at small scales and too small at large scales, and I would like to apply a custom transform to the fill to ensure that it scales properly. Looking at the documentation for polygon-pattern it appears that you can only change the scaling of the polygon being filled, not the actual fill itself.

Is there a way to apply a custom scale to SVG pattern fills? Essentially, it is applying an SVG transform=scale() to the main group of the SVG and changing the root element's height and width.


texts: AGS Tiled Map Services and the .jpg file extension using js api

texts:
I have an old app that is still using esri js api 1.4 (i know, it's dated) but that app works and it can pull tiles from a 10.0 tiled map services. When it does that it calls for the tiles at the url with the \000\111\222 type format. It also appends a '.jpg' to the end of the request url. That works fine - i can even knock off the '.jpg' and still get the tile image.
But now, my tiled service is being updated to 10.2 and I'm having some issues. When the request goes out with the .jpg I am no longer able to get the image - i get a 400 invalid url.
If i lop off the .jpg then i get the image.
The problem is that i'm not the one adding the .jpg to the request. That is the old api interacting with the ags service to produce that call. It's a black box to me. I load the layer and the js api does it's thing.

Am i missing something or is my js api just too old to keep working?


texts: Using QGIS 2.14 to creating a evenly distributed grid of points, spaced 500m apart, that takes terrain into account

texts:
I am using QGIS 2.14 (Essen Version).

I need to create a grid of points within a given area so that each points 500m across the ground from any other point. At these points trail cams will be placed, as part of a research project on biodiversity (i.e. a square with 500 meter sides and cameras placed at each corner of each square).

This terrain is in the mountains of Costa Rica. It is very steep (30 deg slopes and sudden drop-offs are not uncommon), covered in dense vegetation and has many large mudslides sporadically found throughout the area.

We can only reach these sites by foot, so I would like to be as accurate in my plotting as possible to reduce the amount of time spent wandering in the bush.

I have plotted a grid with points using MMQGIS, but it is a “flat-grid” that doesn’t seem to account for the steepness and I’m not sure what to do next.

I am fairly new to GIS and my skill set is limited (particularly when it comes to code). I am also working on a Mac, which can be problematic. I have tried searching the web (as well as this site) for detailed answers, but came up short.

I seek a detailed, step-based answer on how to do this, if it is indeed possible for me to do.

As an aside, the Internet here is unreliable and slow, it took me several days to download QGIS and it failed several times because the internet frequently disconnects. So I prefer an answer that does not involve downloading large files or software.


texts: Exporting partial image using Google Earth Engine

texts:
When exporting an image to an asset in Google Earth Engine, the task seemingly runs fine but outputs only a tiny corner of my image. It is a filtered, cloud-masked Sentinel-2 scene with 16 bands (13+ image corrections). It does not display any error messages so I assume that GEE thinks it is executing the task I have asked for. One thought I had was that I have used up too much space in my GEE assets folder but this is not the case. I have also tried a drive export and the outcome is the same. The pictures show a) what I was expecting the output to be and b) what it actually was.
Why is this happening and how can I fix it?

var geometry: polygon, 5 vertices
var collection: ImageCollection "Sentinel-2 MSI: MultiSpectral Instrument, Level-1C"

 //define the cloud mask layer
function maskS2clouds(collection) {
  var qa = collection.select('QA60');
  var cloudBitMask = 1 <<10;
  var cirrusBitMask = 1 <<11;

  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
  .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return collection.updateMask(mask).divide(10000);
}

//image collection is filtered
var dataset = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(geometry)
  .filterDate('2015-08-26', '2015-08-28')
  .map(maskS2clouds);

//image with mask applied
var image = dataset.map(function(image) { return image.clip(geometry); });
Map.addLayer(image, RGB);

//exported to asset
Export.image.toAsset({
    image:dataset.first(),
    description: '2015_CLOUDFREE',
    region: geometry,
    scale: 10,
    maxPixels: 1e13
});


texts: ArcGIS Server Search Service and mapped client drives

texts:
With ArcGIS Server (10.1), is it possible to register a folder to the Search Service for indexing using one path, and serve it to users with another? how?

For example the server sees \\geodata-server\data and \\geodata-server\maps while the clients use Z:\ and Y:\ respectively.

Under .../arcgiserver:6080/arcgis/manager/site.html >> Date Store I've registered and validated the folders, and successfully published map services with these settings, but the data store settings appear to be ignored by Search Service: search results in the clients show the UNC path.

Data Store:
   Name: Repository
   Publisher Folder Path: Z:\
   Publisher Folder Hostname: geodata-server
   Server Folder Path: \\geodata-server\data

I already tried just using the same mapped drive settings on the arcgis server, but that doesn't work (c.f. Is it possible to reference data using mapped network drives instead of UNC path in arcgis server services?)


texts: Accessing Numpy Array found within an XArray from Satpy

texts:
dask.array<shape=(76, 76), dtype=float32, chunksize=(76, 76)>
Coordinates:
    crs      object +proj=utm +datum=WGS84 +ellps=WGS84 +no_defs +units=m +zone=33.0 +type=crs
  * y        (y) float64 4.006e+06 4.006e+06 4.005e+06 ... 3.951e+06 3.95e+06
  * x        (x) float64 4.094e+05 4.104e+05 4.114e+05 ... 4.832e+05 4.842e+05
Attributes:
    sensor:               slstr
    view:                 n
    polarization:         None
    standard_name:        toa_brightness_temperature
    coordinates:          ('longitude_in', 'latitude_in')
    resolution:           1000
    wavelength:           (10.4, 10.85, 11.3)
    file_type:            esa_l1b_ntir
    long_name:            Gridded pixel brightness temperature for channel S8...
    units:                K
    level:                None
    calibration:          brightness_temperature
    name:                 S8_in
    platform_name:        Sentinel-3A
    modifiers:            ()
    start_time:           2018-07-20 09:12:25.111891
    end_time:             2018-07-20 09:15:24.803164
    area:                 Area ID: LS8_duplicate\nDescription: Patch subset i...
    ancillary_variables:  []

I have this Xarray which resulted from reading and resampling a NetCDF file using satpy.

newscn = scn.resample(area_def, resampler='nearest')

Printing the newscn results in an xarray. How can I extract the (76, 76) array with float32 values as a numpy array?


texts: Write Zonal Statistics as table to an editable table in Python

texts:
I am performing zonal statistics as table on many rasters (different years of data). I would like to write the result to a table that I can edit, so that I can rename the "SUM" column and left join the results to the previous table iteratively. This way, after my loop is done, I have columns: OBJECTID, SUM_2001, SUM_2002, SUM_2003, SUM_2004... as my resulting table.

How do I write my results to a table that lives in memory and is editable, as opposed to a file?


texts: Returning route geometry using OpenTripPlanner API

texts:
I already managed to return geometries of Isochrones using the API as follows:

https://sampleserver.com/otp/routers/myrouterid/isochrone?algorithm=accSampling&mode=WALK&walkSpeed=1.11&cutoffSec=900&precisionMeters=10&date=2017/11/22&time=14:00:00&maxWaitingTime=20&maxWalkDistance=1000&transfers=5&offRoadDistanceMeters=100&fromPlace=44.097088542141,14.5792197160854

Now I'd like to also return requested routes via OpenTripPlanner REST-API as line geometry which can be loaded into a GIS.

https://sampleserver.com/otp/routers/myrouterid/plan?fromPlace=44.170960,14.5311114&toPlace=44.127964,14.603506&time=14:00:00&date=2017/11/22&mode=CAR

returns the route as JSON, but not its geometry.

I already tried different URLs as described in the API-Docs here: http://dev.opentripplanner.org/apidoc/1.0.0/resources.html, but so far I wasn't able to return a route geometry.

Is it possible to return the geometry of a route using OpenTripPlanner API and if so, how?


texts: ArcGIS python script freezes when using ctypes MessageBox

texts:
I have written a simple ArcGIS python script tool that takes only one input, a file, and displays a "Yes/No" message box when the script is run.

The script runs, but the message box pops up behind the Geoprocessing window and does not obtain focus.  Sliding the Geoprocessing window to the side to reveal the message box, and clicking either "Yes" or "No", has no effect: the script appears frozen.  To unfreeze, i must press "Cancel" in the Geoprocessing window and then press either "Yes" or "No" in the message box. The script then cancels/aborts.

The script code is:

import arcpy
import ctypes

inFile = arcpy.GetParameterAsText(0)
arcpy.AddMessage("\nInput file: {0}\n".format(inFile))

# Yes/No mbox = 4, return code Yes = 6, return code No = 7
mBox = ctypes.windll.user32.MessageBoxW
result = mBox(None, u"Just checking...", u"File Check", 4)

if result == 6:
    arcpy.AddMessage("User pressed Yes\n")
elif result == 7:
    arcpy.AddMessage("User pressed No\n")
else:
    arcpy.AddMessage("Unknown return code\n")

arcpy.AddMessage("Complete.\n")

I realize using a message box in this context is a bit of a no-no.  But if my transgression of best practice might be overlooked for a moment, can anyone help me understand the technical reason of this message box freeze behavior?

I think it might be to do with the window handle parameter of the message box (currently set to None), and that creating it with the window handle of the Geoprocessing window might solve it - allowing the Geoprocessing window to release focus to the message box. However, i am still trying to figure out how to the get the window handle of the Geoprocessing window.

Using ArcGIS 10.5, python 2.7, Win 7


texts: Vector menu tools missing (processing plugin broken)

texts:
All of a sudden my vector menu tools are missing. I know this is a common bug but all the instructions on how to fix it didn't work for me.

I checked if the processing plugin is installed but QGIS tells me this plugin is broken.

I tried uninstalling it and deleted the .qgis2 folder but even after reinstalling nothing changed.

I'm using QGIS version 2.18.2 on macOS Sierra.


texts: Does raster interpolation recognize multiple z values for one x, y point?

texts:
I am working with lithological data, where I have multiple z depths for one set of coordinates. How is this read by the raster interpolation toolset? Does it utilize only one or all of these z values?

If it does not recognize multiple values, is there a type of interpolation I can use that can?

My ultimate goal is to create something in ArcScene that has a 3D appearance and can demonstrate the profundity of a certain lithology.


texts: Postgis: Query ST_DWithin in square instead of circle

texts: The ST_DWithin query queries all geometries within a certain radius from another geometry (in my case a point). For performance optimization of my application I want to query in a square. So that the query gives me back all geometries within a certain square of my point. Is this possible? Because I don't find any kind of query to this. My geometries are polygons/multipolygons.


texts: Convert geographic coordinates to UTM zone 29N with Javascript

texts:
I have an issue converting lat longs to UTM with the ArcGIS javscript API.

Lets say that I have a map, based on the UTM zone 29N, also known as WKID 32629.

I have a known set of coordinates from a GPS, in this format:

Latitude:  62.220596
Longitude: -6.564331

The geograpic coordinates are in WKID 4326.

To convert using javascript I have been trying to use the Geometry Service's project method.

If we look at this post in the ArcGIS forum, it is about the same issue. Here, I have tried using the same code, and adjusting the wkid of the output spatialReference to 32629, and you see it returns only NaN, try in the JSfiddle to change the wkid to 4326, and you will see geographic coordinates.

Edit 1

Looking at the HTTP requests made to the server, I see this url is requested (now using the map of the Faroe Islands)

http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer/project?f=json&outSR=32629&inSR=4326&geometries={"geometryType":"esriGeometryPoint","geometries":[{"x":62.01387,"y":-6.7866,"spatialReference":{"wkid":4326}}]}&callback=dojo.io.script.jsonp_dojoIoScript2._jsonpCallback

This looks like a correct url, containing outSR, inSR, geometries, and callback. However, simulating requests to this url gives me either a 500 Internal Server Error, or this string:

dojo.io.script.jsonp_dojoIoScript2._jsonpCallback({"error":{"code":400,"message":"Unable to complete Project operation.","details":["Error projecting geometries"]}});

Any ideas?


texts: Calculate land cover for grid cells in QGIS

texts:
I have a raster containing information on the land cover type in a region.
I also have an ESRI shapefile consisting of a grid composed of 2500 cells. I would like to calculate the percentage cover of each land cover type in each cell.

To do this, I vectorized my raster and dissolved the polygons to obtain a shapefile where each polygon is a land cover type. I then tried to use the intersection tool to obtain a shapefile with polygons for each land cover type in each cell. This procedure takes too much time however (days by my calculation). I have also tried to split the grid into smaller grids but still faced the same problem.

Is there a more effective way of obtaining what I want?


texts:
I have 300 polygons in a single shapefile, I want merge some polygons together, when I merge two polygons the attribute does not come total count of two polygons,

What I need is, when I merge two polygon the attribute value has to add together and the total value has to come in merged polygon attribute.


texts: PostGIS 2.0 st_geomfromgml using LineString fails with invalid GML representation

texts:
I am attempting to insert a geometry into a table using the st_geomfromgml function.

The GML is being sourced from a Geoserver WFS request.

The function call looks like this

st_geomfromgml('<?xml version="1.0" encoding="UTF-8"?><gml:LineString xmlns:gml="http://www.opengis.net/gml" srsName="http://www.opengis.net/gml/srs/epsg.xml#4326" srsDimension="2">
 <gml:posList>146.48296 -18.6222 146.48296 -18.6222</gml:posList>
</gml:LineString>');

This worked on PostGIS version POSTGIS="1.5.3" GEOS="3.3.1-CAPI-1.7.1" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.7.3" USE_STATS

but fails on POSTGIS="2.0.0 r9605" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.8.0, 6 March 2012" LIBXML="2.7.8" LIBJSON="UNKNOWN"

However I can make it work if I specify three dimensions like so;

select st_geomfromgml('<?xml version="1.0" encoding="UTF-8"?><gml:LineString xmlns:gml="http://www.opengis.net/gml" srsName="http://www.opengis.net/gml/srs/epsg.xml#4326" srsDimension="3">
 <gml:posList>146.48296 -18.6222 0 146.48296 -18.6222 0</gml:posList>
</gml:LineString>');

Is there a way to keep the geometry in two dimensions and insert it?


texts: How to convert Distance, Azimuth, Dip to XYZ?

texts:
I have an excel spreadsheet with header and survey drill data. Header data contains hole ID, and location coordinates, and the survey data contains related downhole survey with Distance, Azimuth and Dip values.

Since I know the hole location and surface elevation, I would like to be able to convert the survey table to XYZ coordinates as well.  Does anyone has a function / procedure / example?  (VB & ArcObjects)

Header Data:

Survey Data:


texts: random sampling from spatial points in R

texts:
Can anyone help me out on the best package to sample from a spatial point dataset?

My dataset consists of a list of positions inside a region, and for each point I also have the name/number of the subregion where the point is located. I would like to make my sampling using only this dataset (for example I don't want to load a shape file with the subregions or a probability raster).

I have already made a "random" sampling for each subregion, however I am not taking into account the spatial distribution of the points inside each subregion.

I am looking for a package that would let me do both spatial random sampling as well as stratified random sampling, or some guidelines on how I can create my own code.


texts: Extract raster values to multipoint shapefile and write to new field using gdal/python

texts:
I have a raster with one band and I would like to query a set of overlying points (.shp) to determine how many points as well as which points points fall under a range of certain pixel values associated with the raster.

This is straightforward in ArcGIS but I want to use GDAL/OGR/Python to automate this process.

I can think of a few ways to potentially accomplish this:
(1) Join the raster values to the points, add a new field with the raster values and then query the resulting new_points.shp file.
(2) Write an effective SQL statement that queries the points and the raster to write out a new file with only those points that correspond to the value range I'm interested in. Potentially using a VRT file
(3) Convert the raster to a shapefile using gdal_polygonize (very time intensive!) and then use an SQL statement along the lines of SELECT COUNT(*) FROM points, new_poly WHERE ST_INTERSECTS(new_poly.geometry,points.geometry) -dialect SQLITE input.vrt again using a vrt table.

Number (2) seems like the best and most efficient option, but I am unfamiliar with querying both raster and vector data simultaneously. Is there a standard or accepted way to do this?

EDIT: Another option seems to be the gdallocationinfo utility, but it seems that you can only pass lat/long coordinates and not shapefiles.


texts: Converting aeronautical information coordinates to Google Maps API coordinates

texts:
I'm in the process of creating a webpage to ease the process of visualizing NOTAM navigation warnings. Never done any geospatial software development before, so all the terms like geoid, reference ellipsoid, datums, different coordinate systems are new to me, though I've been doing a lot of reading up on the related material. I'm using the DotSpatial.Projections library to load the projection info either based on the EPSG code or the corresponding proj4 string. The target coordinate system is the one used by Google Maps API - (WGS84 / EPSG 3857).

I'm having trouble figuring out the proper conversion source coordinate system/projection. What I've come up with does not give the proper results: https://epsg.io/transform#s_srs=42310&t_srs=4326

The source coordinate system is described here: http://eaip.eans.ee/2018-02-01/html/eAIP/EE-GEN-2.1-en-GB.html#GEN-2.1-Geodetic

The converted nor the original coordinates don't align with the expected result. For example, the Zone 9 near north of Estonia surrounding Tallinn: Expected result, actual result on my website (large red/green polygon south of Tallinn)

Some proj4 strings I've come up with:
NOTAM:

+proj=lcc +lon_0=24.0 +lat_0=57.51755393 +lat_1=58.0 +lat_2=59.3333333 +ellps=GRS80 +datum=WGS84 +units=degrees +no_defs

Google Maps:

+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

The aforementioned two strings result in a very tiny projection near the reference point of the source coordinate system.

This application is initially for personal use as a drone pilot to visualize if there are any flight restriction zones active at a given time.


texts: External editing of .tiff, retaining georeferencing?

texts:
Is there any way I can reduce the colours of a TIFF to make conversion to polygon more efficient. There's quite a bit of artification on my original and I'd like to reduce it to 2 colours before the conversion (which I have done, however once re-saved it loses the referencing).

Reduced coluor image:

Original TIFF in QGIS

Output polygon with solid white removed


texts: Activate the editing tool "Move feature" for Oracle layers in QGIS?

texts:
In QGIS I have a map layer which is connected to an Oracle database.

After I turn the map layer to editing mode I discover that I can only change the attribute data, create and delete objects, but I cannot move objects on the map since the move feature tool is deactivated.

Does anyone know what is the reason and is there any solution to this problem?

As shown in the picture, a few functions are deactivated in the editing toolbox even though the layer is in editing mode:


texts: Building a private GPS almanac server - practical experience

texts:
I'm trying to provide a source of GPS almanac files (in SEM format) for a non-connected / sometimes connected network.

I know that I can download them from the USCG NAVCEN. However I don't want to have to air-gap them across to my network.

So the plan is to (somehow) get them from a GPS receiver, format them into SEM, and serve them up using a file share / web server. I'm confident about the last step. I have a rough idea of how to download the almanac, but no code. I've seen some online resources that allow writing SEM, but nothing that provides all of

a simple API
open source
(obviously) maintained

So at this stage I'm considering writing the "missing parts". Has anyone actually tried this? Any recommendations for usable (open source) blocks.


texts: Update multiple different layer data sources in mxd with multiple lists

texts:
I need to update 3 layers in a template mxd then save that mxd. It needs to loop through a fGDB and for each feature, update the three layers then save a copy for each loop.

What I have loops through and updates the layers but I can't seem to figure out how to save each mxd consecutively after each loop. It currently only saves a copy for the last feature in the loop.

# Local variables:
newMXD_wks = ("R:\\Projects\\Active\\Inundation_Mapping\\Inundation_Maps\\")

Current_Grid_wks = "R:\\Projects\\Active\\Inundation_Mapping\\MapBook_Grids.gdb\\Template_Grid"
New_Grid_wks = "R:\\Projects\\Active\\Inundation_Mapping\\MapBook_Grids.gdb"

Current_Structures_wks = "R:\\Projects\\Active\\Inundation_Mapping\\Structures_fGDBs\\Structures_Final.gdb\\Template_Structure"
New_Structures_wks = "R:\\Projects\\Active\\Inundation_Mapping\\Structures_fGDBs\\Structures_Final.gdb"

Current_Inundation_wks = "R:\\Projects\\Active\\Inundation_Mapping\\Inundation_Polygons_Final\\Oconee.gdb\\Template_Inundation"
New_Inundation_wks = "R:\\Projects\\Active\\Inundation_Mapping\\Inundation_Polygons_Final\\Oconee.gdb"

arcpy.env.workspace = New_Grid_wks
grids = arcpy.ListFeatureClasses('Oconee*')

for grid in grids:
    mxd = arcpy.mapping.MapDocument("R:\\Projects\\Active\\Inundation_Mapping\\InundationTemplate.mxd")
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.supports("DATASOURCE"):
            if lyr.dataSource == Current_Grid_wks:
                lyr.replaceDataSource(New_Grid_wks, "FILEGDB_WORKSPACE", grid)

            if lyr.dataSource == Current_Structures_wks:
                arcpy.env.workspace = New_Structures_wks
                structures = arcpy.ListFeatureClasses('Oconee*')
                for structure in structures:
                    lyr.replaceDataSource(New_Structures_wks, "FILEGDB_WORKSPACE", structure)

            if lyr.dataSource == Current_Inundation_wks:
                arcpy.env.workspace = New_Inundation_wks
                inundations = arcpy.ListFeatureClasses('*Final')
                for inundation in inundations:
                    lyr.replaceDataSource(New_Inundation_wks, "FILEGDB_WORKSPACE", inundation)

mxd.saveACopy(newMXD_wks + inundation + ".mxd")
del mxd


texts: Calculate the sum of the areas of polygons within bigger polygons from another layer _ QGIS

texts:
I would like to add in the attributes of each of my roof boundaries (red polygons, which sometimes overlap) the sum of the surfaces of the (green) surfaces.
I have already tried "join attributes by locations" but some of the green surfaces are sometimes just partially on the red limits.
It leads to an error:

the area added in the attributes is the whole area of the green
surface and NOT only the partial overlap

What would be the easiest approach for working with Qgis, taking into account that there are a lot of data?


texts: I tiled hillshading RGBA GeoTIFF with gdal2tiles but produced tiles has no transparency and instead it is blended to black. Is there a way to preserve transparency?


texts: Schema locking logistics on ArcSDE for SQL Server?

texts:
I have a set of feature classes that I want to be able to edit schema or what ever I want to do in the heat of the battle and I don't want to be stiffled by schema locks.

I have a interpretable requirement to have all analysis done in the db.  I have them outside of any feature dataset but  I guess that there may be schema locks for any client hooked to the geodb.

Is the answer to house this in a feature dataset?  Or are there other alternatives?


texts: Buffer Different Distances Based on Field Values

texts:
A while back I was able to create buffers around point data of different distances in miles based on a range of field values.  Somehow, I've forgotten how I accomplished it and I'm hoping someone can assist.  I have about 30 points and each point has a range of values from 0-50. What I'd like to do is:

0-20= 50 mile buffer
21-35= 75 mile buffer
36-50= 125 mile buffer

I could always separate them, buffer each one and then merge them, but I'd prefer a quicker solution and possibly one that I could sketch out in ModelBuilder.


texts: What is the range of numbers for coordinates in geographical coordinate systems

texts:
Up to now I thought the range of values in the geographic coordinate systems was:

Latitude (S-N):   -90 to  +90
Longitude (W-E): -180 to +180

I have set the coordinate system of my ArcGIS data frame in an empty map to WGS 1984 and I'm getting numbers like the following when I move the mouse over the window:

What does it mean?


texts: Can't view PostGIS rasters after changing schemas

texts:
I have loaded a few rasters into the Public schema of a PostGIS 2.0 database. Using the db manager in QGIS 1.8.0 I can add the rasters to the map canvas. After moving the rasters to another schema in the same PostGIS database, I am unable to view the rasters in QGIS. However, moving the raster back to the public schema I am again able to view the rasters in QGIS.

Why is this the case? and
is it possible to fix?

I have tried to load the raster directly into another schema using the someschema.sometable syntax in raster2psql but this returns an error stating that my target schema does not exist, however it does exist in the database. Therefore, I have left out the someschema.sometable syntax resulting in the rasters being stored in the public schema and manually moving them to another schema.

EDIT:

Some additional information.
When loading a raster to PostGIS a rid and -i column are created. Only after changing the -i column to rast does QGIS recognise the raster layer (layer is still in public schema). However, the QGIS error "not a valid layer" returns after moving the raster table to another schema.

EDIT2:

I have figured out why the raster column heading was named -i. I used the -f parameter without specifying the colummn name. The -I constrained directly followed the -f constrianed therefore resulting in the raster column name being -i.


texts: Adding pixel value data from a raster layer to a number of buffer polygons ArcGIS

texts: I have created a layer of polygons with buffers around my points of interest, all of the same size, in ArcGIS. I also have a raster layer of nighttime lights with pixel values between 0 and 63 clipped to the country of interest. I would like to know how to add the average pixel value of the raster within each buffer area to the attribute table of the buffers. The attribute table of the night lights has 3 columns, the object ID, value and count.


texts:
We were using Bing maps to place pins on a map when a user loads a page.

The number of pins has grown significantly so we started to cluster the pins to reduce the load on the client. Now even clustering is placing too much load on the client.

I'm asking a general question of what to do next.  It would seem that a server-side solution is required as expecting the client to do all of the processing with an ever-increasing (over 50000) number of pins seems to be a road to nowhere.

The website is written in ASP.Net Mvc. Can anyone point me in the right direction?