texts:
I am working on a project in which I require road network.

I have downloaded data through OpenStreetMap but it is not updated.

Is it possible to download Google Maps data with attributes?


texts: How to convert a higher resolution shapefile to the resolution of another shapefile?

texts:
I have a shapefile with ~110,000 grid cells which I want to convert into the same grid cell sizes as in another shapefile with ~27,000 grid cells? Is there a way to do this in QGIS?
Details about the two shapefiles:

Higher resolution shapefile: ~110,000 rows in attribute table:

Lower resolution shapefile: ~27,000 rows in attribute table:

But both have the same country boundary. - India:

I basically want both to have the same number of rows which in a way means same resolution.


texts: Removing disconnected cells in a stream network raster using ArcGIS

texts:
I've got a predicted stream network raster with cells values 0 (not a stream) and 1 (is a stream). My end goal is to create a network of connected polylines, but unfortunately there are some areas that are disconnected that we want to include in the final file, and some we do not - we have decided to get rid of any segments disconnected from the larger network by more than 2 cells.

Therefore, I would only like to keep the cells with value 1 that are either 1. connected to the larger network or 2. within 2 cells (~40ft) of this main network. The areas circled in green below are what we're hoping to maintain, while removing the areas circled in red.

I have tried using Spatial Analyst --> Region Group then reclassification (using a cutoff of count = 7) to remove the "orhpans" as we call we them, but that also gets rid of some areas that shouldn’t be cut, like a segment that’s short but only disconnected by 1 or 2 cells. See image below - areas in red should be maintained but are not.

I understand why using the cutoff of 7 wouldn't necessarily work, as it would cut out small segments that are disconnected by 1 or 2 cells. I'm not sure what approach I SHOULD use, however.


texts:
I am trying to iterate through a set of 50 points of GPS data. Each have their own unique Estimated Horizontal Accuracy. What I am looking to do in ModelBuilder is create a model so that it creates a buffer for each unique value (using Iterate Row Selection or some other iterator), and then combine them all based on collection site at the end.

I have tried iterating based on multiple values (site name and horizontal accuracy, with the horizontal accuracy at the top of the list), but it did not read properly. I also did just the horizontal accuracy, but that also did not work. Every time I run it, it gets to the buffer part, and then ArcMap crashes with an error report. I think this is because of a runtime error, and, therefore, I have a parameter that is wrong in my buffer tool. The horizontal accuracy value is a Double value in the original table/data layer.

I have "Get Values" set to read the Sitename attribute.

Any ideas?


texts:
I've been given Ground control points for job and the vertical co-ordinates are in AHD (Australian Height Datum). Does anyone out there have any idea how to transform vertical heights in AHD to Pix4d accepted vertical coordinate system?

e.g. MSL emg88 or 94 or08.

Also I have access to ArcGIS, QGIS etc for transformation etc..


texts: Polygon to Polyline tabular data transfer

texts:
I am working with a Divisional database. The division is divided by districts and then the district is further sub-divided by its micro units called sub-districts. Its topological relationship was error free and in the table the polygon data were well defined in its sub-district, district and divisional columns. Now I have drawn its road network throughout the division. There are some highway with unique name that covers more than 6-7 sub-district polygon.

Now I want to create a road database that will break the road at the sub-district boundary line and join the polygon data (Sub-district, district and divisional data) to the clipped road with its existing database. I think it is possible but I couldn't find any solution. Can anybody help me?


texts: Raster image changes colour after import

texts:
after raster .tiff image import in project, it changes colour.

Imported

Original

What could be the options in QGIS 2.18.16 to get back to original colours?
This kind of problem is quite new, appeard in this week:)


texts: NDVI analysis of LANDSAT 8 imagery in QGIS

texts:
I am trying to run an NDVI classification of 2015 multispectral LT8 imagery in QGIS to assess the quality of vegetation cover in regards to habitat suitability for a primate species in northern Vietnam.

Can anyone tell me if I need to make a reflectance conversion of the bands that I will use (band 4 and 5) or can I just assign the bands and use the (NIR-RED)/(NIR+RED) equation through the raster calculator?

How do reflectance conversions such as DOS1 affect imagery analysis?


texts: What is the right solution (and map projection) to compute distances between points located all over the World?

texts:
After my study on map projection I realized that it is very difficult to choose a projection that allows to compute simultaneously the distances between a great deal of points if these points are located all over the world. Equidistant projections allows only to compute distance between couple of points.
Then I decided to use the following solution. My data set can be split in smaller datasets:

United States territory;
the whole EURASIA;
the territory encompassing the African continent and the Europe.

In ArcGIS 9.3 (ArcInfo license) I projected each dataset with the Lambert Azimuthal Equal Area projection properly “centered” for each case. For each dataset I loaded a new layer with the Tissot's Indicatrix shapefile and computed the areas of each ellipse. I considered as "not distorted" (as regards the area, and the distances) all that territories over which the ellipses’ areas do not differ more than the 5% from the ellipse nearest to the center of the projection and when the ellipses' shape do not reasonably differ from a circular one. After having  seen that my points were enclosed in this “not distorted” territory, I computed the point distance.
Do you think my procedure was reasonably correct or have some advice about a better map projection for my purpose?


texts:
I'm trying to figure out an automated way to identify mountainslopes from where a rockfall could hit lakes using the hydrological tools of ArcGIS 10.7. As rockfalls do not behave like runoff water, some improvisation is needed to get only those slopes from where rocks would hit the lake directly.

(I've asked a previous question about this, see Create line feature from aspect value in ArcGIS - but while the answer is correct, the computation time is enormous and not suitable for larger applications.)

Let me try to clarify my question with the picture below. The slopes with the blue runoff arrows drain into the lake, but rockfalls from there would not necessarily reach the lake. From the slopes with the yellow runoff arrows, however, they would hit the lake.

So my new idea would be to identify the main inlet of the lake, make it the pour point and exclude all slopes that drain into this pixel (blue arrows) while keeping all slopes that drain into the outlet of the lake (yellow arrows). This way, I would get more or less only the slopes directly adjacent to the lake.

My question now is: How do I identify the "main inlet" of a lake?

I've tried masking the flow accumulation raster with the lake outline line feature and then take the pixel with the second highest value. Unfortunately, the runoff inside the lake area sometimes overlaps with the lake border, so I don't get the inlet but rather a point further downstream.


texts: I am trying to create an index page using a customized index polygon featureclass that I created. Is there a way to do that without using the grid index features tool, which is more for set sized grids, but for a customized polygon set?


texts:
I'm looking for a rasterized world map using the EPSG4326 projection, but haven't been able to find any ones online (or at least decent quality ones that appear that the top left does correspond to -180, 90). Perhaps I don't know what to search for.

Can anyone point me to a resource with such a pre-rasterized image?

If not, how would I go about producing such an image?


texts: Error: Coordinate values are out of range [-180 -90, 180 90] for GEOGRAHY type on Postgis?

texts:
when I run a query on Postgis

select st_astext(ST_SetSRID(f.geom,4326) ),st_astext(st_centroid(ST_Transform(b.wkb_geometry,4326))) from a,b where ST_Distance_Sphere( ST_Transform(a.geom,4269),ST_Transform(b.wkb_geometry,4269))<10 limit 1;

I got this message on my screen print in infinite loop

 Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY

to give you an idea about what I am getting when I fetch a.geom and wkb_geometry ?

select st_astext(ST_SetSRID(f.geom,4326) ),st_astext(st_centroid(ST_Transform(b.wkb_geometry,4326))) from a,b limit 1;

I got this :

POINT(-73.9798264975362 40.679599746172) | MULTIPOLYGON(((989878.855200246 186914.37490882,989831.404507741 186822.841299236,989812.018719748 186832.772053644,989859.469740331 186924.306319401,989878.855200246 186914.37490882))) |

what does that mean ? Does this affect results I am getting ?


texts: Constructing spatial weight matrix from triangulation of point coordinates in R?

texts:
I want to know how to construct spatial weight matrix from point coordinates in R,such as using the Lucas County, OH, housing point data set in spdep package.

I read here,it says'the list of neighbours provided with the data set in spdep is a sphere of influence
graph constructed from a triangulation of the point coordinates of the houses
after projection to the Ohio North NAD83 (HARN) Lambert Conformal Conical
specification (EPSG:2834).'

I only know how to construct spatial weight matrix from polygons,no idea about points.

There is my code:

library(spdep)
library(rgdal)
data(house)
writeOGR(obj=house, dsn="test", layer="house", driver="ESRI Shapefile")
temp=readOGR("./test/house.shp")


texts:
I'd like to set up a webmap that clients can drop into their existing website.
I'm thinking Geoserver, OpenLayers & GeoExt.
It would be very simple slippy map kind of thing. Do I need Geoserver at all then?
Just want to be able to use WMS/pan/zoom/turn layers on and off.
If I ditch GeoServer, then do I need Apache?

And where would be the best place to host this?
Could most hosting companies handle it if there was no Geoserver since OL & GeoExt are just JavaScript? I don't want to host it myself.

Anyone who could point me in the right direction for a checklist or roadmap to offer up simple webmaps is greatly appreciated. I'm looking for info on setting up the web side of it, since I'm fairly good with GIS already.


texts: Python 2.7 stand alone script for compounded with statement. Errors 'column specified does not exist - but it does

texts:
I'm writing a standalone Python script to use for labeling ArcMap 10.7 data. (Want to avoid label expressions for portability, etc). I'm attempting to use ID attribute to label based on ADMIN_ORG attribute data.
I have found from SO, that nested 'with' statements are better for this than nested 'for' statements. I tried to replicated the answers from this post ArcPy nested loop problem  still having issues....
The script currently DOES print the data from the ID column as expected but errors at the end at
"for label in sc_label:"
Runtime Error: A column was specified that does not exist. But the column absolutely does and was recognized in earlier attempts. What am I missing?
def newlabel():
    fields = ['ID', 'ADMIN_ORG']
    label = fields[0]
    region = fields[1]
    with arcpy.da.SearchCursor(fc, region) as sc_region:
        for region in sc_region:
            with arcpy.da.SearchCursor(fc, label) as sc_label:
                for label in sc_label:#errors here but still prints the attribute data
                    print str(label[0])


texts: When to add impedances to network dataset?

texts:
I am using ArcGIS 10.2 for desktop and will be using point and click to digitise a cycle network shapefile and network data sets then use model builder or python.
My aim is to assess the accessibility of my city by bicycle. To show impedances such as distance and slope, identify barriers such as intersections, and gaps in connectivity.

I am hoping to create a multi-modal network data set including an off road cycle path network, on road cycle path network and complete road network. I was wondering:

A - if it's best to keep these separate particularly the off road and on road cycle paths or merge them before building the network dataset?

B - whether to add surface information such as slope to these layers before building the network dataset?

C - how to produce a map with a bicycle level of service for each segment of the network, giving a classification e.g network segment that is flat, has an off-road cycle path would be identified as "Excellent"?

Would this best be done in ModelBuilder or python using the select tool?

Aim to visualise a bicycle suitability map and eventually an accessibility and bicycle friendliness.


texts:
I'm trying to replicate something similar to the ArcGIS WPF Runtime SDK (C:\ArcGIS\WPF1.0\SDK\Samples\Sample Application\Sample Application.exe) "Routing with Barriers Online" sample with two major modifications ...
1. I have my own network dataset (created from a feature class using ArcObjects classes)
2. I want it to run locally and not on an ArcGIS server instance. Can I exploit LocalServer for this?

First, can this be done? If so, what is the best approach? I have both the WPF runtime and ArcObjects SDK to work with.


texts: Creating non-spatial Many to One Joins/Relationships between tables

texts:
I may have missed a trick here but I can't seem to find a solution.

I have a table of ~1.6 million postcodes, also containing co-ordinates, and a separate file detailing ~6000 a population with postcodes attached. I am trying to attach the co-ordinates from the postcode file to the population file. There is, but not always, two or more individuals associated to one postcode.
The theory, in my head, is simple. Join the postcode table to the population table on the postcode attribute. To start with this idea in ArcGIS 10.1 seemed simple, using a basic join, but what seems to be happening is that a 1:1 relationship is being enforced. After researching further I found the Make Query Table tool which claimed to solve this M:1 problem. Though that doesn't seem to be working.
I am currently in a two-step process to see if I can solve this. First is looking through some, don't particularly want to look through all files in either table, to see if there is a formatting error that would explain this error. The second is looking into using MySQL to see if that can be used to join the tables. Or at least create a query of the two tables that will produce the desired output table.

First off if you have any advice on making this work in ArcGIS that'd be brilliant, I'm sure it can do it, just not sure how to make it work. Alternatively if you know of another method that would easily join these tables, be that QGIS, MySQL or something else, then please do say.

Hope you can help, and if there's anything I've missed out that would help you understand my problem, let me know.


texts:
I have used the qgis2web plugin successfully before but now my exported maps all have a fixed canvas size that is larger than the display screen.

Here is the Head section of my HTML file as I suspect it is here where the problem is:

<head>
    <meta charset="utf-8">
    <meta http-equiv="X-UA-Compatible" content="IE=edge">
    <meta name="viewport" content="initial-scale=1,user-scalable=no,maximum-scale=1,width=device-width">
    <meta name="mobile-web-app-capable" content="yes">
    <meta name="apple-mobile-web-app-capable" content="yes">
    <link rel="stylesheet" href="css/leaflet.css">
    <link rel="stylesheet" href="css/qgis2web.css"><link rel="stylesheet" href="css/fontawesome-all.min.css">
    <style>
    #map {
        width: 1617px;
        height: 834px;
    }
    </style>
    <title></title>
</head>

What is setting the  "width: 1617px; height: 834px;" statement?

Thanks.


texts: GIS "hack": ways arbitarily assigning coordinate in order to match the raster and vector layer in QGIS

texts:
Before I raise my question, I have to admit that I am a total newer on GIS, I don't have much background on GIS. But suddenly I need to take over the GIS job from a GIS specialist who just left the office. And since time is a luxury for me right now, I am looking for some "hack" to fix my immediate problem.

The task is below:
I need to spot out which bed in an elderly home belongs to a male or female resident, and present it visually on a map.

The materials are below:
- A raw map from the elderly home, marked the bed no for each bed
- A spreadsheet contain the bed no and the gender of the resident on that bed

What I have done till now:
- "digitize" the raw map: I have used MapWindows and created a vector layer consist of beds, something like creating a new layer on photoshop, and then trace the bed on the new layer
- created a new project in QGIS, import the raw map as raster layer, the bed layer (vector layer created by MapWindows), and the spreadsheet as csv
- in QGIS, I have joint the bed layer with the csv file correctly
- made adjustment in the "Style" tab of the bed layer (vector) so that the gender of each bed can be shown

I think I am almost there, but problem arise:
- in Mapwindows, both the raw map and the vector layer matched perfectly (i.e. I draw the trace the bed on the new layer and both the raster and vector layer align to each other perfectly)
- this is not the case in QGIS

I am thinking of something like "georeferencing" issue, which the GIS specialist keep saying before, but I don't know what the hack it is.

I need to "realign" the raster and vector layer, how can I do it?


texts: Use GeoTIFF masks to optimise combine layers algorithm

texts:
My app is a C# program that takes in an extent from the user, then loads GeoTIFFs from a folder that overlap the extent. It layers them according to a simple number on each GeoTIFF, then exports the result based on the user-supplied extent and saves that cropped version to a new file. I use the Global Mapper SDK for this and it works great.

I want to speed up the app by not loading any GeoTIFF which is already "covered" by other data. I notice that the GeoTIFF format has a mask. For example, I can see in the thumbnails for the GeoTIFF files in Windows Explorer a black background with a white foreground that is the actual data.

Is there a way to read this mask into a .NET image or similar? Then I would be maintain a growing bitmap of all masks ORed together, and able to skip loading any GeoTIFF that did not increase the white area. That would greatly reduce the time taken for the GM SDK to do the export because there would be less total layers loaded.


texts: QGIS 2.18.10 Sentinel2 RGB clipped raster issue

texts:
I have been using QGIS with Sentinel2A imagery, and i just noticed something wierd to the extent of my knowledge which should have an explanation.

Consider the 3 bands of Sentinel2, 4,3,2 that are the Red, Green and Blue ones and the following steps:

Lets make an RGB composite raster with QGIS GDAL Merge tool with the 3 full bands from the Sentinel imagery.
Now crop the RGB output image to the extent of the desired region of interest (ROI) with GDAL Clip Raster by mask layer, inserting the ROI polygon and saving the output: "clippedFromWholeNDVI.tif" (default options are left intact).
Now lets crop the 3 Sentinel bands to the extent of the ROI and then construct in the same way the RGB composite.

The results are different in regard to their values and visually.

Can you explain why this happens?


texts: How to slice a long track into 100km segments?

texts:
I'd like to slice a single-segment, 1,000km-track into 100km segments.

I tried the following applications, to no avail:

GpsTrackEditor 1.15 beta build 141: N.A.
Gpx Editor 1.6.18.1649: N.A.
GpsPrune 18.6 Build 343: Track > Split track into segments > Distance limits = 100km : "The track could not be split"

Do you know of any (Windows) application that can do this?


texts: What projection to use in POSTGIS for LAT/LONG and 'contains' queries

texts:
I would like to store latitude / longitude information in a POSTGIS database. I will need to query that data with simple contains query (for display on a map, so the bounding will be simple box) and lookup by lat/long.

As I'm new to all this GIS, what would be the recommend storage/column format? Anything that needs to be considered for indexes?


texts: Is there any methodology or technology to post process the Mobile/Tablet GPS data

texts:
We have developed an application which is running on android OS for collecting the Lat/Long value of any feature and it runs in any type of android Tablet/Mobile but the problem is that the accuracy is very poor (because the accuracy is depends upon the GPS receiver of that device).

Is there any methodology  or technology to post process the GPS data?


texts: Extract Values to Point "RASTERVALU" rename field

texts: Is it possible to rename the RASTERVALU field on the attributes table, perhaps somewhere in the code?  I need to loop thru this command for several raster files due to "Extract MultiValue to Point" tool being buggy on version 10.4.  I know it is possible to do this with the multivalue tool in the arcpy code.  Or is there another solution?  Maybe add a blank field ahead of time and copy raster data to that field.  I'm willing to try anything at this point.


texts: Counting how many times a point goes through a loop and write the count to a text

texts:
I have a several feature classes that go through a loop to calculate points.  I want to add up how many go through it each time so I can archive that number.  I'm using ArcGIS 10.1.   Lets say I have 100 shapefiles with anywhere between 1 - 50 points per shapefile. I want to run a count on the total number of points in those 100 shapefiles. Right now it's printing out the total number of points in each shapefile loop. Is there a way to add the 'count' up to create a total?

SMON = arcpy.ListFeatureClasses("SectMon_PLSS_*")
for fc in SMON:
        filename =  "'" + arcpy.env.workspace + os.sep + fc + "'"
        arcpy.env.overwriteOutput = True
        arcpy.AddField_management(fc, 'shpname','text')
        arcpy.AddField_management(fc, 'GIS_ID_PAM','text')
        arcpy.AddField_management(fc, 'LATDMS','text')
        arcpy.AddField_management(fc, 'LONGDMS','text')
        arcpy.AddField_management(fc, 'DESC_','text')
        arcpy.AddField_management(fc, 'StatePlane','text')
        arcpy.AddField_management(fc, 'State','text')
        for fc in SMON:
            result = arcpy.GetCount_management(fc)
            count = int(result.getOutput(0))
            x.append(count)
            c = sum(x for x in x)
            print "The cumulative sum for this iteration is %s" %c

This is what the printouts are giving me:

 P:\Scripts\SystemData\TestData\TO PAM\140306_BCR_PNTS102034-102035&103003-103005 (1).shp
 This is located in Colorado North
 This is located in Colorado North
 This is located in Colorado North
 This is located in Colorado North
 This is located in Colorado North
 Fields Calculated
 The cumulative sum for this iteration is 5
 #############COLORADO NORTH FINISHED################
 P:\Scripts\SystemData\TestData\5pm\140226_BCR_PNTS100463-100465.shp
 This is located in North Dakota North
 This is located in North Dakota North
 This is located in North Dakota North
 Fields Calculated
 The cumulative sum for this iteration is 3
 P:\Scripts\SystemData\TestData\5pm\140311_BCR_PNT100466.shp
 This is located in North Dakota North
 Fields Calculated
 The cumulative sum for this iteration is 1
 P:\Scripts\SystemData\TestData\5pm\140312_BCR_PNTS100467-100469.shp
 This is located in North Dakota North
 This is located in North Dakota North
 This is located in North Dakota North
 Fields Calculated
 The cumulative sum for this iteration is 3
 #############NORTH DAKOTA FINISHED################

I want to be able write to a log file each time it goes through a loop.  I'm trying this, but it's giving me an error.

            txtFile = open(count,"a")
            txtFile.write("This is the total amount " + str(c))

It's telling me "coercing to Unicode: need string or buffer, int found"


texts: Devising measure of vertical distance from spring line?

texts:
I am building a model (using Logistic Regression) to explore the relation between a binary categorical dependent variable (parcels agricultural production quality) and several predictors, among which slope, curvature, distance from the coastline, etc.

One of the predictors I wish to include in the model (the significance of its contribution pending) is distance from spring line. I have digitized the spring line as a polyline.

The spring line is usually rather high in elevation (in the territory under study) and it lies where two geological layers meet up.

Now, I would like to devise a distance measure which allows to discriminate between areas below the spring line and above it.

The rationale of it would be that, due to gravity, areas (and related parcels) above the spring line could possibly be more difficult to irrigate, whereas areas (and related parcels) below the spring line would be easier to irrigate.

I hope my poor English manages to make my point.


texts: Raster to point - not in the centre of the cell

texts:
The raster to point does the following:

The raster will be converted to a point feature class using the ArcGIS Raster to Point tool. This tool creates a point at the center of each raster cell, except noData cells.

My questions is if there is a way to convert a raster to point but not creating a point at the centre of the raster rather than in one corner? This is to do an experiment of how different distributions across the cell would influence my analysis.

ArcMAP 10.2


texts: How to mosaic 3 different layers that do not overlap and make them into 1 layer

texts:
I have 3 different independent layers which I want to merge into one layer:

For creating those layers, I have calculted NDVI, calculted the mean and then I have created one layer that is lower than the mean value, one layer that is higher and one ehich is the mean. I want now to take those 3 layers and merge them into one.

I have read all the section of the mosaicking but as I understand, I can't just merge those 3 layers.

Is there any way to merge those layers into one ? I thought to create visualization based ont their values but I saw I can't use this calssify.

This is how I created the 3 different layers:

var MyImage=ee.Image(withNDVI.first());

Map.addLayer(MyImage,colorizedVis,'ndvi');

//var std2 = ee.Number(STDDictionary.get("NDVI")).divide(2);
// var mean1 = ee.Number(meanDictionary.get("NDVI"));

///////the good function I want to work with:///////
for (var i in MyImage){

  var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});

var meanDictionary = MyImage.reduceRegion({
  reducer: reducers,
  bestEffort: true,
});

}

print(meanDictionary);

//calculation
var std2 = ee.Number(meanDictionary.get('NDVI_stdDev')).divide(2);
var mean1 = ee.Number(meanDictionary.get("NDVI_mean"));

// the classes borders
var negBorder=mean1.subtract(std2);
var posBorder=mean1.add(std2);

//create the layers
var imageNDVI=MyImage.select('NDVI');
var gtPOS=MyImage.gt(posBorder).selfMask().rename('PositiveBorder');
var ltNEG=MyImage.lt(negBorder).selfMask().rename('NegativeBorder');
var betMEAN=MyImage.gt(negBorder).and(imageNDVI.lt(posBorder)).selfMask().rename('MeanBorder');

var PositiveCOL = {
  palette: [
  '006622'
  ],
};

var NegativeCOL = {
  palette: [
  'ff0000'
  ],
};

var MeanCOL = {
  palette: [
  '00ff00'
  ],
};

Map.addLayer(gtPOS,PositiveCOL,'High');
Map.addLayer(ltNEG,NegativeCOL,'Low');
Map.addLayer(betMEAN,MeanCOL,'Mean');

I have try mosaic like this but got null image:

var mosaic = ee.ImageCollection.fromImages([gtPOS, ltNEG,betMEAN]).mosaic();
var PalMos={
  palette:[
    '006622','ff0000','00ff00']
};

Map.addLayer(mosaic,PalMos,'Mosaic');

but got this error:

  Mosaic: Tile error: Expected a homogeneous image collection, but an
  image with incompatible bands was encountered:   First image type: 1
  bands ([PositiveBorder]). Current image type: 1 bands
  ([NegativeBorder]).
            Image ID: 1 Some bands might require explicit casts.

My end goal: to merge those 3 different classifications into 1 layers, right now they are 3 independent layers.


texts: Calculating if moon is visible in mountainous terrain with GRASS GIS?

texts:
My goal is to produce a binary map indicating whether the moon is visible at a certain time of day with GRASS GIS. (i.e. A mountain could be blocking the view of the moon.) r.sunmask would be the perfect module. Pretend the moon is the sun and just provide the azimuth and altitude. However the module runs too slowly. r.sun is then the recommended substitute. However r.sun does not provide for direct azimuth/altitude input, just the time of the day. I could somehow do a funky backwards calculation of the time when the sun would match the azimuth/altitude of the moon, but I'm not sure that would work 100% of the time.

What is a good way to accomplish this task?

p.s. Is web Mercator (EPSG:3857) suitable for this task or do I need a different projection?


texts: Getting a Bounding Box in Latitude / Longitude using PostGIS shp2pgsql

texts:
Using the PostGIS shp2pgsql -s 4326 ... I generated the sql file and created my spatially enabled table.

When executing my query:

select box2d(ST_extent(the_geom)) from feature_data;

I'm returned with:

BOX(490755.1875 5460528,504640.4375 5466997)

My result appears to be in Easting and Northing, but I would like to return the bounding box in latitude / longitude. I have limited experience with GIS. How can I convert this result to the latitude / longitude coordinate system?


texts: Dates are off when opening netcdf file in R

texts:
I have a netcdf file with the time dimension that goes from 12/2007 until 12/2010 with monthly time steps. I can see the correct dates in both arcgis and panoply. Panoply tells me this file is a NetCDF-3/CDM, and reports

time double time(time=37);
      :units = "hour since 1950-1-1T00:00:00Z";
      :long_name = "forecast time";

But when I open the dataset in R using either brick(), ncdf, or ncdf4 I get these values in the time variable

>  v1<-get.var.ncdf(nc, "time")
> v1

 [1] 508008 508752 509472 510192 510936 511656 512400 513120 513864 514608
[11] 515328 516072 516792 517536 518256 518952 519696 520416 521160 521880
[21] 522624 523368 524088 524832 525552 526296 527016 527712 528456 529176
[31] 529920 530640 531384 532128 532848 533592 534312

and these is the values I get from the brick function

names       : X508008, X508752, X509472, X510192, X510936, X511656, X512400, X513120, X513864, X514608, X515328, X516072, X516792, X517536, X518256, ...
z-value     : 508008, 534312 (min, max)

Is there some parameter that needs to be changed that deals with time. I can't find any information rather than some sparse posts around the web.


texts: QGIS point displacement with features that each have 2 variables

texts:
My goal is to use point disclacement to have a view like this without overlapping symbols:

However, the result with using point displacement this (QGIS displays only the first attribute):

If I remove the cross in the first three categories (R1), the symbols that point in the other direction (R2) are shown on my map. Hence, the existing rules work if I can trick point displacement in displaying both directions.


texts: Creating custom CRS for India toposheets using QGIS?

texts: Indian toposheets (pre-2000 series) of 1:50k scale are in Polyconic projection, with Everest 1956 spheroid and India-nepal datum. How do I create a custom CRS in QGIS to represent this, so that I can then reproject the topo corners (tic points) to UTM-WGS 84 and then geo-reference the toposheet.


texts: Fill in empty space in raster image QGIS

texts:
I need a special projection to show the results of my research, but this projection is probably not defined in some points. How to fill them? I expect to fill at least the blank line at 180 degree meridian. Maybe be there is a way to somehow stretch the imagery or scale the projection?
CRS parameters:
+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
My results:


texts: Creating custom form that requires authentication when opening new project in QGIS?

texts:
I have a project in Qgis, and the layers I have stored in a Postgis database, I have created the credentials from the authentication menu and I have chosen the basic authentication method, I have only put user but not password to make it easier , but I do not accept the user (the credentials), you know what can be?

I introduce the data for a basic authentication, (user and password), of the PostGis connection. Then I save the changes and close the project. I reopen the project and then when I try to open the connection it asks me the master password and it gives me error and I get this message.


texts:
How do I limit spatial extent of features that is returned by arcpy.da.SearchCursor?  Can I use arcpy.env.extent?  Or somehow come up with sql statement to specify spatial relationship?   Or do I have to arcpy.SelectLayerByLocation_management to do intersection of my area of interest.

EDIT:
Here is what I want to accomplish.  Want to visit features from a feature class within a set extent.  I was hoping that second use of arcpy.da.SearchCursor respects arcpy.env.extent, but apparently it ignores it.

How can I do this (visit spatial subset of features)?  I want to have method with decent efficiency (e.g. avoid copying of features).

import arcpy
def proc():

    # I/O files
    fo = open('c:/users/yosuke/documents/log.txt','w')
    fc = arcpy.CopyFeatures_management('c:/Users/Yosuke/documents/ArcGIS/polytical_boundaries.mdb/tl_2010_us_state10', r'in_memory\fc')

    # grab everything, works fine
    arcpy.env.extent = 'MAXOF'
    n = int(arcpy.GetCount_management(fc).getOutput(0))
    m = len(list(arcpy.da.SearchCursor(fc, ('SHAPE@',))))
    extent = arcpy.Describe(fc).extent
    # got 52 from both, as expected
    fo.write('n: %s\n' % n)
    fo.write('m: %s\n' % m)
    fo.write('extent: %s\n' % extent)
    fo.write('\n')

    # want to grab only western USA
    arcpy.env.extent = arcpy.Extent(-120, 0, -90, 80)
    # this counts 27 western states, I suppose.
    n = int(arcpy.GetCount_management(fc).getOutput(0))
    # this still marches through all states, I want this to be spatial subset of USA
    m = len(list(arcpy.da.SearchCursor(fc, ('SHAPE@',))))
    extent = arcpy.Describe(fc).extent
    fo.write('n: %s\n' % n)
    fo.write('m: %s\n' % m)
    fo.write('extent: %s\n' % extent)
    fo.write('\n')

    # For now I am making copy of feature class
    # turned out that this is not very bad for my application, other processing
    # tasks I am applying on cursor rows are more computationally expensive.
    # I still want to know if this is standard way to accomplish this kind of
    # tasks (subset cursor by extent)
    fc2 = arcpy.CopyFeatures_management(r'in_memory\fc', r'in_memory\fc2')
    n = int(arcpy.GetCount_management(fc2).getOutput(0))
    m = len(list(arcpy.da.SearchCursor(fc2, ('SHAPE@',))))
    extent = arcpy.Describe(fc2).extent
    fo.write('n: %s\n' % n)
    fo.write('m: %s\n' % m)
    fo.write('extent: %s\n' % extent)
    fo.write('\n')


texts: Calculate new band for satellite image after predict values with RF

texts:
I have satellite image which contains 3 different bands.
I'm using Python (jupyter notebook) in order to calculate new band by applying random forest regression.
My problem is that after I have predict all the values for the new pixels, I don't know how to take it back to the original dataframe with the original bands in order to create in the end new image.

This is the process I did:

open the 3 bands image with rasterio, the band has this shape: (3, 869, 1202)
create pandas df when each row represents a pixel and each column is a band :

train the data and fit to random forest:

#split the data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

#import the algorithm
rf=RandomForestRegressor()

#reshape the y_train to fit the the model
y_train=y_train.values.ravel()

#fit the model
rf.fit(X_train,y_train)

rf_pred=rf.predict(X_test)

after checking the results, apply it to the full dataset and not only to training and set in order to predict the new band:

#create the data
data=df.iloc[:,1:]

pred_all=rf.predict(data)

#reshape to one column:
pred_all.reshape(1006560,1)

So after this, I don't know how to take this predicted values back to my table or to "link" it with the original pixel entities.

My end goal is to be able in the end to have this predicted values as new band so I can create image with the new predicted values.


texts: Controlling Overwrite Behavior in ArcGIS Python Script Tool Validation?

texts:
I'm creating a script tool to read the USGS earthquake GeoRSS feeds.

It works fine as a script with hardcoded values. The script also works fine as a tool for creating a new feature class for a given feed (hourly, daily, weekly possible choices).

The problem I run into is when I want to control the overwrite behavior. I want to give the user a choice: overwrite the selected feature class or append the selected feed to that feature class.

My tool has 3 parameters:

Output File
Feed Type: Hour, Day, Week
Overwrite: Boolean

In the tool validation script, I disable the overwrite until the output file parameter selects a feature class that already exists. The problem is my script won't run because of the default behavior of the tool validation.

My question is how best to overcome this problem. I want it to default to letting the user just select a possibly existing feature class and append records to that feature class. If the user checks the now enabled overwrite parameter, then it should delete that feature class and build a new one at its location.

Even if allowing the overwrite possibility seems like bad programming practice, I'm still facing the problem of having an automated updating to the feature class, which is what I'm trying to make default: e.g., if they select the week feed to an hour feed they just ran, then it should give the entire week, including the hourly records it already possesses, without duplicating them--my script already handles this process; the validation is the problem.

I'm thinking that using arcpy.env.Overwrite in the validation might help?


texts: Enlarging polygons only on two sides using ArcGIS Desktop?

texts: I want to enlarge a line of polygons only on the sides but I only found tools like "scale" or "buffer" which extend the whole polygon and not only widen them. The problem with scaling or buffering is that there's an overlapping of the polygons.


texts: Recast data from one geographic standard to another

texts:
I'm working right now a project looking at medical labour supply in rural areas in Australia. I aim to use STATA for most of this analysis, but a change in geographic standards by the Australian Bureau of Statistics has created a problem. To be more specific:

The program whose effect I am focussed on uses the old Australian ASGC 2006 Statistical Local Area (SLA) as a base
The ABS makes most contemporary data available using the newer ASGS Statistical Area level 2 (SA2) standard

Digital Boundaries are available for both of these on ABS Boundaries Online. I need to recast some of this contemporary data (Building Approvals, population, etc.) to the older standard. I am thinking of using a population grid (ABS Population Grid 2011) as a base on which to recast the data from one set of boundaries (SA2) to another (SLA2006). i.e. Look at the fraction of population in one SA2 that fits into the SLA, and transfer that much of the SA2's data to the SLA.

Any ideas on how I could do this in ArcGIS or QGIS?


texts:
I have two geoservers:
A has 1 workspace, 1 store (postgis) and 1 layer.
B is empty.
I copy the full data_dir from A to B, and reload config on B. But after reload only the workspace shows on B, not the store or the layer. If I restart B the store and the layer shows up.
Am I correct in copying the full data_dir?
Am I missing a part somewhere?


texts:
I try to display a wms-server in qgis. It serves several coordinate systems, but only few are working. Here is the link:
http://www.webatlasde.de/arcgis/services/Maps4BW/MapServer/WMSServer
For example EPSG:3857 is displayed, 31467 (which I need) is not
In ArcGis, the WMS is displayed correctly, so the server should be ok.
Can anyone confirm this problem or help me?


texts: Removing white space from field using field calculator

texts:
I am trying to remove extra spaces from values in a field programatically using a python script. This is the script I'm using:

field1 = "FULLST_NAME"
calc1 = "s.join(!" + field1 + "!.split())"
s = "\' \'"
arcpy.CalculateField_management(c_lines, field1, calc1, "PYTHON_9.3")

The expression used in calc1 works fine when I manually input it in the field calculator, i.e.

"   S STATE    ROUTE " is changed to "S STATE ROUTE"

However, when executing the script, it outputs:

"S' 'STATE' 'ROUTE"

with the ' ' between each string. I'm guessing there is some syntax I'm missing when storing the calculation in the variable and using it as a parameter in CalculateField_management(). How can I fix this?


texts: Add Road IDs to Parcel Layer (Or Parcel IDs to road layer)

texts:
I am trying to create a table that links a road segment with all the parcels that would be "touching" that road.
A road has an average of about 10 parcels, but most parcels would only be on one road. However, a parcel on a corner would be "touching" 2 road segments and a rural parcel may actually be on 3 road segments (front, back and one side).
My road is a road centreline, so they don't overlap, but a 15-20m buffer would be perfect, I just don't know how to add\join the fields after that. It doesn't matter whether the road id is added to the parcel or the parcel ID is added to the road, in the end I just want a table with both IDs.

I am using ArcGIS 10.3.1

For the example above, the table I want would be:

P   R
1          A
2          A
3          A
4          A
4          B
5          A
6          A
7          A
8          A
8          B


texts:
I have a quadkey function to try to reach a non-Bing Source which uses a quadkey but it returns 404s everytime although the request looks correct. I wondered instead if I should try to pass a custom url to an ol.source.BingMaps layer?
I've looked at the api and found the getTileUrlFunction/setTileUrlFunction but can't find any examples to demonstrate how to use it. My original question was kindly answered by Alvin Lindstam but I haven't yet had any success.

I'm wondering if the getTileUrlFunction is the right approach or whether I should try to use the uncompiled version of Openlayers and try to extend the base class?
My javascript isn't really up to this and I wouldn't really know where to start but if necessary could I just duplicate the BingSource, give it a new name and insert a different URL?

From the base directory, the path is as follows: /50k/1202/0031/3033/12020031303302.ipic
My function therefore returns this:

return baseURL + z + "/" + quadkey(x, y, z) + ".ipic";

I'm sure this may not make much sense so I've created a fiddle to try to demonstrate. These tiles all work perfectly in an old Silverlight app but I just can't get them to load and my client is non-technical to say the least!


texts: Exporting accuracy in Google Earth Engine

texts:
I am trying to export the Number trees with the corresponding accuracy, here is the code I am using:
https://code.earthengine.google.com/33277d8909f0cead0d1ef641564cd5d5
I am trying to achieve something like this:
id: 0
geometry: null
properties: Object (2 properties)
Accuracy: 0.8813559322033898
numTrees: 5
id: 1
geometry: null
properties: Object (2 properties)
Accuracy: 0.8713559322033898
numTrees: 10


texts: How to access ArcGISDynamicMapServiceLayer with credentials (token)

texts:
I know I simply need to append a token=XXX to the export query string but how, in code, do we initialize "esri/layers/ArcGISDynamicMapServiceLayer" to use credentials?  In my case I only want to provide a token value.

Here's my code to show how I've hacked around this...probably good enough for my needs but seems wrong.  Should I be creating a Credentials instance and associating it with the layer?

function createDynamicLayer(options, map) {
    var dataset = options.dataset;
    var query = ioQuery.queryToObject(acme.right(dataset.CacheURL, "?"));
    var imageParameters: esri.layers.ImageParameters = new esriImageParameters();
    var credentials: esri.Credential;

    imageParameters.format = query.format || "png";

    var url = acme.left(dataset.CacheURL, "/export");
    var layerName = acme.left(acme.right(dataset.CacheURL, "/export/", true), "?");
    if (layerName) {
        layerName = layerName.split(",");
        layerName = array.map(layerName, "return parseInt(item, 10)");
        lang.mixin(imageParameters, {
            layerIds: layerName,
            layerOption: esriImageParameters.LAYER_OPTION_SHOW
        });
    }

    if (query.token) {
        credentials = new IdentityManagerBase.Credential({ token: query.token });
        url += "?" + ioQuery.objectToQuery({ token: query.token });
    }

    var layerOptions = lang.mixin({}, lang.getObject("layers.esri.dynamic", esriConfig), {
        id: options.id,
        imageParameters: imageParameters
    });

    //Takes a URL to a non cached map service.
    var dynamicMapServiceLayer = new ArcGISDynamicMapServiceLayer(url, layerOptions);
    return dynamicMapServiceLayer;
}


texts: How can I display a spatialLite view in a custom data sheet with print composer

texts:
I created a spatiaLite Database with a point layer (every point is a tree). Additional I created several tables to fill the columns of the spatial layer with attributes using the editing widget “value relations”. I need this because I want to speed up data entry! For some other attribute fields the data entry occurs via “line edit”. I also designed a custom form like Nathan Woodrow describes it in his blog. This works pretty good.
However, for output I want to create several custom data sheets with the print composer - one for every tree. Because I want to display the values from my linked tables and not the common keys I created a view that is composed of all attributes I need inclusive the geometry column (geom) of the spatial layer.

 CREATE VIEW tree_attrib1 AS SELECT baumkataster.ROWID AS ROWID,
 baumkataster.baum_ID, baumkataster.geom AS geom, gattung.gattung AS gattung, art.art AS
 art FROM baumkataster JOIN gattung ON baumkataster.gattung = gattung.gattung_ID
 JOIN art ON baumkataster.art = art.art_ID;

To use the atlas generator in the print composer I need to “convert” my view into a point layer.

My questions now are:
1: How can I display my view as an point layer in qgis?
2: A view is stored only temporary. Will this new point layer be updated when I insert new features (trees)?


texts: Auto-Increment Integer Field in a Feature Class

texts:
Does anyone know of a good way of auto-incrementing an integer field within a SDE feature class when the feature is created?  I don't want to have a script run afterwards to populate the field and I'd like this to be run when a feature is created (insert).

Is there a trigger that could be run or can the Identity property be used?

Just wondering if anyone has a good solution for this?

ArcGIS Version = 10.2.1

ArcSDE = 10.1 sp1

RDBMS = SQL Server 2008 R2

Thanks in advance!


texts: Create polygons (representing boundaries) from lines

texts:
I received a file of municipality boundaries, where the boundaries are represented as lines.  I would like to convert this to a set of polygons, each representing a municipality.  In other words, from a set of lines, I would like to create the highest number of non-overlapping polygons.

I would like to do this in R (currently using the sf package, but could use another).  I believe the correct command is st_polygonize, however the data does not seem to be extremely high quality, so some of the lines do not exactly touch (for one pair, I calculated the distance between the endpoints as 1.77e-15 meters).  There are also holes in the municipalities.

My questions:

Is st_polygonize the correct command?  I only have one observation for each line segment, but a single line segment can belong to multiple municipalities.  Whoever created this file seems to have done it in a way so that each line segment belongs to at most 2 municipalities, but I do not know how to verify this.
How do I deal with the fact that not all line segments will touch?  st_polygonize does not seem to have a tolerance parameter, is there some other command I can use?

This is the file as I received it:

And here is a closeup, with each line segment in a different color to give an idea:

EDIT:

I solved the tolerance issue via the following:

Use st_line_sample(line, sample = c(0,1)) to find the endpoints of each line.  I then calculated pairwise distances between all endpoints using st_distance and filtered for distances that were greater than 0 and less than a tolerance number.  The next part requires some manual work, but using the filtered distances, one can identify which line segments are misaligned.  I then used st_snap to snap these line segments to each other.  With the new snapped line segments, I can then create the polygons as per Spaceman's answer below.


texts: Create a Dynamic Table in your Map Layout - how to use Graphic table element?

texts: Is this Solution (Create a Dynamic Table in your Map Layout using  graphic table element)  is part of this Extension ESRi production-mapping extension ? how can i use graphic table element ?


texts: How can I convert Military Grid Reference System (MGRS) to lat/long or UTM?


texts: Question about Datum not recognised by GRASS

texts:
I am now projecting the dataset to CGCS2000 / 3-degree Gauss-Kruger CM 120E.
But the dataset could not be recognised by GRASS when I trying to import the DEM data.  The warning message kept popping up:
WARNING: Datum <China_2000> not recognised by GRASS and no parameters found
And the output from g.region -p is:
projection: 99 (CGCS2000 / 3-degree Gauss-Kruger CM 120E)
zone:       0
datum:      ** unknown (default: WGS84) **
ellipsoid:  grs80
north:      1
south:      0
west:       0
east:       1
nsres:      1
ewres:      1
rows:       1
cols:       1
cells:      1
how can I fix it?


texts: Can Draw and Edit mode be enabled at the same time on Leaflet?

texts:
Imagine you want to draw several shapes (marker, linestring, polygon) and also edit them after the shape is drawn or whenerver you want having the possibility to continue adding stuff.

at the moment i have to draw everything and then switch to edit mode then switch to draw if i missed something and so on. Are they mutually exclusive?

I know for sure that this is possible with openlayers.


texts: Is there a way to use old Print composer templates in 2.0?

texts:
As the title says, the old templates do not work in 2.0.  It started with having to redraw boxes, now it causes it to crash completely.

Is there an easy fix.

Can i download templates for use??


texts: How to cast Geometry to Geography in PostGIS/Postgres with sqlalchemy + geoalchemy

texts:
Some PostGIS functions behave differently based on the arguments types.

Example: ST_Distance will calculate projected units distances for geometry arguments and meters distances for geography arguments.

In PostgreSQL you can easily cast with the :: syntax, but what if you are using sqlalchemy objects or geoalchemy?


texts: Selecting elements from network with shortest path in QGIS?

texts:
I have a network of vector elements (let´s say 70 elements) that are connected in each node. Actually, it is a pipe network built with the help of GHydraulics (http://epanet.de/ghydraulics/index.html.en).

I want to select the list of elements of the network layer using a shortest path tool like the one of road graph. The order of the list matters - it should follow the path.
Any clues?

Even the order could not be achieved; I tried to use spatial query with:

source: network layer
reference: shortest path layer
where features: overlaps

There should be at least 20 elements from source on the shortest path, but it does not select anything. With a buffer it works, but I have to review the buffer distance to check that no small falls within the buffer by mistake.


texts:
First of all I am not a Database specialist. I am using MSSQL to store my SDE 9.3 Data which is served with AGS 10.0. I have Vector and raster data stored in the databases and have very little room left on my SQL Server because the logfiles are just very big. One Database which contains 314GB of Data has a logfile which has a size of 28GB!

What can I do (if anything) to reduce the size of this logfile?


texts: Foreign key value in destination table is nullified when origin feature is deleted from ArcGIS relationship - workaround?

texts:
I have an SDE geodatabase in which I have an origin feature class with relationship classes to several destination feature classes and tables. When a feature is deleted from the origin table, the foreign keys in the all of the destination feature classes and tables are nullified. This has been annoying for the destination feature classes because the foreign keys have to be reestablished, however, it has resulted in lost data for the destination tables because we have no way to figure out which feature they originally belonged to (because there are no spatial data attached) and this has created serious issues.

Is there a workaround for this or a different setting that I need to consider when setting up my relationship classes?


texts: QGIS shapefile layers to database geometry type

texts:
I have opened Arcgis shapefile with QGIS. And I used QGIS database manager tool to create table and import data to postgresql database. Database manager is creating and importing success. But there is an issue about geometry types.

Shapefile Geometry         Postgis Table Geometry

 - Polygon         ==>     MultiPolygon
 - Polyline        ==>     MultiPolyline
 - Point           ==>     MultiPoint

Polygon type is creating as Multi on postgis table. What do you think about this?


texts: Disconnecting ArcGIS Pro from PostgreSQL database?

texts:
What is the equivalent to the ArcCatalog's Disconnect option in ArcGIS Pro?  I have successfully connected to a Postgresql database, but cannot figure out how to disconnect.

This is what I'm seeing on the right-click of my database connection:

I know the ArcGIS Pro client is still connected to the postgres database:

I also tried:

Removing the connection.
When that didn't work, I opened a blank
project.

Neither closed the connection.

Update 2:
I don't think this would matter to the connection behavior, but there's no sde/gdb in the postgresql database--just straight 9.3 postgresql database.  Thought I would mention in case others stumble upon this question.


texts: 2dm data disappears after rotating layout in QGIS

texts:
I ran into a bug (I suppose) recently. After loading a *.2dm format mesh (triangulated irregular network with elevation data), with the Layer --> Add Layer --> Add Mesh Layer command, I put some data on it (2D water flow vectors and absolute speed to colorize it) and after created a Layout with a Map to export as pdf.
However, when I alter the Map Rotation property on a Map item in the Layout, the Mesh data disappears. So basically any rotation other than 0 degree makes the vectors, and colour disappear (it still shows the triangulation).
It shows again immediately after I put it back to 0 rotation.
Any ideas what's happening here?


texts: Calculating viewing angle (2D) between polygons and points

texts:
I am looking for a solution to the following problem.

I have several points (for example houses) and several polygons (eg forests).
Seen from the houses, I can see the forests within a certain angle. This angle is different for each house. To calculate this angle, I need the points on the polygon where the angle of view is maximum. ie the extreme left point of the forest and the extreme right point of the forest. I could calculate the angle with these two points.
This is about 1,500 points (houses), so that the drawing of lines and the coordination of the endpoints is unfortunately not an option.

I hope the drawing shows what I'm looking for.

Is it possible to determine the coordinates of the two outermost visible points of the forest for each building?

How I can determine the angles?

I use ArcGIS Desktop with a Standard license without any extensions.


texts: Is TOA reflectance in floating point or DN?

texts: After  converting at radiance to TOA reflectance,is going to be in DN or floating point.


texts: Geoserver ImageMosaic ignores existing datastore.properties (connect to PostGIS) and automatically generates index using shapefile

texts:
I am new to Geoserver. I am trying to publish a WMS service for a raster time series using the ImageMosaic plugin in Geoserver. I am using the sample dataset and instructions at this link

I want to use postgres/postgis to map the time values to corresponding layers and created datastore.properties

SPI=org.geotools.data.postgis.PostgisNGDataStoreFactory
host=localhost
port=5432
database=postgis_sample
schema=public
user=postgres
passwd=qwer1234
Loose\ bbox=true
Estimated\ extends=false
validate\ connections=true
Connection\ timeout=10
preparedStatements=true

But when I attempt to create a new store -> Raster data - ImageMosaic -> choose data directory and save, Geoserver ignores my datastore.properties file and uses the default settings by creating shapefile and automatically rewrites and generates a new datastore.properties file. This should not happen according to the documentation. The postgres/postgis service is running and I am able to connect and write from the windows commandline. I also grant all permissions to the public schema. Could someone point out what is going wrong here?

I am using Geoserver 2.11.4, PostgreSQL 10 with PostGIS 2.4.4.

indexer.properties

Caching=false
TimeAttribute=ingestion
Schema=*the_geom:Polygon,location:String,ingestion:java.util.Date
PropertyCollectors=TimestampFileNameExtractorSPI[timeregex](ingestion)

timeregex.properties

regex=[0-9]{8}

EDIT:
I have attached the log entries after following the above procedures in this pastebin link.


texts: Google Maps API: How to remove terrain from the satellite view?

texts:
(I asked this in SO General - but then thought it might be better here - sorry for the cross-post!)

Using the API, is it possible to turn off the 3D terrain in the Google Maps satellite view? I'd like to get the imagery without it being warped by the 3D surface.


texts: How to make an opacity slider for a group of layers in Google Earth Engine

texts:
I am building a ui in Earth Engine by following this example.. https://google.earthengine.app/view/forest-change
The example code can be found here.. https://code.earthengine.google.com/b02a58539a6374cdebd9416729a08f47
The snippet of code that influences the output of the opacity slider can be seen below...
opacitySlider.onSlide(function(value) {
  mapPanel.layers().forEach(function(element, index) {
    element.setOpacity(value);
  });
});

This function subsequently influences the opacity of every layer in the map panel.
I have additional functionality in my application that adds other boundaries/imagery that I would like to keep rendered on the map, even if the opacity slider is active.
Is there an easy way to specify only a select group of layers in the map panel to work with the opacity slider?
I am aware that there is a default layer panel with individual opacity sliders, but I'd like to keep that control hidden from the app.


texts: ModuleNotFoundError: No module named '_gdal' on mac machine

texts:
I have installed gdal in the following way
python3 -m pip install --global-option=build_ext --global-option="-I/usr/include/gdal" GDAL=='gdal-config --version'
If I try to call it
from osgeo import gdal

I get the error
ModuleNotFoundError: No module named '_gdal'


texts: How to load GPX tracks in graphical modeler as input

texts:
Is there a way to load tracks from GPX into graphical modeler as vector input?
When you open GPX in QGIS, it comes up with dialog listing layers:

route_points
routes
track_points
tracks
waypoints

It seems to me in graphical modeler only first layer route_points is getting loaded. So I guess in graphical modeler for vector input you need to specify you want to load 'tracks' and not other layers somehow?


texts: converting coordinates from unusual DMS format to decimal

texts:
So I've been shared a dataset with what to me is an unusual format for degree min second data.
In their own words:

    The coordinates provided are longitude and latitude in degrees minutes and seconds.
    So for example, coordinates recorded in the spreadsheet as
    X_Long: -10.0135 Y_Lat:  53.5726, can be taken as Longitude: 10º 01'35" W,
    Latitude: 53º 57'26" N using standard notation.

I need to get this into decimal. Any quick way to do this ?

Edit:

in trying to solve this I made a following code:

    # Fake data
    X_Long <- as.numeric(c(-2.4564, -1.0134, -8.01))

    long.temp <- abs(X_Long - trunc(X_Long))
    long.min <- trunc(long.temp*100)
    long.sec <- trunc(as.integer((long.temp*100 - long.min)*100))

Its giving me weird rounding errors:

    > X_Long
    [1] -2.4564 -1.0134 -8.0100
    > long.temp
    [1] 0.4564 0.0134 0.0100
    > long.min
    [1] 45  1  0
    > long.sec
    [1] 63 34 99

This seems to be coming from the trunc function - compare results here

    > long.temp*100
    [1] 45.64  1.34  1.00
    > trunc(long.temp*100)
    [1] 45  1  0

Anyone know why this is happening ?


texts: Dictionary.set: Parameter 'value' is required

texts:
I am trying to calculate the area of each class of my grid. My code reference is from here.
My code: Earth engine link
var grid = ee.FeatureCollection("users/miaohu828/finalGridexpot1000")
var modisLandcover = ee.ImageCollection("MODIS/006/MCD12Q1")
var bp = ee.FeatureCollection("users/miaohu828/bankment_polygon")

var filtered2018 = modisLandcover.filter(
  ee.Filter.date('2018-01-01', '2018-12-31'))
// The filtered collection has 1 classified image for 2018
// We can call .first() to get that element and cast it to an image
var landcover2018 = ee.Image(filtered2018.first())
// Each year's image contains multiple bands with different classification scheme
// We pick the band 'LC_Type1' for this example and extract it.
var classified2018 = landcover2018.select('LC_Type1')
// There are 17 classes in the scheme as described in https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD12Q1
// We pick the color associated with each class from the documentation and create a palette
var palette = ['05450a', '086a10', '54a708', '78d203', '009900',//forest
               'c6b044', 'dcd159','fbff13', //shrub, savannah
               'b6ff05', //grassland
               '27ff87', //wetland
               'c24f44','ff6d4c',//croplands
               'a5a5a5',//urban
               '69fff8', //snow
               'f9ffa4', //Barren
               '1c0dff'//water
               ].join(',');

var vis_params = {min: 0, max: 17, palette: palette};
// Visualize the classified image
//Map.addLayer(classified, {min:1, max:17, palette: palette}, 'MODIS Landcover 2018')
var bpLandcover2018 = classified2018.clip(bp)
Map.addLayer(bpLandcover2018,{min:1, max:17, palette: palette},'Yangtze Land Cover 2018')

// Add a new property 'grid_name'
var finalGrid = grid.map(function(f) {
  return f.set('grid_name', ee.String('grid_').cat(f.get('system:index')))
})
print(finalGrid.first())

var calculateClassArea = function(feature) {
    var areas = ee.Image.pixelArea().addBands(classified2018).reduceRegion({
      reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'class',
    }),
    geometry: feature.geometry(),
    scale: 500,
    maxPixels: 1e10
    })

    var classAreas = ee.List(areas.get('groups'))
    var classAreaLists = classAreas.map(function(item) {
      var areaDict = ee.Dictionary(item)
      var classNumber = ee.Number(
        areaDict.get('class')).format()
      var area = ee.Number(
        areaDict.get('sum')).divide(1e6).round()
      return ee.List([classNumber, area])
    })

    var result = ee.Dictionary(classAreaLists.flatten())
    var grid_name = feature.get('grid_name')
    return ee.Feature(
      feature.geometry(),
      result.set('grid_name', grid_name))
}

var districtAreas = grid.map(calculateClassArea)
var classes = ee.List.sequence(1, 17)
var outputFields = ee.List(['grid_name']).cat(classes).getInfo()

print(districtAreas)

Export.table.toDrive({
    collection: districtAreas,
    description: 'class_area_by_grid',
    folder: 'earthengine',
    fileNamePrefix: 'class_area_by_grid',
    fileFormat: 'CSV',
    selectors: outputFields
    })

My code runs with an error:
FeatureCollection (Error)
Error in map(ID=000000000000000000ab):
Dictionary. set: Parameter 'value' is required.

I can't find a way to solve it.


texts: Openlayers 3 connecting to a wms feed with a different projection

texts:
I am trying to incorporate a wms feed in the British national grid projection (EPSG:27700) into a openlayers 3 project. I initially had used this wms layer in an openlayers 2 application using the following syntax:

new OpenLayers.Layer.WMS(
        "ukimage",
        "http://www.linktowms",
        {
            "LAYERS": 'ukimage',
            "STYLES": '',
            format: state.image_format
        },
{
    buffer: 0,
    displayOutsideMaxExtent: true,
    isBaseLayer: true
}
);

But I am having trouble converting this over to openlayers 3. I first tried setting the map projection to British national grid (view: new ol.View({ projection: 'EPSG:27700'}) as I had in the openlayers 2 app, but this produced the error

TypeError: h.B is not a function

and a white map.

When I try and use proj4 to convert the projection to EPSG:3857 as described here: http://openlayers.org/en/v3.15.1/doc/tutorials/raster-reprojection.html I get the same error:

TypeError: h.B is not a function

and again a white map.

my wms request is generally a varient of:

new ol.layer.Image({
      source: new ol.source.TileWMS({
        url: 'http://www.linktowms',
        params: {'LAYERS': 'ukimage'},
      })
    }),

(I tried adding projection but that didn't help)

I am trying to work out if the image is even getting loaded from the WMS feed, but the syntax seems to have changed significantly from openlayers 2, instead of the properties and functions that have rather descriptive names, there are a now lots of short functions with names like Ua and Ig, which leaves me a bit in the dark as to what's actually going on when I try and inspect the layer with firebug.

I don't mind what projection the application is in, as I can re-project my other layers to match the basemap, so basically any way of using this EPSG 27700 WMS layer in openlayers 3 would be great.


texts: How to present animated movement of a dam?

texts:
I'm want to know how I can represent planimetric and altimetric movement, measured at different times, with animated cards or video.

Which software do I need? How do I go about doing this?


texts: Mapping a feature from the attribute table in QGIS?

texts:
I currently finished adding all the needed attributes for my research. I would like to map these attributes. Is there any way to do so? The original method I was taught was to select:  on the map and then add the attributes individually.

However, I added the attributes first by mistake. As in when I was editing the attributes table, I simply selected the "add feature" option. So when I select a row and click "zoom to row on map" it says "geometry is NULL". So is there a way for me to choose a row (each row is a location of a city) and map that row on the map?


texts: Is project more Leaflet or OpenLayers or Map Server friendly?

texts:
I have several tables in a database on PostGIS.

I want to create a website where the users can log-in, show the layers they want and filter the elements of the layers based on the fields of the layers.

I want to create a very friendly interface for the users to filter the layers.
I mean just boxes to check and limits of interval to set.

Based on this description, I would say that OpenLayers is better but on the other hand, the numerous community around Leaflet makes me think that the support is better with Leaflet.

Moreover, I just started learning HTML/CSS and I know barely nothing about JavaScript but learning is not trouble for me. I love it and I will need to improve my web development skills anyway (for another project).

I have seen Choosing OpenLayers or Leaflet? but it dates from 3 years ago at least and I am wondering if Leaflet hasn't solved its cons. For example OpenStreetMap changed from OpenLayers to Leaflet.


texts: MapInfo Error - 'Invalid distance units: .'

texts:
I am receiving an error (Invalid distance units) when attempting to import a MapInfo MIF File (.mif) into MapInfo.

Has anyone come across this error before?

Can't seem to find much about in the MapInfo documentation or online.

Opening a similar file in Notepad, where I get the same error, shows the following:


texts:
Is there a way to include animate properties of a vector layer (GeoJSON)?

var highlightStyle = new ol.style.Style({
    stroke: new ol.style.Stroke({
      color: '#fff',
      width: 1
    }),
    fill: new ol.style.Fill({
      color: 'rgba(255,0,0,0.1)'
    })
  });


texts: Calculate class area after classify random Forest image

texts:
After classify using randomForest my image

var imagemClassificada = addBands2018.select(bands).classify(classificadorTreinado)

I would like to calculate the area for each class, vegetation and pasture.

var area_classifyed = imagemClassificada.multiply(ee.Image.pixelArea())
                      .reduceRegion(ee.Reducer.sum(),VEGETACAO,30,null,null,false,1e13)
                      .get('constant')
print('VEGETATION AREA: ', ee.Number(area_classifyed).divide(1e6))

Focus that I use the FeatureCollection of 'VEGETACAO' to calculate this.

But, after all, I get some error:

  Provided visualization parameters are not an object: geometry
      in , line 3515

This is the correct form to calculate the area em Earth Engine?
HERE goes my code


texts: My Google Maps add Things to do, resturants, etc

texts:
Search high and low for how to add "things to do" or restaurants to a custom map I created using Google MyMaps.

Ex: Take the results of this link
https://www.google.com/maps/search/restaurants/@42.2557654,-86.3962214,9.71z

And add it as a layer to MyMap. I can search within mymaps for locations and add one by one, but I want to leverage what the results of google maps to show me all the restaurants in an area.

Backstory: I have 10 or so houses on a mymap, trying to overlay things to do in order to determine which house to buy :)


texts:
I would like to know is it possible to read or import a CGM file using QGISand also export other file format back to CGM using QGIS.(https://en.wikipedia.org/wiki/Computer_Graphics_Metafile)
Is there any existing plugin for this or if we need to implement a new plugin is is it possible using QGIS


texts: Python command for QGIS that iterates over .xyz files doesn't work

texts:
I have a small Python script for QGIS that goes to a folder and iterates over its files and uploads all .xyz files. This is the script:

import glob, os

dir = r"C:\Users\User\Desktop\xyz_daten"

os.chdir(dir)

for file in glob.glob("*.xyz"):
    uri = "file:///" + dir + file + "?delimiter=%s&crs=epsg:4647&xField=%s&yField=%s" % (" ", "field_1", "field_2")
    name = file.replace('.xyz', '')
    lyr = QgsVectorLayer(uri, name, "delimitedtext")
    QgsMapLayerRegistry.instance().addMapLayer(lyr)

unfortunately I get the following error message right under the first line :

> Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 63, in runsource
    code = self.compile(source, filename, symbol)
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 168, in __call__
    return _maybe_compile(self.compiler, source, filename, symbol)
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 99, in _maybe_compile
    raise err1
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 87, in _maybe_compile
    code1 = compiler(source + "\n", filename, symbol)
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 133, in __call__
    codeob = compile(source, filename, symbol, self.flags, 1)
  File "<input>", line 6
    import glob, os
         ^
SyntaxError: invalid Syntax

but where is the syntax error and how would it be correct?


texts: Automatically reloading layer when source is changed

texts: I have a ol.layer.Vector that regularily changes its source. How do I perform a reload of the source everytime I change it?


texts:
I want to combine categorized symbol data for the sake of symbology. For example: Let's say I classified a data column of a shapefile and it gave me 20 different values. Ten commercial building types and ten residential building types. In Arcgis I could very easily combine that to two symbols, one for residential and one for commercial.

Running QGIS 2.18, MacOS Sierra 10.12.2


texts: Performing table look-up using Python in ArcGIS Desktop?

texts:
Using ArcGIS 10.1 the arcpy module has a TableToNumPyArray command, which converts a table into an array.

Unfortunately I don't have access to 10.1 on my work PC. Is there a workaround which will accomplish this using ArcGIS 10.0?

Background

I'm iterating through features in a featureclass, and I need to look-up a price depending on an attribute value.

I need to find a price for multiple fields per feature, so I can't use a table join. I can't create a relationship class as I only have an ArcView license.

A messy workaround could be to load the table into a searchCursor and iterate through each record until I find the desired value. But this would be intensive as I'd need to work through the searchCursor for every feature in the featureclass.


texts: Drape - z value to raster - add a z column

texts:
I'm trying to use the drape set z value to raster function in qgis. I have a shapefile of points with the coordinates x and y and a DEM raster. I would like to give the z value to the points by taking it from the raster. When I apply the drape function, it generates a new shp file but without any new z field.
How can I resolve it?


texts: WMTS in OpenLayers unable to zoom into detail

texts:
Based on the OpenLayers example at
https://openlayers.org/en/latest/examples/wmts-layer-from-capabilities.html?q=wmts

I'm trying to load a WMTS in my own example here
https://codepen.io/vishnu4/pen/EMBZdp?editors=1010 , but it really only seems to show just one large image that I can't zoom into to see more details.  If I load it in QGIS, there are all sorts of other features that I'm not seeing.

How do I zoom into this WMTS to show the QGIS level of detail instead?


texts:
I'm using QGIS 3.4 and Python 3.7.0 on Windows 7.

From the Python console of QGIS, I created a test script ("test.py") which contains the following line:

print('Hello)

and this is saved in C:\OSGeo4W64\bin\test.py

I want to run this script from outside QGIS.

I tried opening "python-qgis.bat" which I found in C:\OSGeo4W64\bin, and typed:

test.py

but got an error saying [NameError: name 'test' is not defined].

The code itself works fine from QGIS, but because I want to run it from outside QGIS, do I need all the import... etc at the beginning of the script?

What do I need to do?


texts: Error trying to upload data using Mapbox API

texts:
I'm using the Mapbox API and I'm trying to upload my GeoJSON but I'm getting the following error:

  bounds west value must be between -360 and 360

My GeoJSON looks like this:

{"type":"FeatureCollection",
  "crs":
    {"type":"name",
      "properties":
       {"name":"EPSG:27700"}
     },
  "features":
    [{"type":"Feature",
      "id":1,
      "geometry":
        {"type":"MultiLineString",
          "coordinates":[[
             [343131.34379999992,566499.875],
             [343129.625,566502.3125],
             [343126.1875,566505.9375],
             [343124.6875,566507.3125],
             [343122.875,566508.125],
        ...

So I noticed that this is in BNG and I assume that Mapbox needs that in lat long values.
I did a bit of research and I think that I've to convert from EPSG:27700 to EPSG:4326.

Is there anyway to convert the file in a batch process? My file is about 500mb and contains few millions of records.


texts: OpenLayers : how to dynamically change projection of a Tile Layer?

texts:
I've added a tile base map and and a couple of vector layers using OpenLayers. I'm implementing a change map projection functionality. For that, I need to change map projection and all added layers .

changeLayerProjection(layer, oldProj, newProj) {
    if (layer instanceof olLayerGroup) {
      layer.getLayers().forEach(function(subLayer) {
        this.changeLayerProjection(subLayer, oldProj, newProj);
      });
    } else if (layer instanceof olTileLayer) {
      var tileLoadFunction = layer.getSource().getTileLoadFunction();
      layer.getSource().setTileLoadFunction(tileLoadFunction);
    } else if (layer instanceof olLayerVector) {
      var currProj = layer.getSource().getProjection();
      var features = layer.getSource().getFeatures();
      if (currProj === null) {
        for (var i = 0; i < features.length; i += 1) {
          features[i].getGeometry().transform(oldProj, newProj);
        }
      }
    }
  }

But for the tile layer, I'm doing like this:

 if (layer instanceof olTileLayer) {
          var tileLoadFunction = layer.getSource().getTileLoadFunction();
          layer.getSource().setTileLoadFunction(tileLoadFunction);
        }

getTileLoadFunction
will only return that function. I want to change tile layer projection as well. What is best way to do this?


texts: Counting consecutive pixel values for set of rasters using ArcGIS Spatial Analyst?

texts:
I'm using ArcGIS 10 with Spatial Analyst.

I have a set of rasters (8 in total) that only every contain a 1 or a 0 for each cell. Each raster represents a different years worth of data. For arguments sake year 1 to year 8.

I can add all the rasters together which will give me a final grid with values ranging from 0 to 8. An 8 indicating that the cell was constantly 1 for the set of rasters (all years).

I would like to find out for each cell the longest consecutive number of 1's.

So for example the total grid may record for a single cell a value of say 5 but over the 8 grids that cell has the largest consecutive number of 1's equal to 3. Or another way of expressing this is for 3 years that cell was a 1 then it started to oscillate between zeros and ones.

My raster processing skills is not as hot as my vector processing skill and I have had a good look at the ESRI help file but I can't figure out how one would achieve this using off the shelf geo-processing tools?

Any ideas?


texts:
I get an Runtime Error using my QGIS-Plugin after inserting following code

if mapType == 0:
    QgsMapLayerRegistry.instance().addMapLayer(layer1)
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer2)

The entire Error-Message is:


texts: Calculate distance covered by points in a shapefile

texts:
Is there a way to calculate the total distance covered by thousands of points in a shapefile in ArcMap 10.4 or possibly ArcPro? I am trying to apply a buffer for each point, then dissolve and find the centerline of each dissolved polygon per How to convert points to lines?. Then find the distances of the lines. But I'm getting a Distance Value required error message:

I have also tried this, but I don't have the equivalent fields for ROUTE1 or MEASURE, so I can't get it to run.


texts: Defining coordinate reference system with rotation in GeoServer?

texts:
I am using GeoServer and have a layer in EPSG:900913 ("Google Mercator").

I need to "rotate" the map around certain point (say, 1500000, 7000000) by certain degree (say, 30 degrees clockwise). How could I define such a coordinate system based on EPSG:900913?

GeoServer's angle vendor option does not work for my purposes as I need to tile the map later on.

As far as I understand this, my only option is to define an own coordinate system. For GeoServer I'd need to define it in WKT form. The configuration seems to be straightforward, but I have a difficulty defining my rotated CRS in WKT.

I am wondering how to apply a rotation around certain point onto a CRS like Google Mercator:

PROJCS["WGS84 / Google Mercator",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Longitude", EAST],
    AXIS["Latitude", NORTH],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Mercator_1SP"],
  PARAMETER["semi_minor", 6378137.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["central_meridian", 0.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["x", EAST],
  AXIS["y", NORTH],
  AUTHORITY["EPSG","900913"]]

My questions, specifically:

How to write a WKT which transform an existing CRS? My guess would be that I need a new PROJCS wrapping an existing one and adding a PROJECTION clause.
How would I found out the projection id (like Mercator_1SP above) and the required parameters (the PARAMETER clauses)?
Can I "reference" EPSG:900913 in CRS WKT instead of copy-pasting the whole PROJCS clause?


texts: Opening GeoTIFF (georeferenced in ArcGIS for Desktop) in QGIS?

texts:
I tried to open a georeferenced geotiff created in ArcGIS for Desktop in QGIS.

I could open them, but they are not in the right spot.

These raster datasets contains .aux .rrd .tfwx  .tif.aux.xml  .tif.xml and .tif files.
I can open them in ArcMap without problem.

Am I missing something?


texts: How to retrieve altitude from pixel value in a PNG file requested from GeoServer

texts:
I'm using a GeoServer 2.18 from Unity for terrain rendering. I'm asking the GeoServer this kind request:
{0,1}wms?service=WMS&version=1.1.0&request=GetMap&layers=Elevation&bbox={3,1},{4,1},{5,1},{6,1}&width={2,1}&height={1,1}&srs=EPSG:900913&format=image/png

on elevation data. I've got a nice black and white PNG file but how can I convert the pixel value into a real altitude?


texts: Problem with aes and ggplot2 while plotting gwr/aes length

texts:
Having trouble with the following code while doing a Geo-regression

MF <- read.csv("MF.csv", sep=",")
modelo1 <- lm(MF$MFALC_15_2 ~ MF$MFEFF2015+MF$agro1315+MF$MEDALC+MF$density)
summary(modelo1)
plot(modelo1, which=3)
residuo <- residuals(modelo1)
cores <- c("green","darkgreen","blue","darkblue")
map.residuos <- SpatialPointsDataFrame(data=data.frame(residuo),
coords=cbind(MF$COORD_X,MF$COORD_Y))
spplot(map.residuos, cuts=quantile(residuo), col.regions=cores, cex=1)
GWRBandwith <- gwr.sel(MF$MFALC_15_2 ~
MF$MFEFF2015+MF$agro1315+MF$MEDALC+MF$density, data=MF,
coords=cbind(MF$COORD_X,MF$COORD_Y), adapt = TRUE)
gwr.model = gwr(MF$MFALC_15_2 ~
MF$MFEFF2015+MF$agro1315+MF$MEDALC+MF$density, data=MF,
coords=cbind(MF$COORD_X,MF$COORD_Y), adapt = GWRBandwith, hatmatrix = TRUE,
se.fit = TRUE)
gwr.model
resultadogeoregress <- as.data.frame(gwr.model$SDF)
head(resultadogeoregress)
MF$coefMFEFF2015 <- resultadogeoregress$MF.MFEFF2015
MF$coefagro1315 <- resultadogeoregress$MF.agro1315
MF$coefMEDALC <- resultadogeoregress$MF.MEDALC
MF$coefDENSITY <- resultadogeoregress$MF.density
workdir <- "C:\\Users\\pedro\\Documents\\Time Series\\Luciano"
cidades <- readShapePoly("pr2016.shp")
cidadesDF <- fortify(cidades, region = "nome")
head(cidadesDF)

ggplot(data=cidadesDF, aes(y=lat, x=long, group=group,
fill=id))+geom_polygon()+guides(fill=FALSE)
plot(cidades) #checking the .shp structure

#this is where things start to go wrong
paranacidade <- ggplot(MF,
aes(x=x,y=y))+geom_point(aes(colour=MF$coefMFEFF2015))+
scale_colour_gradient2(low="red",
mid="white",high="blue",midpoint=0,space="rgb",
na.value="grey50",guide="colourbar",guide_legend(title="coefs"))
paranacidade+geom_path(data=cidadesDF, aes(group=group),
colour="grey")+coord_equal()

It returns the error "Aesthetics must be either length 1 or the same as the data (399): colour, x, y"
The fortify function breaks the polygons of my shapefile into millions of rows, but it works for plotting the map by itself. When I try to plot the regression (MF$coefs...) within the shape the lenght of the regression (399) is not compatible with the lenght of rows in the shapefile dataframe. I've tried changing the arguments for aes (tried x=long, group=group, group=id,y=lat etc...)

The data:
shp+other: https://drive.google.com/open?id=1OAtXe9rkDUJVVqRSuO-V2Okb5b8x9aoO

csv: https://drive.google.com/open?id=1T1zm9BXyyFYJmqZphaRg72OWkCk1KpAh separated by commas


texts: ST_AsMVT vector tiles distort at low zoom levels

texts:
I'm serving up a data from PostGIS by converting it into a tile and then bytes and then onto a leaflet map. It all works, with one exception, the latitudes stretch ever greater towards the poles when zoomed out, but this effect disappears when zooming in. Both pink and blue lines are from the same data set.

When zoomed in the dataset come together/align suggesting there isn't a problem with the initial dataset.

When zooming out to levels of 5 and below the tiles produced seem to distort and stretch ever more towards the poles. The tiles near the equator are effected a lot less.

This is the code python that is used to produce the tiles from the server:

#

def tile_ul(x, y, z):
    n = 2.0 ** z
    lon_deg = x / n * 360.0 - 180.0
    lat_rad = arctan(sinh(pi * (1 - 2 * y / n)))
    lat_deg = degrees(lat_rad)
    return  lon_deg,lat_deg
############################################################
def get_tile(z,x,y):
    xmin, ymin = tile_ul(x, y, z)
    xmax, ymax = tile_ul(x + 1, y + 1, z)
    query = """SELECT ST_AsMVT(tile) FROM (
               SELECT id, ST_AsMVTGeom(geom, ST_MakeEnvelope( %s, %s, %s, %s ,4326), 4096, 256, false ) geom
               FROM seaarea ) as tile"""
    cursor = db.connection.cursor()
    cursor.execute(query,(xmin,ymin,xmax,ymax))
    tile = bytes(cursor.fetchone()[0])
    cursor.close()
    return tile

the client code which is used before adding the tile layer to a leaflet map:

var url = "http://localhost:5000/tiles/{z}/{x}/{y}"
var mapillaryLayer = L.vectorGrid.protobuf(url, {vectorTileLayerStyles: {
            fillColor: 'red'}}).addTo(mymap);

I believe the issue may lay in the PostGIS query or possibly the bounds created for the tiles when at low zoom levels. A similar issue is report here

Any suggestions on how to make the tiles not distort at low zoom levels?

Stay Safe All.


texts: Creating a popup when hovering a line in Mapbox

texts:
I have a custom line in mapbox that displays a popup when you click it, how can I make it so it shows only when you hover it ? I changed from "click" to "mouseenter" but it did not closed the popup. What needs changing so when I hover over this line it shows the popup and then it closes when leaving the line/popup?
Here is a jsfiddle example : https://jsfiddle.net/andrre/jb4czmew/4/
JS
mapboxgl.accessToken = 'pk.eyJ1IjoibWFya2V0aW5nYnNvIiwiYSI6ImNrYnYwZmk3YjAxZjgyem1wY2Zmc3F4Y2EifQ.gMF-eCCaAHHgWIUoRcnfkg';
    var map = new mapboxgl.Map({
        container: 'map',
        style: 'mapbox://styles/mapbox/streets-v11',
        center:  [-0.066985, 51.503363],
        zoom: 9
    });

map.on('load', function() {
        map.addSource('routeThree', {
            'type': 'geojson',
            'data': {
                'type': 'Feature',
                'properties': {},
                'geometry': {
                    'type': 'LineString',
                    'coordinates': [
                        [-0.066985, 51.503363],
                        [-3.550610, 40.390555 ],
                        [-77.140217, 38.801481],
                        [-118.175979, 34.008447]
                    ]
                }
            }
        });
        map.addLayer({
            'id': 'routeThree',
            'type': 'line',
            'source': 'routeThree',
            'layout': {
                'line-join': 'round',
                'line-cap': 'round'
            },
            'paint': {
                'line-color': '#04A5BD',
                'line-width': 8
            }
        });
   map.on('click', 'routeThree', function(e) {
            new mapboxgl.Popup()
            .setLngLat(e.lngLat)
            .setHTML("
This is the third line that will explain something
")
            .addTo(map);

         });

    });


texts: Append a placeholder value to a list wherever a value is not present in a table

texts:
Basically my tables have a field called EOS with four possible values. I dissolve my table on the EOS field, and what I need to do now is append the EOS values to a list. The tricky part is that I want to append a value to the list for each possible EOS value, whether or not the value exists in the table.

So if my table looks like this:

EOS       Shape Length
  1           56.7777
  3           9589.88

This is how I need the list to output:

EOS_list = ['56.777', **'0'**, 9589.88', **'0'**]

I can't figure out how to append a zero value to the list where a value does not exist without getting things tangled up.

workspace = "C:\MyDocuments\Documents\Update.gdb"
testing = r"C:\MyDocuments\Documents\Update.gdb\Dissolve"
arcpy.env.workspace = testing
arcpy.env.overwriteOutput = True

fields = ["EOS", "Shape_Length"]
EOS_list = []
cursor = arcpy.da.SearchCursor("Dissolve_EOS", fields)
arcpy.management.MakeFeatureLayer("Dissolve_EOS", "fc_EOS")
for f1, f2 in cursor:
    print f1
    print f2
    if arcpy.Exists(f1) and f1 == "1":
        EOS_list.append(f2)
        print EOS
    elif arcpy.Exists(f1) and f1 == "2":
        EOS_list.append(f2)
    elif arcpy.Exists(f1) and f1 == "3":
        EOS_list.append(f2)


texts:
I'm writing a simple script to export postgis data to geojson.
The ogr2ogr command works well, but I'm having problems when I'm passing the parameters.

#!/bin/bash
echo 'enter country code'
read export
echo 'exporting data from $export'
ogr2ogr -f GeoJSON export.geojson PG:'dbname=alamedaok user=postgres host=localhost' -sql "SELECT id_alameda, pais, geom FROM poblaciones_def WHERE pais = `$export`"

The problem comes in the WHERE pais = **$export"** I don't know how "insert" the variable in the sql sentence.
I've tried with '$' and with `, with no luck.

Any help?


texts:
This question may be a little vague, but I'm not entirely sure the right question to ask here...

I have a MrSID file that I want to display in NASA worldwind. To do this, would I need to convert the MrSID region of interest to a jpg (or something else that worldwind could display)?


texts: Adding WMS Service for Silverlight Viewer

texts: Is there a way to add a WMS image service in basemap gallery?


texts: Merging Rows in Table using ArcGIS Desktop?

texts:
I am fairly new to GIS and I am using ArcGIS Desktop 10.4.

The wildlife monitoring company I work for uses camera pods to sample fisheries population and uses side scan sonar to collect habitat data. Together this data is entered into ArcMap and analyzed for further study where Camera Sites = point and Habitat = Polygon.

What I want to do for each camera site (point) is create a buffer (50 meters). After creating a buffer I would like to calculate how much area (sq. meters) of what type of habitat fall within each buffer.

To do this I have:

Created a buffer around Camera sites
Performed an Intersect between Buffer and Habitat layers

Originally I had 184 records (184 camera sites). After I run the Intersect I have close to 700 records resulting in many rows for each camera site. This makes sense because there is often more than one polygon per habitat type within each buffer

Reference = Camera Site

Below is an example of the resulting attribute table where Reference = GMM120304501:

What I want is to have one line for each camera site showing area for each habitat. To do this I used the Pivot Table tool to make Geoforms column headings and the Area_50 field to fill the values.

The Pivot table result is below. My issue is that there are multiple lines for individual camera sites.

Example image of my issue where Reference is same as above ex. (GMM120304501).

Is there a way to Merge these values in the table above based on Reference and have the habitat types add together into one row?

Once this problem is resolved I plan on repeating the process with 3 buffers total and joining all the tables together based on the Reference Field.


texts: I cannot start up QGIS and get an error message

texts:
When I start up my QGIS, the following appears on the screen and it closes down. Please help me; I really need my files that are there and do my work as well.  The message says:

The ordinal 284 could not be located in the dynamic link library SSLEAY32.dll.


texts:
I try to add a point xy to feature class and I want to select the latest record in table (recently added). The script runs and I don't have any error message but the tool doesn't select any thing. How can I fix this?

import arcpy
arcpy.env.workspace = r'D:\Geodatabases\BD_VOCATION.gdb '
points = r'D:\Geodatabases\BD_VOCATION.gdb\Loc_vocation'
try:
    arcpy.MakeFeatureLayer_management(points,'points_layer')
    rows = arcpy.SearchCursor(points,'','',"OBJECTID", "OBJECTID D")
    lst = []
    for row in rows:
        m = row.getValue("OBJECTID")
        lst.append(m)
    v = lst[0]
    r= str(v)
    query = """{0} = {1}""".format(arcpy.AddFieldDelimiters(points, "OBJECTID"), int(r))
    arcpy.SelectLayerByAttribute_management("points_layer", "NEW_SELECTION", query)
except:
    print "An error occured during selection"


texts: Is it possible to override default QGIS error messages?

texts: From my QGIS plugin I would like to override/replace one of the the default QGIS error message ("The layer layer_url is not a valid layer and can not be added to the map") with my own message. I have looked at the API and can't find anything along these lines. Is it possible to achieve this?


texts:
I am writing a plugin that includes interaction with the map. Therefore, I use a QgsMapTool. It is possible for a user to mark a point on the map. When the user clicks a red X is drawn on the mapCanvas via QgsVertexMarker.
vertex_marker = QgsVertexMarker(self.canvas)
vertex_marker.setCenter(QgsPoint(map_coordinates['x'], map_coordinates['y']))
vertex_marker.setColor(QColor(255, 0, 0))
vertex_marker.setIconSize(7)
vertex_marker.setIconType(QgsVertexMarker.ICON_X)  # ICON_BOX, ICON_CROSS, ICON_X
vertex_marker.setPenWidth(2)

This is working fine. And the red X's are displayed on the map. But I can't delete the VertexMarker. I am trying to delete the VertexMarker from the scene. I've used something like this:
vertex_items = [ i for i in iface.mapCanvas().scene().items() if issubclass(type(i), qgis.gui.QgsVertexMarker)]

for ver in vertex_items:
    if ver in iface.mapCanvas().scene().items():
        iface.mapCanvas().scene().items().remove(ver)

iface.mapCanvas().refresh()

Like this I get the used VertexMarker that are visible on the scene. But the remove function somehow doesn't remove the marker from the scene.
Is there another possibility to remove the QgsVertexMarker from the scene?


texts: Why is the database connection lost on every query?

texts:
Every time I started a query the connection with the database is lost. Anyone familiar with this problem?

This problem occurs with various query, for instance:

CREATE TABLE table
SELECT a.*, b.column, b.column
FROM table_a a JOIN table b
ON ST_intersects (a.geom, b.the_geom)
WHERE st_IsValid (a.geom);

I just updated from postgis 1.5.3 to 1.5.5 but that doesn't seem to help.

Anyone an idea?

UPDATE

I solved the problem.

There was something wrong with the topology of one of the layers. After editing the query worked without disconnecting.

I find it rather strange that the program and service crashes when the topology is not correct. Can someone shed a light on this issue.


texts: How can I select a polygon in MS SQL 2008

texts:
I am very new to MS SQL 2008. I got a task to select a polygon from the polygon layer by giving the X-Y coordinates. I already got the script for PostgreSQL for the same thing, but I need it for MS SQL 2008.
the code in PostgreSQL is:

SELECT *
FROM "Kar_bbmp_198"
WHERE ST_Within
  ((ST_GeomFromText('POINT(773824 1442012)',32643)), the_geom) = TRUE;`

Can anybody help me from this?


texts:
I cross posted this in the Google Earth forum but they changed the format there and I'm not sure I like the format.

In my dataset I sometimes have multiple locations at a single point.

GE will represent this as a single dot until you click on the dot and then you get a scattering of all the points at that location.

Is there anyway to script into the KML a sizing on the initial dot viewed to represent how many dots are on that location? (ex. 1-3 small dot, 4-6 med dot, 7-9 large dot)?

Thanks


texts:
I've created a VRT from some 6-band TIFF tiles. The VRT displays perfectly in QGIS 3.6 and I can see manipulate the band arrangement and create different composites.

However, when I open the raster calculator, it only shows one band available. If I write the calculation anyway with the band numbers ("image@2" + "image@3" ..). it does not work.